Re: [Help-gsl] Passing data to multidim minimization routines

[David Pleydell]

I only recently started using gsl and this is my first post at this forum.

I was having the same difficulty as Michael Braun who started this thread 
(http://lists.gnu.org/archive/html/help-gsl/2006-12/msg00034.html). While 
Paulo's response was helpful it has taken me most of a day to get a working 
concrete example. In the interest of helping others in the same boat here is a 
concrete example using the Nelder Mead simplex algorithm to perform maximum 
likelihood estimation for a simple linear regression problem. 

Hint, for those from an applied statistics background who might be confused 
"void *params" refers to constants and not the parameters being optimised which 
go in "const gsl_vector *v". OK, here's the program.


/* An example of passing data to gsl_multimin_fminimizer_nmsimplex
   for MLE of a linear regression. */
#include <stdio.h> 
#include <stdlib.h>
#include <gsl/gsl_multimin.h> 
#include <gsl/gsl_errno.h> 
#include <gsl/gsl_randist.h>

struct data {
  int         n;    /* Number of data points */
  gsl_vector *y;    /* Response data */
  gsl_vector *x;    /* Covariate */
};
  
double my_f (const gsl_vector *v, void *params) { 
  /* Returns -2 * log-likelihood */ 
  struct data *p = (struct data *) params;
  int i, n;
  double loglik = 0;
  double xi, yi;
  double c = gsl_vector_get (v, 0);
  double m = gsl_vector_get (v, 1);
  double d = gsl_vector_get (v, 2);
  for (i=0;i<p->n;i++) {
    yi = gsl_vector_get(p->y, i);
    xi = gsl_vector_get(p->x, i);
    loglik += log(gsl_ran_gaussian_pdf(yi - (c + m * xi), d));  
  }
  return -2 * loglik;
}


int main (void) {
  /* Initialise data and model parameters */
  int i, n=10, nparas=3;
  double x[10] = {1,2,3,4,5,6,7,8,9,10};
  double y[10] = {-0.73, 2.00, 2.40, 3.70, 3.80, 
                   7.80, 6.70, 6.40, 9.20, 10.00};
  double initial_p[3] = {0.5,0.5,0.5};
  struct data myStruct;
  myStruct.n = n;
  myStruct.x = gsl_vector_alloc(myStruct.n); 
  myStruct.y = gsl_vector_alloc(myStruct.n); 
  for (i=0;i<myStruct.n;i++) { 
    gsl_vector_set (myStruct.x, i, x[i]);
    gsl_vector_set (myStruct.y, i, y[i]);
  }

  /* 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, 1.0);
  
  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_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, 1e-2);
    
    if (status == GSL_SUCCESS)
    printf ("converged to minimum at\n");
    
    printf ("%5d %10.3e %10.3e %10.3e   f() =%7.3f   size = %.3f\n",
        iter,
        gsl_vector_get (s->x, 0),
        gsl_vector_get (s->x, 1),
        gsl_vector_get (s->x, 2),
        s->fval, size);
  } 
  while (status == GSL_CONTINUE && iter < 100);
  
  gsl_multimin_fminimizer_free (s);  
  gsl_vector_free(myStruct.x);
  gsl_vector_free(myStruct.y);
  gsl_vector_free(paras);
  gsl_vector_free(ss);
  
  return status;
}

/* For debugging compile this with
  gcc mle_with_data.c -ggdb -lgsl -lgslcblas -lm -o mle_with_data
  or otherwise use
  gcc mle_with_data.c -lgsl -lgslcblas -lm -o mle_with_data
*/

/* and run this using
  ./mle_with_data
*/