#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;
}