help-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

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

 From: David Pleydell Subject: Re: [Help-gsl] Passing data to multidim minimization routines Date: Wed, 16 Dec 2009 08:15:43 -0800 (PST)

```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 = {1,2,3,4,5,6,7,8,9,10};
double y = {-0.73, 2.00, 2.40, 3.70, 3.80,
7.80, 6.70, 6.40, 9.20, 10.00};
double initial_p = {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
*/

```

reply via email to

 [Prev in Thread] Current Thread [Next in Thread]
• Re: [Help-gsl] Passing data to multidim minimization routines, David Pleydell <=