Структурные указатели в Cython с минимизацией GSL Monte-Carlo
Я застрял в этом упражнении и не достаточно хорош, чтобы решить его. В основном я пишу алгоритм максимального правдоподобия Монте-Карло для распределения Бернулли. Проблема заключается в том, что мне нужно передать данные в качестве параметра алгоритму минимизации GSL (однонаправленного), а также передать размер данных (поскольку во внешнем цикле используются разные размеры выборки "наблюдаемых" данных).). Поэтому я пытаюсь передать эти параметры в виде структуры. Тем не менее, я сталкиваюсь с ошибками сегмента и уверен, что это происходит из той части кода, которая касается структуры, и рассматривает ее как указатель.
[РЕДАКТИРОВАТЬ: я исправил для распределения структуры и ее компонентов]
%%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)
Я отдельно протестировал распределение структуры, и она работает, делая то, что должна делать. Однако во время работы с Монте-Карло ядро убивается вызовом прерывания (согласно выводу на моем Mac), и вывод Jupyter на моей консоли выглядит следующим образом:
gsl: fsolver.c:39: ERROR: computed function value is infinite or NaN
Вызван обработчик ошибок GSL по умолчанию.
Теперь кажется, что решатель не работает. Я не знаком с пакетом GSL, так как использовал его только один раз для генерации случайных чисел из распределения gumbel (минуя команды scipy).
Буду признателен за любую помощь в этом! Спасибо
[РЕДАКТИРОВАТЬ: Изменить нижнюю границу]
Повторяя упражнение с экспоненциальным распределением, функция логарифмического правдоподобия которого содержит только один журнал, я справился с проблемой, с которой был gsl_min_fminimizer_set Первоначальная оценка по нижней границе a в 0 дает результат -INF (так как он оценивает проблему до решения, чтобы сгенерировать f(нижний уровень), f(верхний уровень), где f - моя функция для оптимизации). Когда я устанавливаю нижнюю границу к чему-то, кроме 0, но очень маленькому tol переменная моего определенного допуска) алгоритм решения работает и дает правильные результаты.
Большое спасибо @DavidW за подсказки, которые помогут мне добраться туда, куда мне нужно было пойти.
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с, чтобы найти A, Нахождение места, где производная равна 0, дает p=A, (Вы захотите перепроверить мою работу!). Благодаря этому вы можете избежать использования процедур минимизации GSL.