Cython 中的结构指针与 GSL Monte-Carlo 最小化

Question

我被困在这个练习上，还不足以解决它。基本上我正在为伯努利分布编写蒙特卡罗最大似然算法。问题是我必须将数据作为参数传递给 GSL 最小化（one-dim）算法，并且还需要传递数据的大小（因为外循环是“观察到”数据的不同样本大小）。所以我试图将这些参数作为结构传递。但是，我遇到了段错误，我确信它来自与结构相关的代码部分并将其视为指针。

[编辑：我已经更正了结构及其组件的分配]

%%cython

#!python
#cython: boundscheck=False, wraparound=False, nonecheck=False, cdivision=True   

from libc.stdlib cimport rand, RAND_MAX, calloc, malloc, realloc, free, abort
from libc.math cimport log

#Use the CythonGSL package to get the low-level routines
from cython_gsl cimport *

######################### Define the Data Structure ############################

cdef struct Parameters:
    #Pointer for Y data array
    double* Y
    #size of the array
    int* Size

################ Support Functions for Monte-Carlo Function ##################

#Create a function that allocates the memory and verifies integrity
cdef void alloc_struct(Parameters* data, int N, unsigned int flag) nogil:

    #allocate the data array initially
    if flag==1:
        data.Y = <double*> malloc(N * sizeof(double))
    #reallocate the data array
    else:
        data.Y = <double*> realloc(data.Y, N * sizeof(double))

    #If the elements of the struct are not properly allocated, destory it and return null
    if N!=0 and data.Y==NULL:
        destroy_struct(data)
        data = NULL     

#Create the destructor of the struct to return memory to system
cdef void destroy_struct(Parameters* data) nogil:
    free(data.Y)
    free(data)

#This function fills in the Y observed variable with discreet 0/1
cdef void Y_fill(Parameters* data, double p_true, int* N) nogil:

    cdef:
        Py_ssize_t i
        double y

    for i in range(N[0]):

        y = rand()/<double>RAND_MAX

        if y <= p_true:
            data.Y[i] = 1 
        else:
            data.Y[i] = 0
#Definition of the function to be maximized: LLF of Bernoulli
cdef double LLF(double p, void* data) nogil:

    cdef:
        #the sample structure (considered the parameter here)
        Parameters* sample

        #the total of the LLF
        double Sum = 0

        #the loop iterator
        Py_ssize_t i, n

    sample = <Parameters*> data

    n = sample.Size[0]

    for i in range(n):

        Sum += sample.Y[i]*log(p) + (1-sample.Y[i])*log(1-p)

    return (-(Sum/n))

########################## Monte-Carlo Function ##############################

def Monte_Carlo(int[::1] Samples, double[:,::1] p_hat, 
                Py_ssize_t Sims, double p_true):

    #Define variables and pointers
    cdef:
        #Data Structure
        Parameters* Data

        #iterators
        Py_ssize_t i, j
        int status, GSL_CONTINUE, Iter = 0, max_Iter = 100 

        #Variables
        int N = Samples.shape[0] 
        double start_val, a, b, tol = 1e-6

        #GSL objects and pointer
        const gsl_min_fminimizer_type* T
        gsl_min_fminimizer* s
        gsl_function F

    #Set the GSL function
    F.function = &LLF

    #Allocate the minimization routine
    T = gsl_min_fminimizer_brent
    s = gsl_min_fminimizer_alloc(T)

    #allocate the struct
    Data = <Parameters*> malloc(sizeof(Parameters))

    #verify memory integrity
    if Data==NULL: abort()

    #set the starting value
    start_val = rand()/<double>RAND_MAX

    try:

        for i in range(N):

            if i==0:
                #allocate memory to the data array
                alloc_struct(Data, Samples[i], 1)
            else:
                #reallocate the data array in the struct if 
                #we are past the first run of outer loop
                alloc_struct(Data, Samples[i], 2)

            #verify memory integrity
            if Data==NULL: abort()

            #pass the data size into the struct
            Data.Size = &Samples[i]

            for j in range(Sims):

                #fill in the struct
                Y_fill(Data, p_true, Data.Size)

                #set the parameters for the GSL function (the samples)
                F.params = <void*> Data
                a = tol
                b = 1

                #set the minimizer
                gsl_min_fminimizer_set(s, &F, start_val, a, b)

                #initialize conditions
                GSL_CONTINUE = -2
                status = -2

                while (status == GSL_CONTINUE and Iter < max_Iter):

                    Iter += 1
                    status = gsl_min_fminimizer_iterate(s)

                    start_val = gsl_min_fminimizer_x_minimum(s)
                    a = gsl_min_fminimizer_x_lower(s)
                    b = gsl_min_fminimizer_x_upper(s)

                    status = gsl_min_test_interval(a, b, tol, 0.0)

                    if (status == GSL_SUCCESS):
                        print ("Converged:\n")
                        p_hat[i,j] = start_val

    finally:
        destroy_struct(Data)
        gsl_min_fminimizer_free(s)

