C++ RETURN_ARRAYS_INCREMENT函数代码示例

/** negative log likelihood of negative binomial with mean and tau 
\brief Negative binomial with mean=mu and variance = mu*tau
\author Mollie Brooks
\param x observed counts
\param mu is the predicted mean
\param tau is the overdispersion parameter like in the quasi-poisson. should be >1
\return negative log likelihood \f$ -( \ln(\Gamma(x+k))-\ln(\Gamma(k))-\ln(x!)+k\ln(k)+x\ln(\mu)-(k+x)\ln(k+\mu) )\f$
where \f$ k=\mu/(10^{-120}+\tau-1.0) \f$
\ingroup STATLIB
**/
dvariable dnbinom_tau(const dvector& x, const dvar_vector& mu, const dvar_vector& tau)
{
	//the observed counts are in x
	//mu is the predicted mean
	//tau is the overdispersion parameter
	RETURN_ARRAYS_INCREMENT();
	int i,imin,imax;
	imin=x.indexmin();
	imax=x.indexmax();
	dvariable loglike;
	loglike=0.;

	for(i = imin; i<=imax; i++)
	{
		if (value(tau(i))<1.0)
		{
			cerr<<"tau("<<i<<") is <=1.0 in dnbinom_tau()";
			return(0.0);
		}
		
		loglike += log_negbinomial_density(x(i), mu(i), tau(i));
	}
	RETURN_ARRAYS_DECREMENT();
	return(-loglike);
}

/**
 * Description not yet available.
 * \param
 */
dvariable ghk(const dvar_vector& lower,const dvar_vector& upper,
  const dvar_matrix& Sigma, const dmatrix& eps,int _i)
{
  RETURN_ARRAYS_INCREMENT();
  int n=lower.indexmax();
  dvar_matrix ch=choleski_decomp(Sigma);
  dvar_vector l(1,n);
  dvar_vector u(1,n);

  ghk_test(eps,_i);

  dvariable weight=1.0;
  int k=_i;
  {
    l=lower;
    u=upper;
    for (int j=1;j<=n;j++)
    {
      l(j)/=ch(j,j);
      u(j)/=ch(j,j);
      dvariable Phiu=cumd_norm(u(j));
      dvariable Phil=cumd_norm(l(j));
      weight*=Phiu-Phil;
      dvariable eta=inv_cumd_norm((Phiu-Phil)*eps(k,j)+Phil+1.e-30);
      for (int i=j+1;i<=n;i++)
      {
        dvariable tmp=ch(i,j)*eta;
        l(i)-=tmp;
        u(i)-=tmp;
      }
    }
  }
  RETURN_ARRAYS_DECREMENT();
  return weight;
}

/**
 * Description not yet available.
 * \param
 */
dvariable ghk_m(const dvar_vector& upper,const dvar_matrix& Sigma,
  const dmatrix& eps)
{
  RETURN_ARRAYS_INCREMENT();
  int n=upper.indexmax();
  int m=eps.indexmax();
  dvariable ssum=0.0;
  dvar_vector u(1,n);
  dvar_matrix ch=choleski_decomp(Sigma);

  for (int k=1;k<=m;k++)
  {
    dvariable weight=1.0;
    u=upper;
    for (int j=1;j<=n;j++)
    {
      u(j)/=ch(j,j);
      dvariable Phiu=cumd_norm(u(j));
      weight*=Phiu;
      dvariable eta=inv_cumd_norm(1.e-30+Phiu*eps(k,j));
      for (int i=j+1;i<=n;i++)
      {
        dvariable tmp=ch(i,j)*eta;
        u(i)-=tmp;
      }
    }
    ssum+=weight;
  }
  RETURN_ARRAYS_DECREMENT();
  return ssum/m;
}

