hmmer-3.3/easel/esl_gev.c
/* gev_gradient()
 * Computes the gradient of the negative log likelihood of a complete
 * or censored GEV sample; in the API of the CG optimizer.
 */
gev_gradient()
Computes the gradient of the negative log likelihood of a complete
or censored GEV sample; in the API of the CG optimizer.
gev_gradient () Computes the gradient of the negative log likelihood of a complete or censored GEV sample ; in the API of the CG optimizer .
static void
gev_gradient(double *p, int nparam, void *dptr, double *dp)
{
  double mu, w, lambda, alpha;
  struct gev_data *data;
  double *x;
  int    i; 
  double dmu, dw, dalpha;
  double y, ay, ay1, lay1;
  mu     = p[0];
  w      = p[1];   
  lambda = exp(w);
  alpha  = p[2];
  data   = (struct gev_data *) dptr;
  x      = data->x;
  dmu    = 0.;
  dw     = data->n; 
  dalpha = 0.;
  for (i = 0; i < data->n; i++)
    {
      y    = lambda * (x[i]-mu);
      ay   = alpha*y;
      ay1  = 1+ay;		
      lay1 = log(ay1);
      dmu += (alpha+1) / ay1;
      if (fabs(ay) < 1e-12) dmu -= exp(-y);
      else                  dmu -= exp(-(1+1/alpha) * lay1);
      dw -= y*(1+alpha) / ay1;
      if (fabs(ay) < 1e-12) dw += y*exp(-y);
      else                  dw += y*exp(-(1+1/alpha) * lay1);
      dalpha -= (1 + 1/alpha) * y/ay1;
      if (fabs(ay) < 1e-12) {
	dalpha += y/alpha;
	dalpha += y*exp(-y) / (alpha*ay1);
	dalpha -= y*exp(-y) / alpha;
      } else {
	dalpha += lay1 / (alpha*alpha);
	dalpha += y    * exp(-lay1/alpha) / (alpha*ay1);
	dalpha -= lay1 * exp(-lay1/alpha) / (alpha*alpha);
      }
    }
  dmu *= lambda;
  if (data->is_censored)
    {
      y    = lambda * (data->phi - mu);
      ay   = alpha * y;
      ay1  = 1 + ay;
      lay1 = log(ay1);
      if (fabs(ay) < 1e-12) 
	{	
	  dmu    -= data->z * lambda * exp(-y) / ay1;
	  dw     += data->z * y      * exp(-y) / ay1;
	  dalpha -= data->z * exp(-y) * y/alpha * ay/ay1;
	}
      else 
	{	
	  dmu    -= data->z * lambda * exp(-lay1/alpha) / ay1;
	  dw     += data->z * y      * exp(-lay1/alpha) / ay1;
	  dalpha -= data->z * exp(-lay1/alpha) *
	    (lay1/(alpha*alpha) - y/(alpha*ay1));
	}
    }
  dp[0] = -dmu;
  dp[1] = -dw;
  dp[2] = -dalpha;
  return;
}
static void gev_gradient ( double * p , int nparam , void * dptr , double * dp ) { double mu , w , lambda , alpha ; struct gev_data * data ; double * x ; int i ; double dmu , dw , dalpha ; double y , ay , ay1 , lay1 ; mu = p [ 0 ] ; w = p [ 1 ] ; lambda = exp ( w ) ; alpha = p [ 2 ] ; data = ( struct gev_data * ) dptr ; x = data -> x ; dmu = 0. ; dw = data -> n ; dalpha = 0. ; for ( i = 0 ; i < data -> n ; i ++ ) { y = lambda * ( x [ i ] - mu ) ; ay = alpha * y ; ay1 = 1 + ay ; lay1 = log ( ay1 ) ; dmu += ( alpha + 1 ) / ay1 ; if ( fabs ( ay ) < 1e-12 ) dmu -= exp ( - y ) ; else dmu -= exp ( - ( 1 + 1 / alpha ) * lay1 ) ; dw -= y * ( 1 + alpha ) / ay1 ; if ( fabs ( ay ) < 1e-12 ) dw += y * exp ( - y ) ; else dw += y * exp ( - ( 1 + 1 / alpha ) * lay1 ) ; dalpha -= ( 1 + 1 / alpha ) * y / ay1 ; if ( fabs ( ay ) < 1e-12 ) { dalpha += y / alpha ; dalpha += y * exp ( - y ) / ( alpha * ay1 ) ; dalpha -= y * exp ( - y ) / alpha ; } else { dalpha += lay1 / ( alpha * alpha ) ; dalpha += y * exp ( - lay1 / alpha ) / ( alpha * ay1 ) ; dalpha -= lay1 * exp ( - lay1 / alpha ) / ( alpha * alpha ) ; } } dmu *= lambda ; if ( data -> is_censored ) { y = lambda * ( data -> phi - mu ) ; ay = alpha * y ; ay1 = 1 + ay ; lay1 = log ( ay1 ) ; if ( fabs ( ay ) < 1e-12 ) { dmu -= data -> z * lambda * exp ( - y ) / ay1 ; dw += data -> z * y * exp ( - y ) / ay1 ; dalpha -= data -> z * exp ( - y ) * y / alpha * ay / ay1 ; } else { dmu -= data -> z * lambda * exp ( - lay1 / alpha ) / ay1 ; dw += data -> z * y * exp ( - lay1 / alpha ) / ay1 ; dalpha -= data -> z * exp ( - lay1 / alpha ) * ( lay1 / ( alpha * alpha ) - y / ( alpha * ay1 ) ) ; } } dp [ 0 ] = - dmu ; dp [ 1 ] = - dw ; dp [ 2 ] = - dalpha ; return ; }
gev_gradient()
 Computes the gradient of the negative log likelihood of a complete
 or censored GEV sample; in the API of the CG optimizer.
gev_gradient () Computes the gradient of the negative log likelihood of a complete or censored GEV sample ; in the API of the CG optimizer .
/* Unpack what the optimizer gave us */ /* w is a c.o.v. to allow unconstrained opt of lambda>0 */ /* d/dw, term1 */ /* 1+ay=1, for ay < DBL_EPSILON */ /* d/dmu, term1. (will become 1, for small alpha.) */ /* d/dmu, term2. For tiny ay, use log(1+x) ~ x to simplify. */ /* d/dw, term2. converges to -y, for small alpha. */ /* d/dw, term2. For tiny ay, use log(1+x) ~ x to simplify. */ /* d/dalpha, term1
       */ /* d/dalpha, terms 2,3,4: for tiny ay, simplify.
       * d/dalpha will go to +/-inf for alpha ~ 0, so watch out.
       */ /* Add the terms that come from the censored data gradient,
   * if it's a censored dataset.
   */ /* special case of small alpha, converging towards Gumbel */ /* normal case */ /* Return the negative gradient, because we're minimizing NLL,
   * not maximizing LL.
   */