Hi,
Thank you for his response. I attach the program in c for details.
El 17/06/12 10:12, Juan Pablo Amorocho D. escribió:
Hola Claudio,
Without the code and concrete problem is very difficult to say
anything about your question.
Cheers,
Juan Pablo
2012/6/17 Claudio Tablada <[email protected]
<mailto:[email protected]>>
Hi,
I am trying to maximize a -log likelihood function of a
Birnbaum-Saunders distribution with two parameters using BFGS and
Monte Carlo. I have the samples generates using
gsl_ran_gaussian_ziggurat and later transforming for the
Birnbaum-Saunders distribution. The problem is that for some
values of the sample the algorithm converges and for other values
it results in NaN. Thanks in advanced.
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_math.h>
#include <math.h>
#define N 100 //tamaño de la muestra
#define B 1.0 //parámetro beta fijado en 1.0
#define NREP 200 //número de replicas Monte Carlo
double logv(const gsl_vector *v , void *params);
void dlogv(const gsl_vector *v , void *params , gsl_vector *df);
void logv_dlogv(const gsl_vector *x , void *params , double *f , gsl_vector *df);
double *mle(double array[]);
double *montecarlo(double *mle(double array1[]) , double alpha);
gsl_rng * r;
gsl_multimin_fdfminimizer *s;
gsl_vector *x;
double z[4];
int main(void)
{
int i , op;
float alpha;
double sigma;
double *a;
const gsl_rng_type *T;
FILE *ofp1; //puntero para archivo
printf("Ingrese un valor del parámetro alpha\n");
scanf("%f",&alpha);
sigma = alpha/2;
ofp1 = fopen("saida1.txt","a");
/*Crea un generador por la variable GSL_RNG_TYPE*/
gsl_rng_env_setup();
T = gsl_rng_default;
r = gsl_rng_alloc (T);
/*Generación de las replicas Monte Carlo*/
a = montecarlo(mle , alpha);
fprintf(ofp1, "Estimación de los parámetros alpha y beta de máxima verosimilitud\n");
fprintf(ofp1 , "Número de replicas Monte Carlo: %d\n" , NREP);
fprintf(ofp1 , "Salida para muestra con N=%i, alpha=%1.2f e beta=1.00\n" , N , alpha);
fprintf(ofp1 , "alpha_MC: % 2.6f\t beta_MC: % 2.6f\n", a[0] , a[1]);
fprintf(ofp1 , "sesgo relativo alpha: % 2.6f\t sesgo relativo beta: % 2.6f\n", (a[0] - alpha)/alpha , (a[1] - B)/B);
fprintf(ofp1 , "EQM alpha: % 2.6f\t EQM beta: % 2.6f\n", (a[0] - alpha)*(a[0] - alpha) + (a[2] - a[0]*a[0]) , (a[1] - B)*(a[1] - B) + (a[3] - a[1]*a[1]));
fprintf(ofp1 , "Var alpha: % 2.6f\t Var beta: % 2.6f\n\n" , a[2] - a[0]*a[0] , a[3] - a[1]*a[1]);
return 0;
}
/*-log-likelihood function*/
double logv(const gsl_vector *v , void *params)
{
int i;
double x1 , y1 , soma1 , soma2;
double *p = (double *)params;
x1 = gsl_vector_get(v , 0);
y1 = gsl_vector_get(v , 1);
for(i = 0 , soma1 = 0.0 ; i < N ; i++)
soma1 += log(sqrt(y1/p[i]) + sqrt((y1/p[i])*(y1/p[i])*(y1/p[i])));
for(i = 0 , soma2 = 0.0 ; i < N ; i++)
soma2 += p[i]/y1 + y1/p[i] - 2;
return N*log(x1*y1) - soma1 + soma2/(2*x1*x1);
}
/*gradient of logv, dlogv=(dlogv/dx,dlogv/dy)*/
void dlogv(const gsl_vector *v , void *params , gsl_vector *df)
{
int i;
double x1 , y1 , soma1 , soma2 , soma3;
double *p = (double *)params;
x1 = gsl_vector_get(v , 0);
y1 = gsl_vector_get(v , 1);
for(i = 0 , soma1 = 0.0 ; i < N ; i++)
soma1 += p[i];
for(i = 0 , soma2 = 0.0 ; i < N ; i++)
soma2 += 1./p[i];
for(i = 0 , soma3 = 0.0 ; i < N ; i++)
{
//soma3 = 1./(p[i] + y1);
soma3 += (1./(sqrt(y1/p[i]) + sqrt((y1/p[i])*(y1/p[i])*(y1/p[i]))))*(1./(2*p[i]))*(sqrt(p[i]/y1) + 3*sqrt(y1/p[i]));
}
gsl_vector_set(df , 0 , (N/x1)*(1. + 2./(x1*x1)) - soma1/(x1*x1*x1*y1) - (y1/(x1*x1*x1))*soma2);
gsl_vector_set(df , 1 , N/(y1) - soma3 - soma1/(2*x1*x1*y1*y1) + soma2/(2*x1*x1));
}
/*logv and dlogv*/
void logv_dlogv(const gsl_vector *x , void *params , double *f , gsl_vector *df)
{
*f = logv(x , params);
dlogv(x , params , df);
}
/*mle function (maximum-likelihood estimators)*/
double *mle(double array[])
{
const gsl_multimin_fdfminimizer_type *V;
gsl_multimin_function_fdf func;
int iter = 0;
int i , status;
func.n = 2; //número de componentes
func.f = &logv;
func.df = &dlogv;
func.fdf = &logv_dlogv;
func.params = array;
/*Punto inicial x=(1,1)*/
x = gsl_vector_alloc (2);
gsl_vector_set(x , 0 , 1.0);
gsl_vector_set(x , 1 , 1.0);
V = gsl_multimin_fdfminimizer_vector_bfgs2;
s = gsl_multimin_fdfminimizer_alloc(V , 2);
gsl_multimin_fdfminimizer_set(s , &func , x , 0.01 , 1e-4);
do{
iter++;
status = gsl_multimin_fdfminimizer_iterate(s);
if(status)
break;
status = gsl_multimin_test_gradient(s -> gradient , 1e-3);
if(status == GSL_SUCCESS)
printf("Minimo hallado en: \n");
printf("%5d %.5f %.5f %10.5f\n" , iter , gsl_vector_get(s->x , 0) , gsl_vector_get(s->x , 1) , s->f);
}
while(status == GSL_CONTINUE && iter < 100);
z[0] = gsl_vector_get(s->x , 0);
z[1] = gsl_vector_get(s->x , 1);
return z;
}
/*montecarlo function*/
double *montecarlo(double *mle(double array1[]) , double alpha)
{
int i , k;
double u , soma1 , soma2 , soma3 , soma4 , sigma;
double vetor_MC[N];
double *a;
sigma = alpha/2;
for(k = 0 , soma1 = 0.0 , soma2 = 0.0 , soma3 = 0.0 , soma4 = 0.0 ; k < NREP ; k++)
{
for(i = 0 ; i < N ; i++)
{
double u = gsl_ran_gaussian_ziggurat (r , sigma);
/*generación de los datos con distribución Birnbaum-Saunders*/
vetor_MC[i] = (double)B*(1 + 2*u*u + 2*u*sqrt(1 + u*u));
}
a = mle(vetor_MC);
soma1 += a[0];
soma2 += a[1];
soma3 += a[0]*a[0];
soma4 += a[1]*a[1];
}
gsl_rng_free (r);
gsl_multimin_fdfminimizer_free(s);
gsl_vector_free(x);
z[0] = soma1/NREP;
z[1] = soma2/NREP;
z[2] = soma3/NREP;
z[3] = soma4/NREP;
return z;
}