/**	
	\author Steven James Dean Martell UBC Fisheries Centre
	\date 2011-06-24
	\param  k observed number
	\param  lambda epected mean of the distribution
	\return returns the negative loglikelihood \f$ -k  \ln( \lambda ) - \lambda  + \ln(k!) \f$
	\sa
**/
dvariable dpois(const prevariable& k, const prevariable& lambda)
{
	RETURN_ARRAYS_INCREMENT();
	dvariable tmp = -k*log(lambda)+lambda + gammln(k+1.);
	RETURN_ARRAYS_DECREMENT();
	return tmp;
}

/**
 * Description not yet available.
 * \param
 */
dvariable old_cumd_norm(const prevariable& x)
{
  RETURN_ARRAYS_INCREMENT();
  const double b1=0.319381530;
  const double b2=-0.356563782;
  const double b3=1.781477937;
  const double b4=-1.821255978;
  const double b5=1.330274429;
  const double p=.2316419;
  if (x>=0)
  {
    dvariable u=1./(1+p*x);
    dvariable y=  ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
    dvariable z=1.0-0.3989422804*exp(-.5*x*x)*y;
    RETURN_ARRAYS_DECREMENT();
    return z;
  }
  else
  {
    dvariable w=-x;
    dvariable u=1./(1+p*w);
    dvariable y=  ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
    dvariable z=0.3989422804*exp(-.5*x*x)*y;
    RETURN_ARRAYS_DECREMENT();
    return z;
  }
}

/**
 * Description not yet available.
 * \param
 */
dvariable ghk_choleski(const dvar_vector& lower,const dvar_vector& upper,
  const dvar_matrix& ch, const dmatrix& eps)
{
  RETURN_ARRAYS_INCREMENT();
  int n=lower.indexmax();
  int m=eps.indexmax();
  dvariable ssum=0.0;
  dvar_vector l(1,n);
  dvar_vector u(1,n);

  for (int k=1;k<=m;k++)
  {
    dvariable weight=1.0;
    l=lower;
    u=upper;
    for (int j=1;j<=n;j++)
    {
      l(j)/=ch(j,j);
      u(j)/=ch(j,j);
      dvariable Phiu=cumd_norm(u(j));
      dvariable Phil=cumd_norm(l(j));
      weight*=Phiu-Phil;
      dvariable eta=inv_cumd_norm((Phiu-Phil)*eps(k,j)+Phil);
      for (int i=j+1;i<=n;i++)
      {
        dvariable tmp=ch(i,j)*eta;
        l(i)-=tmp;
        u(i)-=tmp;
      }
    }
    ssum+=weight;
  }
  RETURN_ARRAYS_DECREMENT();
  return ssum/m;
}

/** Negative bionomial density; variable objects.
A local parameter r is used to make it robust.
\f$ r=\frac{\mu}{10.0^{-120}+\tau-1.0} \f$
\ingroup PDF
\param x
\param mu
\param tau
\return Log of NegativeBinomial density. \f$ \frac{\Gamma(x+r)}{\Gamma(r)x!}(\frac{r}{r+\mu})^r(\frac{\mu}{r+\mu})^x \f$
*/
  dvariable negbinomial_density(double x,const prevariable& mu,
    const prevariable& tau)
  {
    if (value(tau)-1.0<=0.0)
    {
      cerr << "tau <=1 in log_negbinomial_density " << endl;
      ad_exit(1);
    }
    RETURN_ARRAYS_INCREMENT();
    dvariable r=mu/(1.e-120+(tau-1));
    dvariable tmp;
    //tmp=exp(gammln(x+r)-gammln(r) -gammln(x+1)
     // +r*log(r)+x*log(mu)-(r+x)*log(r+mu));
    tmp=gammln(x+r);
    tmp-=gammln(r);
    tmp-=gammln(x+1);
    tmp+=r*log(r);
    tmp+=x*log(mu);
    tmp-=(r+x)*log(r+mu);
    tmp=exp(tmp);


    RETURN_ARRAYS_DECREMENT();
    return tmp;
  }

/**
 * Description not yet available.
 * \param
 */
dvariable mean(const dvar_vector& v)
  {
    dvariable tmp;
    RETURN_ARRAYS_INCREMENT();
    tmp=sum(v)/double(v.size());
    RETURN_ARRAYS_DECREMENT();
    return(tmp);
  }

