Увеличить размер gsl minimizer

Я должен разработать оптимизатор для расчета 14 параметров скрытой цепочки Маркова (я знаю, это много, но возможно :)). Я уже сделал тот же оптимизатор, используя IMSL (который у меня был во время стажировки), и я знаю, что это возможно, но у меня сейчас нет лицензии на IMSL, потому что я делаю студенческий проект и пытаюсь сделать тот же оптимизатор в C ++.

Но я получил эту ошибку (я прокомментировал в коде, где он появляется):

gsl: simplex.c: 288: ОШИБКА: встречается не конечное значение функции

Вызван обработчик ошибок GSL по умолчанию.

Отменено

Я использую GSL libraryболее конкретно, суббиблиотека multimin,

Я показываю вам мой код здесь и после, и я буду рад комментировать его, если вам нужно, но я могу отправить его с данными любому, кто сможет мне помочь.

Большое спасибо за ваше время и помощь.

Benjamin,

Структура нам нужна для оптимизатора:

    struct dataMom {
int         n;    /* Number of data points */
gsl_vector *Rmom;    /* Response data */
gsl_vector *Rm;    /* Covariate */
};

Функция, которая вычисляет значения

double my_f (const gsl_vector *v, void *params) {
struct dataMom *p = (struct dataMom *) params;
int i;

double loglik = 0;

Итак, вы можете видеть, что у нас есть 14 параметров для оптимизации.

    double alphaC = gsl_vector_get (v, 0);
double beta0C =gsl_vector_get (v, 1);
double betaUC =gsl_vector_get (v, 2);
double sigmamomC =gsl_vector_get (v, 3);
double muC =gsl_vector_get (v, 4);
double sigmamC =gsl_vector_get (v, 5);
double probaC =gsl_vector_get (v, 6);

double alphaT= gsl_vector_get (v, 7);
double beta0T =  gsl_vector_get (v, 8);
double betaUT= gsl_vector_get (v, 9);
double sigmamomT= gsl_vector_get (v, 10);
double muT =gsl_vector_get (v, 11);
double sigmamT = gsl_vector_get (v, 12);
double probaT =gsl_vector_get (v, 13);

double epsilonmomC;
double epsilonmC;
double epsilonmomT;
double epsilonmT;
double Rmomentum;
double Rmarket;
double IUT;
double Pc;
double Pt;
double PC;
double PT;
double PI = 3.14159;

double probac = (1-probaT)/(2-probaC-probaT);
double probat = 1-probac;

for (i=0;i<p->n;i++) {

Rmomentum = gsl_vector_get(p->Rmom, i)*0.01;
Rmarket = gsl_vector_get(p->Rm, i)*0.01;
if (Rmarket>0){IUT = 1;}
else {IUT = 0;}

epsilonmomC = (1/sigmamomC)*(Rmomentum - alphaC - (beta0C + IUT*betaUC)*Rmarket);
epsilonmC = (1/sigmamC)*(Rmarket-muC);epsilonmomT = (1/sigmamomT)*(Rmomentum - alphaT - (beta0T + IUT*betaUT)*Rmarket);
epsilonmT = (1/sigmamT)*(Rmarket-muT);Pc = ((1/(sigmamomC*sqrt(2*PI)))*exp(-(epsilonmomC*epsilonmomC)/2)) * ((1/(sigmamC*sqrt(2*PI)))*exp(-(epsilonmC*epsilonmC)/2));
Pt = ((1/(sigmamomT*sqrt(2*PI)))*exp(-(epsilonmomT*epsilonmomT)/2)) * ((1/(sigmamT*sqrt(2*PI)))*exp(-(epsilonmT*epsilonmT)/2));

PC = Pc * (probaC*probac + (1-probaC)*(probat));
PT = Pt * ((1-probaT)*(probac) + probaT*probat);

loglik -=log(PC+PT);

probac = PC / (PC+PT);
probat = PT / (PC+PT);

}
return loglik;
}

Главный:

int main (void) {
/* Initialise data and model parameters */
//these are 2 vectors of return from txt files that i can give you
vector<double> Benchmark;
vector<double> Momentum;

//these are the 2 fonctions i'm using to get the return for the txt files
Benchmark = readBenchmark();
Momentum = readMomentum();

int i, n=Benchmark.size(), nparas=14;
double initial_p[14] = {0.0204,0.41,-0.52,0.0432 , 0.0098 , 0.0362 , 0.97, 0.0402 , -0.26 , -1.28 , 0.1105 , -0.0070 , 0.0904 , 0.92};
struct dataMom myStruct;
myStruct.n = n;
myStruct.Rmom = gsl_vector_alloc(myStruct.n);
myStruct.Rm = gsl_vector_alloc(myStruct.n);

vector<double>::iterator itbenchmark;
vector<double>::iterator itmomentum;
double tempmom;
double tempm;
itbenchmark = Benchmark.begin();
itmomentum = Momentum.begin();

for (i=0;i<myStruct.n;i++) {

tempmom = *itmomentum;
tempm = *itbenchmark;
gsl_vector_set (myStruct.Rmom, i, tempmom);
gsl_vector_set (myStruct.Rm, i, tempm);
itbenchmark++;
itmomentum++;
}

/* Starting point */
gsl_vector *paras;
paras = gsl_vector_alloc (nparas);
for (i=0;i<nparas;i++)
gsl_vector_set (paras, i, initial_p[i]);

/* Initialise algorithm parameters */
size_t iter = 0;
int status;
double size;

/* Set initial STEP SIZES to 1 */
gsl_vector *ss;
ss = gsl_vector_alloc (nparas);
gsl_vector_set_all (ss, 0.5);

const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
gsl_multimin_fminimizer *s = NULL;
gsl_multimin_function minex_func;

/* Initialize method and iterate */
minex_func.n = nparas;
minex_func.f = my_f;
minex_func.params = &myStruct;
s = gsl_multimin_fminimizer_alloc (T, nparas);

Пока здесь нет проблем, но следующая строка генерирует эту ошибку:

gsl: simplex.c: 288: ОШИБКА: встречается не конечное значение функции

Вызван обработчик ошибок GSL по умолчанию.

Прервано (ядро сброшено)

    gsl_multimin_fminimizer_set (s, &minex_func, paras, ss);do {
iter++;
status = gsl_multimin_fminimizer_iterate(s);

if (status)
break;

size = gsl_multimin_fminimizer_size (s);
status = gsl_multimin_test_size (size, 1);

if (status == GSL_SUCCESS)
printf ("converged to minimum at\n");

for (int z = 0 ; z<14 ; z++)
{
cout <<gsl_vector_get (s->x, z)<<endl;
}
}

while (status == GSL_CONTINUE && iter < 10000);

gsl_multimin_fminimizer_free (s);
gsl_vector_free(myStruct.Rmom);
gsl_vector_free(myStruct.Rm);
gsl_vector_free(paras);
gsl_vector_free(ss);

return status;
}