Code:
/*
* likelihood.c
*
*/
#include <unistd.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_sf_erf.h>
#include <gsl/gsl_sf_pow_int.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_sort_vector.h>
#include <gsl/gsl_linalg.h>
#include <ctime>
#include<gsl/gsl_multimin.h>
#include<gsl/gsl_deriv.h>
struct ML_params {
gsl_matrix *vy;
gsl_matrix *vyinit;
}parameters;
gsl_rng * r;
gsl_vector* sumr(gsl_matrix *mX )
{
size_t M=mX -> size1;
size_t N=mX -> size2;
int m, n;
gsl_vector *sum = gsl_vector_calloc(M);
double aux;
for(m=0;m<M;m++)
{
aux=0.0;
for(n=0;n<N;n++)
{
aux +=gsl_matrix_get(mX,m,n);
}
gsl_vector_set(sum,m,aux);
}
return sum;
}
gsl_matrix* compute_mD(const gsl_matrix*my,const gsl_matrix*myinit,const gsl_matrix *mbeta, const gsl_vector *vsigma, int p,int K)
/*
input: data matrix, location and scale parameters and size of data
output:f(y_t|y_t-1,theta)
*/
{
size_t T=my -> size1;
size_t cy=my -> size2;
int i,j;
int P=p;
int k;
int t,ctemp=T+p-1;
gsl_matrix *mD = gsl_matrix_calloc(T,K);
gsl_matrix *mX = gsl_matrix_calloc(T,P);
gsl_matrix *mytemp = gsl_matrix_calloc(ctemp,cy);
gsl_vector *vx = gsl_vector_calloc(P);
gsl_vector *vbeta = gsl_vector_calloc(P);
double aux,xbeta=0;
for(i=0;i<ctemp;i++)
{
for(j=0;j<cy;j++)
{
if(i<p-1)
gsl_matrix_set(mytemp,i,j,gsl_matrix_get(myinit,i, j));
else
gsl_matrix_set(mytemp,i,j,gsl_matrix_get(my,i-p+1,j));
}
}
for(t=0;t<T;t++)
gsl_matrix_set(mX,t,0,1);
for(t=0;t<T;t++)
{
for(i=1;i<P;i++)
{
gsl_matrix_set(mX,t,i,gsl_matrix_get(mytemp,(P-i-1+t),0));
}
}
for(t=0;t<T;t++)
{
for(k=0;k<K;k++)
{
gsl_matrix_get_row(vx, mX, t);
gsl_matrix_get_col(vbeta, mbeta, k);
xbeta=0.0;
for(i=0;i<P;i++)
xbeta+=gsl_vector_get(vx,i)*gsl_vector_get(vbeta,i );
aux =(exp(-0.5*((gsl_pow_int ((gsl_matrix_get(my,t,0)-xbeta),2)*1.0/gsl_vector_get(vsigma,k)))))*pow((2.0*M_PI*gsl_vec tor_get(vsigma, k)),-0.5);
gsl_matrix_set(mD, t, k, aux);
}
}
// deallocate memory
gsl_matrix_free (mX);
gsl_matrix_free (mytemp);
gsl_vector_free (vx);
gsl_vector_free (vbeta);
return mD;
}
gsl_matrix * compute_mF(const gsl_matrix *vy,
const gsl_matrix *vyinit,
const gsl_matrix *mbeta,
const gsl_vector *vsigma,
const gsl_matrix *mP,
int p,
int K)
{
int t,i,k;
size_t T=vy -> size1;
gsl_matrix *mF = gsl_matrix_calloc(T,K);
gsl_matrix *mD = gsl_matrix_calloc(T,K);
gsl_matrix_set(mF, 0, 0, 1);
gsl_matrix_set(mF, T-1,K-1, 1);
mD=compute_mD(vy,vyinit,mbeta,vsigma,p,K);
gsl_vector *vtemp = gsl_vector_calloc(K);
gsl_vector *vF = gsl_vector_calloc(K);
gsl_vector *vD = gsl_vector_calloc(K);
double sum,aux;
for( t=1;t<T-1;t++)
{
gsl_matrix_get_row (vF, mF, t-1);
gsl_matrix_get_row (vD, mD, t);
for(k=0;k<K;k++)
{
aux=0.00;
for(i=0;i<K;i++)
aux+= gsl_matrix_get(mF,t-1,i) * gsl_matrix_get(mP,i,k);
gsl_vector_set(vtemp,k,aux);
}
gsl_vector_mul (vtemp, vD);//vtemp=vtemp.*mD[t][]
sum =0.0;
for(k=0;k<K;k++)
sum+= gsl_vector_get(vtemp,k);
gsl_vector_scale (vtemp, 1.0/sum);
for(k=0;k<K;k++)
gsl_matrix_set(mF, t, k,gsl_vector_get(vtemp,k));
}
// deallocate memory
gsl_matrix_free (mD);
gsl_vector_free (vtemp);
gsl_vector_free (vF);
gsl_vector_free (vD);
return mF;
}
gsl_matrix* compute_mF_likelihood(const gsl_matrix *vy,const gsl_matrix *vyinit,const gsl_matrix *mbeta, const gsl_vector *vsigma, const gsl_matrix *mP,int p, int K)
{
size_t T=vy -> size1;
gsl_matrix *mF=gsl_matrix_calloc(T,K);
mF=compute_mF(vy,vyinit,mbeta, vsigma, mP, p, K);
double aux;
gsl_vector *vtemp=gsl_vector_alloc(K);
gsl_matrix *QF=gsl_matrix_calloc(T,K);
gsl_matrix_set(QF,0,0,1);
gsl_matrix_set(QF,T-1,K-1,1);
for(int t=1;t<T-1;++t)
{
for(int k=0;k<K;k++)
{
aux=0.00;
for(int i=0;i<K;i++)
aux+= gsl_matrix_get(mF,t-1,i) * gsl_matrix_get(mP,i,k);
gsl_vector_set(vtemp,k,aux);
}
for(int k=0;k<K;k++)
gsl_matrix_set(QF, t, k,gsl_vector_get(vtemp,k));
}
gsl_vector_free(vtemp);
gsl_matrix_free(mF);
return QF;
}
double llikelihood_ML(const gsl_vector *vparamML,const gsl_matrix *vy,const gsl_matrix * vyinit,int p,int K)
{
size_t T=vy -> size1;
gsl_vector *vbetastar=gsl_vector_alloc(p*K);
for(int i=0;i<p*K;i++)gsl_vector_set(vbetastar,i,gsl_vecto r_get(vparamML,i));;
gsl_vector *vsigmastar_inv=gsl_vector_alloc(K);
for(int k=0;k<K;++k)
gsl_vector_set(vsigmastar_inv,k,pow(gsl_vector_get (vparamML,k+p*K),-1));
gsl_vector *vprobstar=gsl_vector_alloc(K-1);
for(int i=0;i<K-1;i++)gsl_vector_set(vprobstar,i,gsl_vector_get(vp aramML,i+p*K+K));;
gsl_matrix *mbetastar=gsl_matrix_calloc(p,K);
gsl_vector *vtemp=gsl_vector_calloc(p);
for(int k=0;k<K;++k)
{
for(int i=0;i<p;i++) gsl_vector_set(vtemp,i,gsl_vector_get(vbetastar,i+ k*p));
gsl_matrix_set_col (mbetastar,k,vtemp);
}
gsl_matrix *mP=gsl_matrix_calloc (K,K);
gsl_matrix_set(mP,K-1,K-1,1);
for(int k=0;k<K-1;++k)
{
for(int i=0; i<K-1;++i)
{
if(i==k)
{
gsl_matrix_set(mP,k,i,gsl_vector_get(vprobstar,k)) ;
gsl_matrix_set(mP,k,i+1,1-gsl_vector_get(vprobstar,k));
}
}
}
gsl_matrix *mlik=gsl_matrix_calloc(T,K);
mlik=compute_mD(vy,vyinit,mbetastar, vsigmastar_inv, p, K);
gsl_matrix *mprob=gsl_matrix_calloc(T,K);
mprob=compute_mF_likelihood(vy,vyinit,mbetastar,vs igmastar_inv, mP, p,K);
gsl_matrix *mtemp=gsl_matrix_calloc(T,K);
for(int t=0;t<T;t++)
{
for(int k=0;k<K;k++)
gsl_matrix_set(mtemp,t,k,gsl_matrix_get(mlik,t,k)* gsl_matrix_get(mprob,t,k));
}
gsl_vector *vllikelihood= gsl_vector_alloc(T);
vllikelihood= sumr(mtemp);
for(int t=0;t<T;t++)
gsl_vector_set(vllikelihood,t,log(gsl_vector_get(v llikelihood,t)));
double llik=0.0;
for(int t=0;t<T;++t)
llik+=gsl_vector_get(vllikelihood,t);
gsl_vector_free(vbetastar);
gsl_vector_free(vsigmastar_inv);
gsl_vector_free(vprobstar);
gsl_matrix_free(mbetastar);
gsl_vector_free(vtemp);
gsl_matrix_free(mP);
gsl_matrix_free(mlik);
gsl_matrix_free(mprob);
gsl_matrix_free(mtemp);
gsl_vector_free(vllikelihood);
return llik;
}
/* Compute Gradient by central differences Maheu*/
/* For use with GSL minimization routines */
void gsl_central_diff (const gsl_multimin_function * f, const gsl_vector * x, gsl_vector * g)
{
size_t i, n = f->n;
double h = GSL_SQRT_DBL_EPSILON;
// double h = pow(GSL_DBL_EPSILON,1.0/3.0);
double tmp,dx;
gsl_vector * x1 = gsl_vector_alloc (n);
gsl_vector_memcpy (x1, x);
for (i = 0; i < n; i++)
{
double fl, fh;
double xi = gsl_vector_get (x, i);
// double dx = fabs(xi) * h;
tmp = GSL_MAX_DBL ( fabs(xi)*h , h ) + xi;
dx = tmp - xi;
if (dx == 0.0) dx = h;
gsl_vector_set (x1, i, xi + dx);
fh = GSL_MULTIMIN_FN_EVAL(f, x1);
gsl_vector_set (x1, i, xi);
gsl_vector_set (x1, i, xi - dx);
fl = GSL_MULTIMIN_FN_EVAL(f, x1);
gsl_vector_set (x1, i, xi);
gsl_vector_set (g, i, (fh - fl) / (2.0 * dx));
}
gsl_vector_free (x1);
// return GSL_SUCCESS;
}
double
my_f (const gsl_vector *v, void *params)
{
int *p = (int *)params;
return -(llikelihood_ML(v,parameters.vy,parameters.vyinit, p[0], p[1]));
}
/* The gradient of f, df = (df/dx, df/dy). */
void
my_df (const gsl_vector *v, void *params,
gsl_vector *df)
{
gsl_multimin_function F ;
F.f = my_f;
F.params = params;
F.n = v->size;
gsl_central_diff (&F,v,df);
}
/* Compute both f and df together. */
void
my_fdf (const gsl_vector *x, void *params,
double *f, gsl_vector *df)
{
*f = my_f(x, params);
my_df(x, params, df);
}
int main (void)
{
gsl_rng_env_setup();
// Taus = gsl_rng_default;
r = gsl_rng_alloc (gsl_rng_taus);
gsl_rng_set (r,(unsigned long int)894526599);
int p =2;
int K=9;
int i, j;
// lecture des données
int ctemp=735,cy=1;
gsl_matrix * vytemp = gsl_matrix_calloc(ctemp,cy);
FILE *f = fopen("yext.txt", "rb");
if(!f) { // file couldn't be opened
printf( "Error: file could not be opened" );
}
gsl_matrix_fscanf (f, vytemp);
fclose (f);
int ct=(ctemp-p+1);
parameters.vyinit = gsl_matrix_calloc(ctemp-ct,cy);
parameters.vy = gsl_matrix_calloc(ct,cy);
for(j=0;j<cy;j++)
{
for(i=0;i<p-1;i++)
gsl_matrix_set(parameters.vyinit,i,j,gsl_matrix_ge t(vytemp,i,j));
for(i=0;i<ct;i++)
gsl_matrix_set(parameters.vy,i,j,gsl_matrix_get(vy temp,(i+p-1),j));
}
int R=70;
for(int rep=0;rep<R;rep++)
{
size_t iter = 0;
int status;
const gsl_multimin_fdfminimizer_type *T;
gsl_multimin_fdfminimizer *s;
/* Position of the minimum (1,2), scale factors
10,20, height 30. */
int par[2] = { p, K };
gsl_multimin_function_fdf my_func;
my_func.n = 35;
my_func.f = &my_f;
my_func.df = &my_df;
my_func.fdf = &my_fdf;
my_func.params = ∥
/* Starting point, x = (5,7) */
gsl_vector *m_vPar=gsl_vector_calloc(35);
gsl_vector_set(m_vPar,0,0.021);
gsl_vector_set(m_vPar,1,1.001);
gsl_vector_set(m_vPar,2,0.207);
gsl_vector_set(m_vPar,3,0.908);
gsl_vector_set(m_vPar,4,0.020);
gsl_vector_set(m_vPar,5,1.006);
gsl_vector_set(m_vPar,6,0.201);
gsl_vector_set(m_vPar,7,0.972);
gsl_vector_set(m_vPar,8, 0.215);
gsl_vector_set(m_vPar,9,0.973);
gsl_vector_set(m_vPar,10,0.206);
gsl_vector_set(m_vPar,11,0.973);
gsl_vector_set(m_vPar,12,0.007);
gsl_vector_set(m_vPar,13,0.990);
gsl_vector_set(m_vPar,14,0.167);
gsl_vector_set(m_vPar,15,0.977);
gsl_vector_set(m_vPar,16,0.098);
gsl_vector_set(m_vPar,17, 0.853);
gsl_vector_set(m_vPar,18,44.066);
gsl_vector_set(m_vPar,19,4.193);
gsl_vector_set(m_vPar,20,65.428);
gsl_vector_set(m_vPar,21, 3.918);
gsl_vector_set(m_vPar,22,0.420);
gsl_vector_set(m_vPar,23, 6.440);
gsl_vector_set(m_vPar,24,21.681);
gsl_vector_set(m_vPar,25,58.375);
gsl_vector_set(m_vPar,26,4.695);
gsl_vector_set(m_vPar,27,0.988);
gsl_vector_set(m_vPar,28,0.959);
gsl_vector_set(m_vPar,29,0.980);
gsl_vector_set(m_vPar,30,0.990);
gsl_vector_set(m_vPar,31,0.961);
gsl_vector_set(m_vPar,32, 0.982);
gsl_vector_set(m_vPar,33,0.991);
gsl_vector_set(m_vPar,34,0.968);
T = gsl_multimin_fdfminimizer_vector_bfgs2;
s = gsl_multimin_fdfminimizer_alloc (T, 35);
gsl_multimin_fdfminimizer_set (s, &my_func, m_vPar, 0.005, 0.1);
do
{
iter++;
status = gsl_multimin_fdfminimizer_iterate (s);
if (status)
break;
status = gsl_multimin_test_gradient (s->gradient, 10);
if (status == GSL_SUCCESS)
printf ("Minimum found at:\n");
printf ("%5d %10.5f\n", iter, s->f);
}
while (status == GSL_CONTINUE && iter <35);
printf("Replication %d\n ", rep+1);
if( iter==35)
rep--;
gsl_multimin_fdfminimizer_free (s);
gsl_vector_free (m_vPar);
}
return 0;
}