help-gsl
[Top][All Lists]
Advanced

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

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


From: pleydell
Subject: Re: Re: [Help-gsl] Passing data to multidim minimization routines,
Date: Sat, 19 Dec 2009 12:00:29 +0100
User-agent: Internet Messaging Program (IMP) H3 (4.0.2)

Here is a rewrite of the code in my last post with two points of clarification
1) The objects in my_f were badly named - the constants were confusingly called
parameters. Now constants are called constants and parameters are called
parameters.
2) The clumsy use of pointers has been refined.

Hopefully the code can now help others starting out with gsl without too much
confusion. I'd be happy for it to be added to the documentation.

#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     */
  double *y;        /* Pointer to response data  */
  double *x;        /* Pointer to covariate data */
};

double my_f (const gsl_vector *parameters, void *Constants) {
  /* Returns -2 * log-likelihood */
  struct data *constants = (struct data *) Constants;
  int i, n;
  double loglik = 0;
  double *x = constants->x;
  double *y = constants->y;
  double c = gsl_vector_get (parameters, 0);
  double m = gsl_vector_get (parameters, 1);
  double d = gsl_vector_get (parameters, 2);
  for (i=0;i<constants->n;i++) {
    loglik += log(gsl_ran_gaussian_pdf(y[i] - (c + m * x[i]), d));
  }
  /* printf("%f",-2*loglik); /\* For Testing Only *\/       */
  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 constants;
  constants.n = n;
  constants.x = x;
  constants.y = y;

  /* Starting point */
  gsl_vector *paras;
  paras = gsl_vector_alloc (nparas);
  for (i=0;i<nparas;i++)
    gsl_vector_set (paras, i, initial_p[i]);

  /* my_f (paras,  &constants);                       /\* Test my_f *\/         
         */

  /* Initialise algorithm parameters */
  size_t iter = 0;
  int status = 0;
  double size;

  /* Set initial STEP SIZES to 1 */
  gsl_vector *ss;
  ss = gsl_vector_alloc (nparas); /* same size as starting vecotr above */
  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 = &constants;

  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(paras);
  gsl_vector_free(ss);

  return status;
}






reply via email to

[Prev in Thread] Current Thread [Next in Thread]