用下面的python代码运行上面的函数：

import numpy as np

#Sample Sizes
N = np.array([5,50,500,5000], dtype='i')

#Parameters for MC
T = 1000
p_true = 0.2

#Array of the outputs from the MC
p_hat = np.empty((N.size,T), dtype='d')
p_hat.fill(np.nan)

Monte_Carlo(N, p_hat, T, p_true)

我已经单独测试了结构分配并且它工作正常，做它应该做的。然而，在玩 Monte Carlo 时，内核被一个 abort 调用杀死（根据我 Mac 上的输出），我的控制台上的 Jupyter 输出如下：

gsl: fsolver.c:39: ERROR: computed function value is infinite or NaN

调用了默认的 GSL 错误处理程序。

现在看来求解器不起作用。我不熟悉 GSL 包，只使用过一次从 gumbel 分布中生成随机数（绕过 scipy 命令）。

我将不胜感激任何帮助！谢谢

[编辑：更改 a 的下限]

用指数分布重做练习，它的对数似然函数只包含一个对数我已经磨练了这个问题，gsl_min_fminimizer_set 最初在 a 的下限评估为 0 产生 -INF 结果（因为它在求解以生成 f(lower), f(upper) 之前评估问题，其中 f 是我要优化的函数）。当我将下限设置为 0 以外但非常小的值（比如我定义的容差的 tol 变量）时，解决方案算法起作用并产生正确的结果。

非常感谢@DavidW 提供的提示，让我能够到达我需要去的地方。

Answer 1

这是一个有点投机的答案，因为我没有安装 GSL，所以很难测试它（如果有错误，请道歉！）

---

我认为问题是线

Sum += sample.Y[i]*log(p) + (1-sample.Y[i])*log(1-p)

看起来Y[i] 可以是 0 或 1。当p 位于 0-1 范围的任一端时，它会给出0*-inf = nan。在只有所有“Y”相同的情况下，该点是最小值（因此求解器将可靠地最终到达无效点）。幸运的是，您应该能够重写该行以避免获得nan：

if sample.Y[i]:
   Sum += log(p)
else:
   Sum += log(1-p)

（将生成nan 的情况是未执行的情况）。

---

我发现了第二个小问题：在alloc_struct 中，如果出现错误，您可以使用data = NULL。这只会影响本地指针，因此您在Monte_Carlo 中对NULL 的测试毫无意义。你最好从alloc_struct 返回一个真或假标志并检查它。我怀疑你是否遇到了这个错误。

---

编辑：另一个更好的选择是通过分析找到最小值：A log(p) + (1-A) log (1-p) 的导数是A/p - (1-A)/(1-p)。平均所有sample.Ys 以​​找到A。找到导数为 0 的位置给出p=A。 （您需要仔细检查我的工作！）。这样您就可以避免使用 GSL 最小化例程。