用GSL蒙特卡洛最小化的Cython中的結構指針

Question

我被困在這個練習中，我不夠好解決它。基本上我正在爲伯努利分佈寫一個蒙特卡洛最大似然算法。問題在於我必須將數據作爲參數傳遞給GSL最小化（one-dim）算法，並且還需要傳遞數據的大小（因爲外部循環是「觀察」數據的不同樣本大小）。所以我試圖通過這些參數作爲一個結構。但是，我遇到了seg錯誤，我確定它來自涉及結構的代碼部分，並將其視爲指針。

[編輯：我已校正的結構的配置和它的組成]

%%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) 

我已經單獨測試結構分配和它的工作，做它應該做的。然而，雖然funning蒙特卡洛內核與中止呼叫（每輸出在我的Mac）死亡，我的控制檯上的Jupyter輸出如下：

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

默認GSL錯誤處理程序調用。

現在好像解算器不工作了。我對GSL包並不熟悉，只用它一次從gumbel發佈生成隨機數（繞過scipy命令）。

我將不勝感激任何幫助！由於

[編輯：修改下界的]

重做練習與指數分佈，其數似然函數只包含一個日誌我已經磨練下已經與最初gsl_min_fminimizer_set評估的問題一個0的下界產生-INF結果（因爲它在求解產生f（較低），f（較高）之前評估問題，其中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。這隻會影響本地指針，因此您對NULL的測試Monte_Carlo沒有意義。你最好從alloc_struct返回一個真或假的國旗並檢查它。我懷疑你是否遇到這個錯誤。

---

編輯：另一種更好的選擇是分析找出最小：中A log(p) + (1-A) log (1-p)的衍生物爲A/p - (1-A)/(1-p)。平均所有的sample.Y s找到A。找到導數爲0的地方給出p=A。 （您需要仔細檢查我的工作！）。有了這個，你可以避免使用GSL最小化例程。