humpfunc.c 1/1

/*****************************************************************
 * Fits no-interaction model of biomass - species richness data using
 * non-linear least squares with Levenberg-Marquardt method using
 * Numerical Recipes (C) implementation.
 *
 * For the model, see J.Oksanen, J.Ecol. 85, 293-295 (1996).
 * (Describes the corresponding maximum likelihood model with Poisson
 * error).
 *
 * Model is defined so that Species rihness (spno) is a function of
 * biomass with parameters hump, b, alpha:
 *
 *     x = mass/hump       if (mass < hump)
 *     x = (hump/mass)**2  otherwise
 *
 *     Expected(spno) = alpha*log(1+b*x/alpha)
 *
 * humpfunc() below is in the form required by Numerical Recipes
 * function mrqmin() (even with their perverse 1-offset arrays).  The
 * function computes the fitted values (*fv) and derivatives der[] of
 * the parameters par[].
 *
 * The model uses log(alpha) instead of alpha, since this makes
 * iterations much more stable and prevents alpha getting negative
 * values.
 *****************************************************************/
 
#include <math.h>

void humpfunc(float mass, float par[], float *fv, float der[], int npar)
{
  float a1, a2, a3, a4, a5;
  float hump, alpha, b;
  hump = par[1];
  alpha = exp(par[2]); /* Restore alpha: log-alpha estimated */
  b = par[3];
  if (mass <= hump) {
    a1 = mass/hump;
    a5 = -mass/hump/hump;
  }                    
  else { /* when hump < mass */
    a1 = hump/mass*hump/mass;
    a5 = 2.0*hump/mass/mass;
  }
  a2 = b*a1;
  a3 = 1.0 + a2/alpha;
  a4 = log(a3);
  *fv = alpha*a4;
  der[1] = b*a5/a3;
  der[2] = *fv - a2/a3;
  der[3] = a1/a3;
}

Generated by GNU enscript 1.5.0.