/**
 * Description not yet available.
 * \param
 */
dvariable norm2(const dvar_vector& t1)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable tmp;
    tmp=t1*t1;
    RETURN_ARRAYS_DECREMENT();
    return(tmp);
  }

/**  generalized Ricker function, first parameerization; vectorized
\param x independent variable; data vector
\param x0 ; differentiable vector
\param A ; differentiable vector
\param alpha ; differentiable scalar
\return  \f$ A(\frac{x}{x0}e^{(1.0-\frac{x}{x0})})^{\alpha} \f$
\ingroup ECOL
**/
dvar_vector generalized_Ricker1(const dvector& x,  const dvar_vector& x0,  const dvar_vector& A,  const prevariable& alpha)
{
    RETURN_ARRAYS_INCREMENT();	
	dvar_vector y;
	y=elem_prod(A, pow(elem_prod(elem_div(x, x0), exp(1.0-elem_div(x, x0))), alpha));

	RETURN_ARRAYS_DECREMENT();	
    return (y);
}

/**   monomoleular function; vectorized
\param x independent variable; data vector
\param a ; differentiable vector
\param b ; differentiable vector
\return  \f$ a(1-e^{-bx}) \f$ 
\ingroup ECOL
**/
dvar_vector monomolecular(const dvector& x,  const dvar_vector& a,  const dvar_vector& b)
{
    RETURN_ARRAYS_INCREMENT();	
	dvar_vector y;
	y=elem_prod(a, 1.0-exp(-1.0*elem_prod(b,x))); 

	RETURN_ARRAYS_DECREMENT();	
    return (y);
}

/**  ecologically parameterized logistic function with carrying capacity K;  vectorized
\param t independent variable; data vector
\param K carrying capacity; differentiable vector
\param r growth rate; differentiable vector
\param n0 initial population size at t=0; differentiable scalar
\return  \f$ \frac{K}{1+(\frac{K}{n0}-1)e^{-rt}} \f$ 
\ingroup ECOL
**/
dvar_vector logisticK( const dvector& t,  const dvar_vector& K,  const dvar_vector& r,  const prevariable& n0)
{
    RETURN_ARRAYS_INCREMENT();	
	dvar_vector y;
	y=elem_div(K, 1.0 + elem_prod(K/n0-1.0, exp(-1.0*elem_prod(r,t)))); 

	RETURN_ARRAYS_DECREMENT();	
    return (y);
}

/**  logistic function; vectorized
\param x independent variable; data vector
\param a ; differentiable vector
\param b ; differentiable scalar
\return  \f$ \frac{e^{a+bx}}{(1+e^{a+bx})} \f$
\ingroup ECOL
**/
dvar_vector logistic(const dvector& x,  const dvar_vector& a,  const prevariable& b)
{
    RETURN_ARRAYS_INCREMENT();	
	dvar_vector y;
	y=elem_div(exp(a+b*x), 1.0+exp(a+b*x));

	RETURN_ARRAYS_DECREMENT();	
    return (y);
}

/**  Ricker function; vectorized
\param x independent variable; data vector
\param a ; differentiable vector
\param b ; differentiable vector
\return  \f$ axe^{-bx} \f$
\ingroup ECOL
**/
dvar_vector Ricker(const dvector& x,  const dvar_vector& a,  const dvar_vector& b)
{
    RETURN_ARRAYS_INCREMENT();	
	dvar_vector y;
	y=elem_prod(a, elem_prod(x, exp(-1.0*elem_prod(b, x))));

	RETURN_ARRAYS_DECREMENT();	
    return (y);
}

dvariable dgamma(const prevariable& x, const double& a, const double& b)
  {
    //returns the gamma density with a & b as parameters
    RETURN_ARRAYS_INCREMENT();
    dvariable t1 = 1./(pow(b,a)*mfexp(gammln(a)));
    dvariable t2 = (a-1.)*log(x)-x/b;
    RETURN_ARRAYS_DECREMENT();
    return(t1*mfexp(t2));
  }