qfc_est.cpp

Go to the documentation of this file.

  #include "qfclib.h"




  dvar_vector logitProp(const dvar_vector& a)
  {
    int dim; 
    dim=a.size()+1;
    dvar_vector p(1,dim);
    dvar_vector expa=mfexp(a);
    p(1,dim-1)=expa/(1.+sum(expa));
    //p(dim)=1.-sum(p(1,dim-1));  //original version was buggy for the last term
    p(dim)=1./(1.+sum(expa));
    return p;
  }
  df1b2vector logitProp(const df1b2vector& a)
  {
    int dim; 
    dim=size_count(a)+1;
    df1b2vector p(1,dim);
    df1b2vector expa=mfexp(a);
    p(1,dim-1)=expa/(1.+sum(expa));
    //p(dim)=1.-sum(p(1,dim-1));  //original version was buggy for the last term
    p(dim)=1./(1.+sum(expa));
    return p;
  }




  dvar_vector invLogitProp(const dvar_vector& p)
  {
    int dim; 
    dim=p.size()-1;
    dvar_vector a(1,dim); 
    dvar_vector lp(1,dim+1);
  
    lp=log(p+EPS);//take log for each p
    a=lp(1,dim)-lp(dim+1);//each subtract the last one
    return a;
  }
  df1b2vector invLogitProp(const df1b2vector& p)
  {
    int dim; 
    dim=size_count(p)-1;
    df1b2vector a(1,dim); 
    df1b2vector lp(1,dim+1);
  
    lp=log(p+EPS);//take log for each p
    a=lp(1,dim)-lp(dim+1);//each subtract the last one
    return a;
  }





  dvar_vector normalize_p(dvar_vector& p, dvariable fpen)
  {
    dvariable psum=sum(p);
    p/=psum; // Now the p will sum to 1
    fpen+=1000.*square(log(psum)); //penalty
    return p;
  }
  df1b2vector normalize_p(df1b2vector& p, df1b2variable fpen)
  {
    df1b2variable psum=sum(p);
    p/=psum; // Now the p will sum to 1
    fpen+=1000.*square(log(psum)); //penalty
    return p;
  }




  dvar_vector posfun(dvar_vector& x,const dvector& eps, dvariable& pen) //x and eps both are vector
  { 
    for(int i=x.indexmin();i<=x.indexmax();i++)
    {
      x(i)=posfun(x(i),eps(i),pen);
    }
    return x;
  }
  df1b2vector posfun(df1b2vector& x,const dvector& eps, df1b2variable& pen) //x and eps both are vector
  { 
    for(int i=x.indexmin();i<=x.indexmax();i++)
    {
      x(i)=posfun(x(i),eps(i),pen);
    }
    return x;
  }




  dvar_matrix posfun(dvar_matrix& x,const dmatrix & eps, dvariable& pen)
  {
    for(int i=x.rowmin();i<=x.rowmax();i++)
    {
      for(int j=x.colmin();j<=x.colmax();j++)
      {
        x(i,j)=posfun(x(i,j),eps(i,j),pen);
      }
    }
    return x;
  }
  df1b2matrix posfun(df1b2matrix& x,const dmatrix & eps, df1b2variable& pen)
  {
    for(int i=x.rowmin();i<=x.rowmax();i++)
    {
      for(int j=x.colmin();j<=x.colmax();j++)
      {
        x(i,j)=posfun(x(i,j),eps(i,j),pen);
      }
    }
    return x;
  }




  
 
  dvariable mf_upper_bound2(const dvariable & x,const double fmax, dvariable & fpen)
  {
    if (value(x)<=fmax) return x;
    else
    {
      fpen+=.01*square(x-fmax);
      return x-x/(2.-fmax/x)+fmax;
    }
  }

  df1b2variable mf_upper_bound2(const df1b2variable & x,const double fmax, df1b2variable & fpen)
  {
    if (value(x)<=fmax) return x;
    else
    {
      fpen+=.01*square(x-fmax);
      return x-x/(2.-fmax/x)+fmax;
    }
  }





  dvar_vector boundp(const dvar_vector & x, const double fmin, const double fmax, const dvariable & fpen )
  {
    dvar_vector t(x.indexmin(),x.indexmax());
    for(int i=x.indexmin();i<=x.indexmax();i++)
      t(i)=boundp(x(i),fmin,fmax,fpen);
    return (t);
  }
  df1b2vector boundp(const df1b2vector & x, const double fmin, const double fmax, const df1b2variable & fpen )
  {
    df1b2vector t(x.indexmin(),x.indexmax());
    for(int i=x.indexmin();i<=x.indexmax();i++)
      t(i)=boundp(x(i),fmin,fmax,fpen);
    return (t);
  }







  df1b2vector log_comb(const df1b2vector & n,const df1b2vector & x)
  {
    df1b2vector tmp(n.indexmin(),n.indexmax());
    for(int i=n.indexmin();i<=n.indexmax();i++)
      tmp(i)=log_comb(n(i),x(i));
    return tmp;
  }
  df1b2vector log_comb(const df1b2vector & n,const dvector & x)
  {
    df1b2vector tmp(n.indexmin(),n.indexmax());
    for(int i=n.indexmin();i<=n.indexmax();i++)
      tmp(i)=log_comb(n(i),x(i));
    return tmp;
  }
  df1b2vector log_comb(const dvector & n,const df1b2vector & x)
  {
    df1b2vector tmp(n.indexmin(),n.indexmax());
    for(int i=n.indexmin();i<=n.indexmax();i++)
      tmp(i)=log_comb(n(i),x(i));
    return tmp;
  }  







  //
  //                                    Negative Binomial 
  //

  double  nllNegativeBinomial(const double obs, const double m, const double s )
  {
    double nll=0; 
    nll= -gammln(obs+s)+ gammln(s)- s*log(s/(m+s)+EPS)- obs*log(m/(m+s)+EPS)+ gammln(obs+1.);
    return nll;
  }
  dvariable  nllNegativeBinomial(const dvariable & obs, const double m, const double s )
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0; 
    nll= -gammln(obs+s)+ gammln(s)- s*log(s/(m+s)+EPS)- obs*log(m/(m+s)+EPS)+ gammln(obs+1.);
    RETURN_ARRAYS_DECREMENT();
    return nll;
  } 
  dvariable  nllNegativeBinomial(const double obs, const dvariable & m, const dvariable & s)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0; 
    nll= -gammln(obs+s)+ gammln(s)- s*log(s/(m+s)+EPS)- obs*log(m/(m+s)+EPS)+ gammln(obs+1.);
    RETURN_ARRAYS_DECREMENT();
    return nll;
  } 
  dvariable  nllNegativeBinomial(const dvariable & obs, const dvariable & m, const double s)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0; 
    nll= -gammln(obs+s)+ gammln(s)- s*log(s/(m+s)+EPS)- obs*log(m/(m+s)+EPS)+ gammln(obs+1.);
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable  nllNegativeBinomial(const dvariable & obs, const double m, const dvariable & s)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0; 
    nll= -gammln(obs+s)+ gammln(s)- s*log(s/(m+s)+EPS)- obs*log(m/(m+s)+EPS)+ gammln(obs+1.);
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable  nllNegativeBinomial(const dvariable & obs, const dvariable & m, const dvariable & s)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0; 
    nll= -gammln(obs+s)+ gammln(s)- s*log(s/(m+s)+EPS)- obs*log(m/(m+s)+EPS)+ gammln(obs+1.);
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }

  df1b2variable  nllNegativeBinomial(const df1b2variable & obs, const double m, const double s)
  {
    df1b2variable nll=.0; 
    nll= -gammln(obs+s)+ gammln(s)- s*log(s/(m+s)+EPS)- obs*log(m/(m+s)+EPS)+ gammln(obs+1.);
    return nll;
  }
  df1b2variable  nllNegativeBinomial(const double obs, const df1b2variable & m, const df1b2variable & s)
  {
    df1b2variable nll=.0; 
    nll= -gammln(obs+s)+ gammln(s)- s*log(s/(m+s)+EPS)- obs*log(m/(m+s)+EPS)+ gammln(obs+1.);
    return nll;
  }
  df1b2variable  nllNegativeBinomial(const df1b2variable & obs, const df1b2variable & m, const double s)
  {
    df1b2variable nll=.0; 
    nll= -gammln(obs+s)+ gammln(s)- s*log(s/(m+s)+EPS)- obs*log(m/(m+s)+EPS)+ gammln(obs+1.);
    return nll;
  }
  df1b2variable  nllNegativeBinomial(const df1b2variable & obs, const double m, const df1b2variable & s)
  {
    df1b2variable nll=.0; 
    nll= -gammln(obs+s)+ gammln(s)- s*log(s/(m+s)+EPS)- obs*log(m/(m+s)+EPS)+ gammln(obs+1.);
    return nll;
  }
  df1b2variable  nllNegativeBinomial(const df1b2variable & obs, const df1b2variable & m, const df1b2variable & s)
  {
    df1b2variable nll=.0; 
    nll= -gammln(obs+s)+ gammln(s)- s*log(s/(m+s)+EPS)- obs*log(m/(m+s)+EPS)+ gammln(obs+1.);
    return nll;
  }





  dvariable  nllNegativeBinomial(const dvector & obs, const dvariable & m, const dvariable & s)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0; 
    double n=double(obs.size());
    nll= -sum(gammln(obs+s))+ n*gammln(s)- n*s*log(s/(m+s)+EPS)- sum(obs)*log(m/(m+s)+EPS)+ sum(gammln(obs+1.));
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable  nllNegativeBinomial(const dvar_vector & obs, const double m, const double s)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0; 
    double n=double(obs.size());
    nll= -sum(gammln(obs+s))+ n*gammln(s)- n*s*log(s/(m+s)+EPS)- sum(obs)*log(m/(m+s)+EPS)+ sum(gammln(obs+1.));
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable  nllNegativeBinomial(const dvar_vector & obs, const dvariable & m, const dvariable & s)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0; 
    double n=double(obs.size());
    nll= -sum(gammln(obs+s))+ n*gammln(s)- n*s*log(s/(m+s)+EPS)- sum(obs)*log(m/(m+s)+EPS)+ sum(gammln(obs+1.));
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }

  df1b2variable  nllNegativeBinomial(const df1b2vector & obs, const double m, const double s)
  {
    df1b2variable nll=.0; 
    double n=double(size_count(obs));
    nll= -sum(gammln(obs+s))+ n*gammln(s)- n*s*log(s/(m+s)+EPS)- sum(obs)*log(m/(m+s)+EPS)+ sum(gammln(obs+1.));
    return nll;
  }
  df1b2variable  nllNegativeBinomial(const dvector & obs, const df1b2variable & m, const df1b2variable & s)
  {
    df1b2variable nll=.0; 
    double n=double(size_count(obs));
    nll= -sum(gammln(obs+s))+ n*gammln(s)- n*s*log(s/(m+s)+EPS)- sum(obs)*log(m/(m+s)+EPS)+ sum(gammln(obs+1.));
    return nll;
  }
  df1b2variable  nllNegativeBinomial(const df1b2vector & obs, const df1b2variable & m, const df1b2variable & s)
  {
    df1b2variable nll=.0; 
    double n=double(size_count(obs));
    nll= -sum(gammln(obs+s))+ n*gammln(s)- n*s*log(s/(m+s)+EPS)- sum(obs)*log(m/(m+s)+EPS)+ sum(gammln(obs+1.));
    return nll;
  }




  //
  //                                    Negative Binomial with tau
  //

  double nllNegativeBinomial2(const double obs, const double m, const double tau)
  {
    return -log_density_negbinomial (obs,m,tau); //use admb built in function
  } 
  dvariable nllNegativeBinomial2(const double obs, const dvariable & m, const dvariable & tau)
  {
    return -1.*log_negbinomial_density(obs,m,tau); //use admb built in function
  }   
  dvariable nllNegativeBinomial2(const dvariable & obs, const double m, const double tau)
  {
    return nllNegativeBinomial(obs,m,m/(tau-1.+EPS));
  }   
  dvariable nllNegativeBinomial2(const dvariable & obs, const dvariable & m, const double tau)
  {
    return nllNegativeBinomial(obs,m,m/(tau-1.+EPS));
  } 
  dvariable nllNegativeBinomial2(const dvariable & obs, const double m, const dvariable & tau)
  {
    return nllNegativeBinomial(obs,m,m/(tau-1.+EPS));
  } 
  dvariable nllNegativeBinomial2(const dvariable & obs, const dvariable & m, const dvariable & tau)
  {
    return nllNegativeBinomial(obs, m,m/(tau-1.+EPS));
  } 

  df1b2variable nllNegativeBinomial2(const df1b2variable & obs, const double m, const double tau)
  {
    return nllNegativeBinomial(obs,m,m/(tau-1.+EPS));
  }
  df1b2variable nllNegativeBinomial2(const double obs, const df1b2variable & m, const df1b2variable & tau)
  {
    return nllNegativeBinomial(obs,m,m/(tau-1.+EPS));
  }
  df1b2variable nllNegativeBinomial2(const df1b2variable & obs, const df1b2variable & m, const df1b2variable & tau)
  {
    return nllNegativeBinomial(obs,m,m/(tau-1.+EPS));
  }






  //
  //                                    Multinomial
  //

  double  nllMultiNomial(const dvector & obsN, const dvector & p)
  {
    double nll=0;  
    nll= -1.*(obsN*log(p+EPS) + gammln(sum(obsN)+1.) - sum(gammln(obsN+1.))); //full likelihood
    return nll;
  } 
  dvariable  nllMultiNomial(const dvector & obsN, const dvar_vector & p)
  {
    RETURN_ARRAYS_INCREMENT(); 
    dvariable nll=0;  
    nll= -1.*(obsN*log(p+EPS) + gammln(sum(obsN)+1.) - sum(gammln(obsN+1.))); //full likelihood
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable  nllMultiNomial(const dvar_vector & obsN, const dvector & p)
  {
    RETURN_ARRAYS_INCREMENT(); 
    dvariable nll=0;  
    nll= -1.*(obsN*log(p+EPS) + gammln(sum(obsN)+1.) - sum(gammln(obsN+1.))); //full likelihood
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }   
  dvariable  nllMultiNomial(const dvar_vector & obsN, const dvar_vector & p)
  {
    RETURN_ARRAYS_INCREMENT(); 
    dvariable nll=0;  
    nll= -1.*(obsN*log(p+EPS) + gammln(sum(obsN)+1.) - sum(gammln(obsN+1.))); //full likelihood
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }    

  df1b2variable  nllMultiNomial(const df1b2vector & obsN, const dvector & p)
  {
    df1b2variable nll=.0;  
    nll= -1.*(obsN*log(p+EPS) + gammln(sum(obsN)+1.) - sum(gammln(obsN+1.))); //full likelihood
    return nll;
  }
  df1b2variable  nllMultiNomial(const dvector & obsN, const df1b2vector & p)
  {
    df1b2variable nll=.0;  
    nll= -1.*(obsN*log(p+EPS) + gammln(sum(obsN)+1.) - sum(gammln(obsN+1.))); //full likelihood
    return nll;
  }  
  df1b2variable  nllMultiNomial(const df1b2vector & obsN, const df1b2vector & p)
  {
    df1b2variable nll=.0;  
    nll= -1.*(obsN*log(p+EPS) + gammln(sum(obsN)+1.) - sum(gammln(obsN+1.))); //full likelihood
    return nll;
  }




  //
  //                                    Dirichlet
  //

  double nllDirichlet(const dvector & p, const dvector & shape)
  {
    double nll=0;  
    nll= -1.*(gammln(sum(shape)) - sum(gammln(shape)) + (shape-1.)*log(p+EPS));
    return nll;
  }
  dvariable nllDirichlet(const dvar_vector & p, const dvector & shape)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    nll= -1.*(gammln(sum(shape)) - sum(gammln(shape)) + (shape-1.)*log(p+EPS));
    RETURN_ARRAYS_DECREMENT();
    return nll;
  } 
  dvariable nllDirichlet(const dvector & p, const dvar_vector & shape)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    nll= -1.*(gammln(sum(shape)) - sum(gammln(shape)) + (shape-1.)*log(p+EPS));
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllDirichlet(const dvar_vector & p, const dvar_vector & shape)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    nll= -1.*(gammln(sum(shape)) - sum(gammln(shape)) + (shape-1.)*log(p+EPS));
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }

  df1b2variable nllDirichlet(const df1b2vector & p, const dvector & shape)
  {
    df1b2variable nll=.0;  
    nll= -1.*(gammln(sum(shape)) - sum(gammln(shape)) + (shape-1.)*log(p+EPS));
    return nll;
  }  
  df1b2variable nllDirichlet(const dvector & p, const df1b2vector & shape)
  {
    df1b2variable nll=.0;  
    nll= -1.*(gammln(sum(shape)) - sum(gammln(shape)) + (shape-1.)*log(p+EPS));
    return nll;
  }
  df1b2variable nllDirichlet(const df1b2vector & p, const df1b2vector & shape)
  {
    df1b2variable nll=.0;  
    nll= -1.*(gammln(sum(shape)) - sum(gammln(shape)) + (shape-1.)*log(p+EPS));
    return nll;
  }  

 


  //
  //                                    Gamma
  //

  double nllGamma(const double & x, const double a, const double b)
  {
    return -1.*log_gamma_density(x,a,b);  //call admb built in function
  }
  dvariable nllGamma(const dvariable & x, const double a, const double b)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    nll= gammln(a) - a*log(b) -(a-1.)*log(x+EPS) + b*x;  
    RETURN_ARRAYS_DECREMENT();
    return nll;
  } 
  dvariable nllGamma(const double x, const dvariable & a, const dvariable & b)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    nll= gammln(a) - a*log(b) -(a-1.)*log(x+EPS) + b*x;  
    RETURN_ARRAYS_DECREMENT();
    return nll;
  } 
  dvariable nllGamma(const dvariable & x, const dvariable & a, const double b)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    nll= gammln(a) - a*log(b) -(a-1.)*log(x+EPS) + b*x;  
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllGamma(const dvariable & x, const double a, const dvariable & b)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    nll= gammln(a) - a*log(b) -(a-1.)*log(x+EPS) + b*x;  
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllGamma(const dvariable & x, const dvariable & a, const dvariable & b)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    nll= gammln(a) - a*log(b) -(a-1.)*log(x+EPS) + b*x;  
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }

  df1b2variable nllGamma(const df1b2variable & x, const double a, const double b)
  {
    df1b2variable nll=.0;  
    nll= gammln(a) - a*log(b) -(a-1.)*log(x+EPS) + b*x;  
    return nll;
  }
  df1b2variable nllGamma(const double x, const df1b2variable & a, const df1b2variable & b)
  {
    df1b2variable nll=.0;  
    nll= gammln(a) - a*log(b) -(a-1.)*log(x+EPS) + b*x;  
    return nll;
  }
  df1b2variable nllGamma(const df1b2variable & x, const df1b2variable & a, const double b)
  {
    df1b2variable nll=.0;  
    nll= gammln(a) - a*log(b) -(a-1.)*log(x+EPS) + b*x;  
    return nll;
  }
  df1b2variable nllGamma(const df1b2variable & x, const double a, const df1b2variable & b)
  {
    df1b2variable nll=.0;  
    nll= gammln(a) - a*log(b) -(a-1.)*log(x+EPS) + b*x;  
    return nll;
  }
  df1b2variable nllGamma(const df1b2variable & x, const df1b2variable & a, const df1b2variable & b)
  {
    df1b2variable nll=.0;  
    nll= gammln(a) - a*log(b) -(a-1.)*log(x+EPS) + b*x;  
    return nll;
  }



  double nllGamma(const dvector & x, const double a, const double b)
  {
    double nll=0;  
    double n=double(x.size());
    nll=n*gammln(a) - n*a*log(b) -(a-1.)*sum(log(x+EPS)) + b*sum(x);  
    return nll;
  }
  dvariable nllGamma(const dvar_vector & x, const double a, const double b)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    double n=double(x.size());
    nll=n*gammln(a) - n*a*log(b) -(a-1.)*sum(log(x+EPS)) + b*sum(x);  
    RETURN_ARRAYS_DECREMENT();
    return nll;
  } 
  dvariable nllGamma(const dvector & x, const dvariable & a, const dvariable & b)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    double n=double(x.size());
    nll=n*gammln(a) - n*a*log(b) -(a-1.)*sum(log(x+EPS)) + b*sum(x);  
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllGamma(const dvar_vector & x, const dvariable & a, const dvariable & b)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    double n=double(x.size());
    nll=n*gammln(a) - n*a*log(b) -(a-1.)*sum(log(x+EPS)) + b*sum(x);  
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }

  df1b2variable nllGamma(const df1b2vector & x, const double a, const double b)
  {
    df1b2variable nll=.0;  
    double n=double(size_count(x));
    nll=n*gammln(a) - n*a*log(b) -(a-1.)*sum(log(x+EPS)) + b*sum(x);  
    return nll;
  }  
  df1b2variable nllGamma(const dvector & x, const df1b2variable & a, const df1b2variable & b)
  {
    df1b2variable nll=.0;  
    double n=double(size_count(x));
    nll=n*gammln(a) - n*a*log(b) -(a-1.)*sum(log(x+EPS)) + b*sum(x);  
    return nll;
  }    
  df1b2variable nllGamma(const df1b2vector & x, const df1b2variable & a, const df1b2variable & b)
  {
    df1b2variable nll=.0;  
    double n=double(size_count(x));
    nll=n*gammln(a) - n*a*log(b) -(a-1.)*sum(log(x+EPS)) + b*sum(x);  
    return nll;
  }  




  //
  //                                    Beta
  //
   
  double nllBeta(const double x, const double a, const double b)
  {
    double nll=0;  
    nll=gammln(a)+gammln(b)-gammln(a+b)-(a-1.)*log(x+EPS)-(b-1.)*log(1.-x+EPS);  
    return nll;
  }
  dvariable nllBeta(const dvariable & x, const double a, const double b)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    nll=gammln(a)+gammln(b)-gammln(a+b)-(a-1.)*log(x+EPS)-(b-1.)*log(1.-x+EPS);  
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }   
  dvariable nllBeta(const double x, const dvariable & a, const dvariable & b)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    nll=gammln(a)+gammln(b)-gammln(a+b)-(a-1.)*log(x+EPS)-(b-1.)*log(1.-x+EPS);  
    RETURN_ARRAYS_DECREMENT();
    return nll;
  } 
  dvariable nllBeta(const dvariable & x, const dvariable & a, const double b)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    nll=gammln(a)+gammln(b)-gammln(a+b)-(a-1.)*log(x+EPS)-(b-1.)*log(1.-x+EPS);  
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }  
  dvariable nllBeta(const dvariable & x, const double a, const dvariable & b)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    nll=gammln(a)+gammln(b)-gammln(a+b)-(a-1.)*log(x+EPS)-(b-1.)*log(1.-x+EPS);  
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }  
  dvariable nllBeta(const dvariable & x, const dvariable & a, const dvariable & b)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    nll=gammln(a)+gammln(b)-gammln(a+b)-(a-1.)*log(x+EPS)-(b-1.)*log(1.-x+EPS);  
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }  

  df1b2variable nllBeta(const df1b2variable & x, const double a, const double b)
  {
    df1b2variable nll=.0;  
    nll=gammln(a)+gammln(b)-gammln(a+b)-(a-1.)*log(x+EPS)-(b-1.)*log(1.-x+EPS);  
    return nll;
  }   
  df1b2variable nllBeta(const df1b2variable & x, const df1b2variable & a, const double b)
  {
    df1b2variable nll=.0;  
    nll=gammln(a)+gammln(b)-gammln(a+b)-(a-1.)*log(x+EPS)-(b-1.)*log(1.-x+EPS);  
    return nll;
  } 
  df1b2variable nllBeta(const df1b2variable & x, const double a, const df1b2variable & b)
  {
    df1b2variable nll=.0;  
    nll=gammln(a)+gammln(b)-gammln(a+b)-(a-1.)*log(x+EPS)-(b-1.)*log(1.-x+EPS);  
    return nll;
  } 
  df1b2variable nllBeta(const double x, const df1b2variable & a, const df1b2variable & b)
  {
    df1b2variable nll=.0;  
    nll=gammln(a)+gammln(b)-gammln(a+b)-(a-1.)*log(x+EPS)-(b-1.)*log(1.-x+EPS);  
    return nll;
  } 
  df1b2variable nllBeta(const df1b2variable & x, const df1b2variable & a, const df1b2variable & b)
  {
    df1b2variable nll=.0;  
    nll=gammln(a)+gammln(b)-gammln(a+b)-(a-1.)*log(x+EPS)-(b-1.)*log(1.-x+EPS);  
    return nll;
  } 



  dvariable nllBeta(const dvar_vector & x, const double a, const double b)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    double n=double(x.size());
    nll=n*gammln(a)+n*gammln(b)-n*gammln(a+b)-(a-1.)*sum(log(x+EPS))-(b-1.)*sum(log(1.-x+EPS));  
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllBeta(const dvector & x, const dvariable & a, const dvariable & b)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    double n=double(x.size());
    nll=n*gammln(a)+n*gammln(b)-n*gammln(a+b)-(a-1.)*sum(log(x+EPS))-(b-1.)*sum(log(1.-x+EPS));  
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllBeta(const dvar_vector & x, const dvariable & a, const dvariable & b)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    double n=double(x.size());
    nll=n*gammln(a)+n*gammln(b)-n*gammln(a+b)-(a-1.)*sum(log(x+EPS))-(b-1.)*sum(log(1.-x+EPS));  
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  
  df1b2variable nllBeta(const df1b2vector & x, const double a, const double b)
  {
    df1b2variable nll=.0;  
    double n=double(size_count(x));
    nll=n*gammln(a)+n*gammln(b)-n*gammln(a+b)-(a-1.)*sum(log(x+EPS))-(b-1.)*sum(log(1.-x+EPS));  
    return nll;
  }  
  df1b2variable nllBeta(const dvector & x, const df1b2variable & a, const df1b2variable & b)
  {
    df1b2variable nll=.0;  
    double n=double(size_count(x));
    nll=n*gammln(a)+n*gammln(b)-n*gammln(a+b)-(a-1.)*sum(log(x+EPS))-(b-1.)*sum(log(1.-x+EPS));  
    return nll;
  }   
  df1b2variable nllBeta(const df1b2vector & x, const df1b2variable & a, const df1b2variable & b)
  {
    df1b2variable nll=.0;  
    double n=double(size_count(x));
    nll=n*gammln(a)+n*gammln(b)-n*gammln(a+b)-(a-1.)*sum(log(x+EPS))-(b-1.)*sum(log(1.-x+EPS));  
    return nll;
  } 




  //
  //                                    Normal
  //

  double nllNormal(const double x, const double mu, const double sigma)
  {
    double nll=0;  
    nll=log(sigma*sqrt(2.*M_PI)+EPS) +0.5*square((x-mu)/(sigma+EPS));  
    return nll;
  }
  dvariable nllNormal(const dvariable & x, const double mu, const double sigma)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    nll=log(sigma*sqrt(2.*M_PI)+EPS) +0.5*square((x-mu)/(sigma+EPS));  
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllNormal(const double x, const dvariable & mu, const dvariable & sigma)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    nll=log(sigma*sqrt(2.*M_PI)+EPS) +0.5*square((x-mu)/(sigma+EPS));  
    RETURN_ARRAYS_DECREMENT();
    return nll;
  } 
  dvariable nllNormal(const dvariable & x, const dvariable & mu, const double sigma)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    nll=log(sigma*sqrt(2.*M_PI)+EPS) +0.5*square((x-mu)/(sigma+EPS));  
    RETURN_ARRAYS_DECREMENT();
    return nll;
  } 
  dvariable nllNormal(const dvariable & x, const double mu, const dvariable & sigma)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    nll=log(sigma*sqrt(2.*M_PI)+EPS) +0.5*square((x-mu)/(sigma+EPS));  
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllNormal(const dvariable & x, const dvariable & mu, const dvariable & sigma)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    nll=log(sigma*sqrt(2.*M_PI)+EPS) +0.5*square((x-mu)/(sigma+EPS));  
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }

  df1b2variable nllNormal(const df1b2variable & x, const double mu, const double sigma)
  {
    df1b2variable nll=.0;  
    nll=log(sigma*sqrt(2.*M_PI)+EPS) +0.5*square((x-mu)/(sigma+EPS));  
    return nll;
  }
  df1b2variable nllNormal(const df1b2variable & x, const df1b2variable & mu, const double sigma)
  {
    df1b2variable nll=.0;  
    nll=log(sigma*sqrt(2.*M_PI)+EPS) +0.5*square((x-mu)/(sigma+EPS));  
    return nll;
  }
  df1b2variable nllNormal(const df1b2variable & x, const double mu, const df1b2variable & sigma)
  {
    df1b2variable nll=.0;  
    nll=log(sigma*sqrt(2.*M_PI)+EPS) +0.5*square((x-mu)/(sigma+EPS));  
    return nll;
  }
  df1b2variable nllNormal(const double x, const df1b2variable & mu, const df1b2variable & sigma)
  {
    df1b2variable nll=.0;  
    nll=log(sigma*sqrt(2.*M_PI)+EPS) +0.5*square((x-mu)/(sigma+EPS));  
    return nll;
  }
  df1b2variable nllNormal(const df1b2variable & x, const df1b2variable & mu, const df1b2variable & sigma)
  {
    df1b2variable nll=.0;  
    nll=log(sigma*sqrt(2.*M_PI)+EPS) +0.5*square((x-mu)/(sigma+EPS));  
    return nll;
  }



  double nllNormal(const dvector & x, const double mu, const double sigma)
  {
    double nll=0;  
    double n=double(x.size());
    nll=n*log(sigma*sqrt(2.*M_PI)+EPS) +0.5*norm2((x-mu)/(sigma+EPS));  
    return nll;
  }
  dvariable nllNormal(const dvar_vector& x, const double mu, const double sigma)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    double n=double(x.size());
    nll=n*log(sigma*sqrt(2.*M_PI)+EPS) +0.5*norm2((x-mu)/(sigma+EPS));  
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllNormal(const dvector & x, const dvariable & mu, const dvariable & sigma)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    double n=double(x.size());
    nll=n*log(sigma*sqrt(2.*M_PI)+EPS) +0.5*norm2((x-mu)/(sigma+EPS));  
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllNormal(const dvar_vector & x, const dvariable & mu, const dvariable & sigma)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    double n=double(x.size());
    nll=n*log(sigma*sqrt(2.*M_PI)+EPS) +0.5*norm2((x-mu)/(sigma+EPS));  
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
 
  df1b2variable nllNormal(const df1b2vector & x, const double mu, const double sigma)
  {
    df1b2variable nll=.0;  
    double n=double(size_count(x));
    nll=n*log(sigma*sqrt(2.*M_PI)+EPS) +0.5*norm2((x-mu)/(sigma+EPS));  
    return nll;
  }
  df1b2variable nllNormal(const dvector & x, const df1b2variable & mu, const df1b2variable & sigma)
  {
    df1b2variable nll=.0;  
    double n=double(size_count(x));
    nll=n*log(sigma*sqrt(2.*M_PI)+EPS) +0.5*norm2((x-mu)/(sigma+EPS));  
    return nll;
  }
  df1b2variable nllNormal(const df1b2vector & x, const df1b2variable & mu, const df1b2variable & sigma)
  {
    df1b2variable nll=.0;  
    double n=double(size_count(x));
    nll=n*log(sigma*sqrt(2.*M_PI)+EPS) +0.5*norm2((x-mu)/(sigma+EPS));  
    return nll;
  }
 


  dvariable nllNormal(const dvar_vector & x, const dvector & mu, const double sigma)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    double n=double(x.size());
    nll=n*log(sigma*sqrt(2.*M_PI)+EPS) +0.5*norm2((x-mu)/(sigma+EPS));  
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }   
  dvariable nllNormal(const dvector & x, const dvar_vector & mu, const dvariable & sigma)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    double n=double(x.size());
    nll=n*log(sigma*sqrt(2.*M_PI)+EPS) +0.5*norm2((x-mu)/(sigma+EPS));  
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllNormal(const dvar_vector & x, const dvar_vector & mu, const dvariable & sigma)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    double n=double(x.size());
    nll=n*log(sigma*sqrt(2.*M_PI)+EPS) +0.5*norm2((x-mu)/(sigma+EPS));  
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }

  df1b2variable nllNormal(const df1b2vector & x, const dvector & mu, const double sigma)
  {
    df1b2variable nll=.0;  
    double n=double(size_count(x));
    nll=n*log(sigma*sqrt(2.*M_PI)+EPS) +0.5*norm2((x-mu)/(sigma+EPS));  
    return nll;
  }   
  df1b2variable nllNormal(const dvector & x, const df1b2vector & mu, const df1b2variable & sigma)
  {
    df1b2variable nll=.0;  
    double n=double(size_count(x));
    nll=n*log(sigma*sqrt(2.*M_PI)+EPS) +0.5*norm2((x-mu)/(sigma+EPS));  
    return nll;
  } 
  df1b2variable nllNormal(const df1b2vector & x, const df1b2vector & mu, const df1b2variable & sigma)
  {
    df1b2variable nll=.0;  
    double n=double(size_count(x));
    nll=n*log(sigma*sqrt(2.*M_PI)+EPS) +0.5*norm2((x-mu)/(sigma+EPS));  
    return nll;
  } 






  //
  //                                    Normal with tau
  //

  dvariable nllNormal2(const dvariable & x, const double mu, const double tau)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    nll = -0.5*log(tau+EPS)+0.5*tau*square(x-mu)+0.5*log(2.*M_PI);//full likelihood, can drop last term
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllNormal2(const double x, const dvariable & mu, const dvariable & tau)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    nll = -0.5*log(tau+EPS)+0.5*tau*square(x-mu)+0.5*log(2.*M_PI);//full likelihood, can drop last term
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllNormal2(const dvariable & x, const dvariable & mu, const dvariable & tau)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    nll = -0.5*log(tau+EPS)+0.5*tau*square(x-mu)+0.5*log(2.*M_PI);//full likelihood, can drop last term
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }

  df1b2variable nllNormal2(const df1b2variable & x, const double mu, const double tau)
  {
    df1b2variable nll=.0;  
    nll = -0.5*log(tau+EPS)+0.5*tau*square(x-mu)+0.5*log(2.*M_PI);//full likelihood, can drop last term
    return nll;
  }
  df1b2variable nllNormal2(const double x, const df1b2variable & mu, const df1b2variable & tau)
  {
    df1b2variable nll=.0;  
    nll = -0.5*log(tau+EPS)+0.5*tau*square(x-mu)+0.5*log(2.*M_PI);//full likelihood, can drop last term
    return nll;
  }
  df1b2variable nllNormal2(const df1b2variable & x, const df1b2variable & mu, const df1b2variable & tau)
  {
    df1b2variable nll=.0;  
    nll = -0.5*log(tau+EPS)+0.5*tau*square(x-mu)+0.5*log(2.*M_PI);//full likelihood, can drop last term
    return nll;
  }




  dvariable nllNormal2(const dvar_vector & x, const double mu, const double tau)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    double n=double(x.size());
    nll = -0.5*n*log(tau+EPS)+0.5*tau*norm2(x-mu)+0.5*n*log(2.*M_PI);//full likelihood, can drop last term
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllNormal2(const dvector & x, const dvariable & mu, const dvariable & tau)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    double n=double(x.size());
    nll = -0.5*n*log(tau+EPS)+0.5*tau*norm2(x-mu)+0.5*n*log(2.*M_PI);//full likelihood, can drop last term
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllNormal2(const dvar_vector & x, const dvariable & mu, const dvariable & tau)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    double n=double(x.size());
    nll = -0.5*n*log(tau+EPS)+0.5*tau*norm2(x-mu)+0.5*n*log(2.*M_PI);//full likelihood, can drop last term
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }

  df1b2variable nllNormal2(const df1b2vector & x, const double  mu, const double tau)
  {
    df1b2variable nll=.0;  
    double n=double(size_count(x));
    nll = -0.5*n*log(tau+EPS)+0.5*tau*norm2(x-mu)+0.5*n*log(2.*M_PI);//full likelihood, can drop last term
    return nll;
  }
  df1b2variable nllNormal2(const dvector & x, const df1b2variable &  mu, const df1b2variable & tau)
  {
    df1b2variable nll=.0;  
    double n=double(size_count(x));
    nll = -0.5*n*log(tau+EPS)+0.5*tau*norm2(x-mu)+0.5*n*log(2.*M_PI);//full likelihood, can drop last term
    return nll;
  }
  df1b2variable nllNormal2(const df1b2vector & x, const df1b2variable &  mu, const df1b2variable & tau)
  {
    df1b2variable nll=.0;  
    double n=double(size_count(x));
    nll = -0.5*n*log(tau+EPS)+0.5*tau*norm2(x-mu)+0.5*n*log(2.*M_PI);//full likelihood, can drop last term
    return nll;
  }



  dvariable nllNormal2(const dvar_vector & x, const dvector & mu, const double tau)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    double n=double(x.size());
    nll = -0.5*n*log(tau+EPS)+0.5*tau*norm2(x-mu)+0.5*n*log(2.*M_PI);
    RETURN_ARRAYS_DECREMENT();
    return nll;
  } 
  dvariable nllNormal2(const dvector & x, const dvar_vector & mu, const dvariable & tau)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    double n=double(x.size());
    nll = -0.5*n*log(tau+EPS)+0.5*tau*norm2(x-mu)+0.5*n*log(2.*M_PI);
    RETURN_ARRAYS_DECREMENT();
    return nll;
  } 
  dvariable nllNormal2(const dvar_vector & x, const dvar_vector & mu, const dvariable & tau)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    double n=double(x.size());
    nll = -0.5*n*log(tau+EPS)+0.5*tau*norm2(x-mu)+0.5*n*log(2.*M_PI);
    RETURN_ARRAYS_DECREMENT();
    return nll;
  } 

  df1b2variable nllNormal2(const df1b2vector & x, const dvector & mu, const double tau)
  {
    df1b2variable nll=.0;  
    double n=double(size_count(x));
    nll = -0.5*n*log(tau+EPS)+0.5*tau*norm2(x-mu)+0.5*n*log(2.*M_PI);
    return nll;
  } 
  df1b2variable nllNormal2(const dvector & x, const df1b2vector & mu, const df1b2variable & tau)
  {
    df1b2variable nll=.0;  
    double n=double(size_count(x));
    nll = -0.5*n*log(tau+EPS)+0.5*tau*norm2(x-mu)+0.5*n*log(2.*M_PI);
    return nll;
  } 
  df1b2variable nllNormal2(const df1b2vector & x, const df1b2vector & mu, const df1b2variable & tau)
  {
    df1b2variable nll=.0;  
    double n=double(size_count(x));
    nll = -0.5*n*log(tau+EPS)+0.5*tau*norm2(x-mu)+0.5*n*log(2.*M_PI);
    return nll;
  } 







  //
  //                                    Binomial
  //

  double nllBinomial(const double x, const double n,const double p)
  {   
    double nll=0;
    nll=-x*log(p+EPS)-(n-x)*log(1.-p+EPS)-log_comb(n,x); //full likelihood, can drop the last term
    return nll;
  }
  dvariable nllBinomial(const dvariable & x, const double n,const double p)
  {   
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;
    nll=-x*log(p+EPS)-(n-x)*log(1.-p+EPS)-log_comb(n,x); //full likelihood, can drop the last term
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllBinomial(const double x, const dvariable & n,const dvariable & p)
  {   
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;
    nll=-x*log(p+EPS)-(n-x)*log(1.-p+EPS)-log_comb(n,x); //full likelihood, can drop the last term
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllBinomial(const dvariable & x, const dvariable & n,const dvariable & p)
  {   
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;
    nll=-x*log(p+EPS)-(n-x)*log(1.-p+EPS)-log_comb(n,x); //full likelihood, can drop the last term
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }

  df1b2variable nllBinomial(const df1b2variable & x, const double n,const double p)
  {   
    df1b2variable nll=.0;
    nll=-x*log(p+EPS)-(n-x)*log(1.-p+EPS)-log_comb(n,x); //full likelihood, can drop the last term
    return nll;
  }
  df1b2variable nllBinomial(const double x, const df1b2variable & n,const df1b2variable & p)
  {   
    df1b2variable nll=.0;
    nll=-x*log(p+EPS)-(n-x)*log(1.-p+EPS)-log_comb(n,x); //full likelihood, can drop the last term
    return nll;
  }
  df1b2variable nllBinomial(const df1b2variable & x, const df1b2variable & n,const df1b2variable & p)
  {   
    df1b2variable nll=.0;
    nll=-x*log(p+EPS)-(n-x)*log(1.-p+EPS)-log_comb(n,x); //full likelihood, can drop the last term
    return nll;
  }




  dvariable nllBinomial(const dvar_vector & x, const dvector & n,const double p)
  {   
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;
    nll=-sum(x)*log(p+EPS)-sum(n-x)*log(1.-p+EPS)-sum(log_comb(n,x)); //full likelihood, can drop the last term
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllBinomial(const dvector & x, const dvar_vector & n,const dvariable & p)
  {   
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;
    nll=-sum(x)*log(p+EPS)-sum(n-x)*log(1.-p+EPS)-sum(log_comb(n,x)); //full likelihood, can drop the last term
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllBinomial(const dvar_vector & x, const dvar_vector & n,const dvariable & p)
  {   
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;
    nll=-sum(x)*log(p+EPS)-sum(n-x)*log(1.-p+EPS)-sum(log_comb(n,x)); //full likelihood, can drop the last term
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }

  df1b2variable nllBinomial(const dvector & x, const df1b2vector & n,const df1b2variable & p)
  {   
    df1b2variable nll=.0;
    nll=-sum(x)*log(p+EPS)-sum(n-x)*log(1.-p+EPS)-sum(log_comb(n,x)); //full likelihood, can drop the last term
    return nll;
  }
  df1b2variable nllBinomial(const df1b2vector & x, const dvector & n,const double p)
  {   
    df1b2variable nll=.0;
    nll=-sum(x)*log(p+EPS)-sum(n-x)*log(1.-p+EPS)-sum(log_comb(n,x)); //full likelihood, can drop the last term
    return nll;
  }
  df1b2variable nllBinomial(const df1b2vector & x, const df1b2vector & n,const df1b2variable & p)
  {   
    df1b2variable nll=.0;
    nll=-sum(x)*log(p+EPS)-sum(n-x)*log(1.-p+EPS)-sum(log_comb(n,x)); //full likelihood, can drop the last term
    return nll;
  }








  //
  //                                    Lognormal
  //

  double nllLognormal(const double x, const double mu, const double sigma)
  { 
    double nll=0;     
    nll=log(sigma*sqrt(2.*M_PI)+EPS)+log(x+EPS)+0.5*square((log(x+EPS)-mu)/(sigma+EPS));//full likelihood
    return nll;
  }
  dvariable nllLognormal(const dvariable& x, const double mu, const double sigma)
  { 
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;    
    nll=log(sigma*sqrt(2.*M_PI)+EPS)+log(x+EPS)+0.5*square((log(x+EPS)-mu)/(sigma+EPS));//full likelihood
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllLognormal(const double x, const dvariable & mu, const dvariable & sigma)
  { 
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;    
    nll=log(sigma*sqrt(2.*M_PI)+EPS)+log(x+EPS)+0.5*square((log(x+EPS)-mu)/(sigma+EPS));//full likelihood
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllLognormal(const dvariable & x, const dvariable & mu, const double sigma)
  { 
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;    
    nll=log(sigma*sqrt(2.*M_PI)+EPS)+log(x+EPS)+0.5*square((log(x+EPS)-mu)/(sigma+EPS));//full likelihood
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllLognormal(const dvariable & x, const double mu, const dvariable & sigma)
  { 
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;    
    nll=log(sigma*sqrt(2.*M_PI)+EPS)+log(x+EPS)+0.5*square((log(x+EPS)-mu)/(sigma+EPS));//full likelihood
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllLognormal(const dvariable & x, const dvariable & mu, const dvariable & sigma)
  { 
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;    
    nll=log(sigma*sqrt(2.*M_PI)+EPS)+log(x+EPS)+0.5*square((log(x+EPS)-mu)/(sigma+EPS));//full likelihood
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }

  df1b2variable nllLognormal(const df1b2variable & x, const double mu, const double sigma)
  { 
    df1b2variable nll=.0;     
    nll=log(sigma*sqrt(2.*M_PI)+EPS)+log(x+EPS)+0.5*square((log(x+EPS)-mu)/(sigma+EPS));//full likelihood
    return nll;
  }
  df1b2variable nllLognormal(const double x, const df1b2variable & mu, const df1b2variable & sigma)
  { 
    df1b2variable nll=.0;     
    nll=log(sigma*sqrt(2.*M_PI)+EPS)+log(x+EPS)+0.5*square((log(x+EPS)-mu)/(sigma+EPS));//full likelihood
    return nll;
  }
  df1b2variable nllLognormal(const df1b2variable & x, const df1b2variable & mu, const double sigma)
  { 
    df1b2variable nll=.0;     
    nll=log(sigma*sqrt(2.*M_PI)+EPS)+log(x+EPS)+0.5*square((log(x+EPS)-mu)/(sigma+EPS));//full likelihood
    return nll;
  }
  df1b2variable nllLognormal(const df1b2variable & x, const double mu, const df1b2variable & sigma)
  { 
    df1b2variable nll=.0;     
    nll=log(sigma*sqrt(2.*M_PI)+EPS)+log(x+EPS)+0.5*square((log(x+EPS)-mu)/(sigma+EPS));//full likelihood
    return nll;
  }
  df1b2variable nllLognormal(const df1b2variable & x, const df1b2variable & mu, const df1b2variable & sigma)
  { 
    df1b2variable nll=.0;     
    nll=log(sigma*sqrt(2.*M_PI)+EPS)+log(x+EPS)+0.5*square((log(x+EPS)-mu)/(sigma+EPS));//full likelihood
    return nll;
  }




  dvariable nllLognormal(const dvar_vector & x, const double mu, const double sigma)
  { 
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    double n=double(x.size());
    nll=n*log(sigma*sqrt(2.*M_PI)+EPS) +sum(log(x+EPS))+0.5*norm2((log(x+EPS)-mu)/(sigma+EPS));
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }  
  dvariable nllLognormal(const dvector & x, const dvariable & mu, const dvariable & sigma)
  { 
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    double n=double(x.size());
    nll=n*log(sigma*sqrt(2.*M_PI)+EPS) +sum(log(x+EPS))+0.5*norm2((log(x+EPS)-mu)/(sigma+EPS));
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }  
  dvariable nllLognormal(const dvar_vector & x, const dvariable & mu, const dvariable & sigma)
  { 
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    double n=double(x.size());
    nll=n*log(sigma*sqrt(2.*M_PI)+EPS) +sum(log(x+EPS))+0.5*norm2((log(x+EPS)-mu)/(sigma+EPS));
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }  

  df1b2variable nllLognormal(const df1b2vector & x, const double mu, const double sigma)
  { 
    df1b2variable nll=.0;   
    double n=double(size_count(x));
    nll=n*log(sigma*sqrt(2.*M_PI)+EPS) +sum(log(x+EPS))+0.5*norm2((log(x+EPS)-mu)/(sigma+EPS));
    return nll;
  }  
  df1b2variable nllLognormal(const dvector & x, const df1b2variable & mu, const df1b2variable & sigma)
  { 
    df1b2variable nll=.0;   
    double n=double(size_count(x));
    nll=n*log(sigma*sqrt(2.*M_PI)+EPS) +sum(log(x+EPS))+0.5*norm2((log(x+EPS)-mu)/(sigma+EPS));
    return nll;
  }  
  df1b2variable nllLognormal(const df1b2vector & x, const df1b2variable & mu, const df1b2variable & sigma)
  { 
    df1b2variable nll=.0;   
    double n=double(size_count(x));
    nll=n*log(sigma*sqrt(2.*M_PI)+EPS) +sum(log(x+EPS))+0.5*norm2((log(x+EPS)-mu)/(sigma+EPS));
    return nll;
  }  



  dvariable nllLognormal(const dvar_vector & x, const dvector & mu, const double sigma)
  { 
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    double n=double(x.size());
    nll=n*log(sigma*sqrt(2.*M_PI)+EPS) +sum(log(x+EPS))+0.5*norm2((log(x+EPS)-mu)/(sigma+EPS));
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }  
  dvariable nllLognormal(const dvector & x, const dvar_vector & mu, const dvariable & sigma)
  { 
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    double n=double(x.size());
    nll=n*log(sigma*sqrt(2.*M_PI)+EPS) +sum(log(x+EPS))+0.5*norm2((log(x+EPS)-mu)/(sigma+EPS));
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }  
  dvariable nllLognormal(const dvar_vector & x, const dvar_vector & mu, const dvariable & sigma)
  { 
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    double n=double(x.size());
    nll=n*log(sigma*sqrt(2.*M_PI)+EPS) +sum(log(x+EPS))+0.5*norm2((log(x+EPS)-mu)/(sigma+EPS));
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }  

  df1b2variable nllLognormal(const df1b2vector & x, const dvector & mu, const double sigma)
  { 
    df1b2variable nll=.0;   
    double n=double(size_count(x));
    nll=n*log(sigma*sqrt(2.*M_PI)+EPS) +sum(log(x+EPS))+0.5*norm2((log(x+EPS)-mu)/(sigma+EPS));
    return nll;
  }  
  df1b2variable nllLognormal(const dvector & x, const df1b2vector & mu, const df1b2variable & sigma)
  { 
    df1b2variable nll=.0;   
    double n=double(size_count(x));
    nll=n*log(sigma*sqrt(2.*M_PI)+EPS) +sum(log(x+EPS))+0.5*norm2((log(x+EPS)-mu)/(sigma+EPS));
    return nll;
  } 
  df1b2variable nllLognormal(const df1b2vector & x, const df1b2vector & mu, const df1b2variable & sigma)
  { 
    df1b2variable nll=.0;   
    double n=double(size_count(x));
    nll=n*log(sigma*sqrt(2.*M_PI)+EPS) +sum(log(x+EPS))+0.5*norm2((log(x+EPS)-mu)/(sigma+EPS));
    return nll;
  } 






  //
  //                                    Lognormal with tau
  //

  dvariable nllLognormal2(const double x, const dvariable & mu, const dvariable & tau)
  { 
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    nll= -0.5*log(tau+EPS)+log(x+EPS)+0.5*tau*square(log(x+EPS)-mu)+0.5*log(2.*M_PI);//full likelihood
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllLognormal2(const dvariable & x, const double mu, const double tau)
  { 
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    nll= -0.5*log(tau+EPS)+log(x+EPS)+0.5*tau*square(log(x+EPS)-mu)+0.5*log(2.*M_PI);//full likelihood
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllLognormal2(const dvariable & x, const dvariable & mu, const dvariable & tau)
  { 
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    nll= -0.5*log(tau+EPS)+log(x+EPS)+0.5*tau*square(log(x+EPS)-mu)+0.5*log(2.*M_PI);//full likelihood
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }

  df1b2variable nllLognormal2(const df1b2variable x, const double mu, const double tau)
  { 
    df1b2variable nll=.0;   
    nll= -0.5*log(tau+EPS)+log(x+EPS)+0.5*tau*square(log(x+EPS)-mu)+0.5*log(2.*M_PI);//full likelihood
    return nll;
  }
  df1b2variable nllLognormal2(const double x, const df1b2variable & mu, const df1b2variable & tau)
  { 
    df1b2variable nll=.0;   
    nll= -0.5*log(tau+EPS)+log(x+EPS)+0.5*tau*square(log(x+EPS)-mu)+0.5*log(2.*M_PI);//full likelihood
    return nll;
  }
  df1b2variable nllLognormal2(const df1b2variable x, const df1b2variable & mu, const double tau)
  { 
    df1b2variable nll=.0;   
    nll= -0.5*log(tau+EPS)+log(x+EPS)+0.5*tau*square(log(x+EPS)-mu)+0.5*log(2.*M_PI);//full likelihood
    return nll;
  }
  df1b2variable nllLognormal2(const df1b2variable x, const double mu, const df1b2variable & tau)
  { 
    df1b2variable nll=.0;   
    nll= -0.5*log(tau+EPS)+log(x+EPS)+0.5*tau*square(log(x+EPS)-mu)+0.5*log(2.*M_PI);//full likelihood
    return nll;
  }
  df1b2variable nllLognormal2(const df1b2variable x, const df1b2variable & mu, const df1b2variable & tau)
  { 
    df1b2variable nll=.0;   
    nll= -0.5*log(tau+EPS)+log(x+EPS)+0.5*tau*square(log(x+EPS)-mu)+0.5*log(2.*M_PI);//full likelihood
    return nll;
  }




  dvariable nllLognormal2(const dvar_vector & x, const double mu, const double tau)
  { 
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;
    double n=double(x.size());
    nll=-0.5*n*log(tau+EPS)+sum(log(x+EPS))+0.5*tau*norm2(log(x+EPS)-mu)+0.5*n*log(2.*M_PI);
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllLognormal2(const dvector & x, const dvariable & mu, const dvariable & tau)
  { 
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;
    double n=double(x.size());
    nll=-0.5*n*log(tau+EPS)+sum(log(x+EPS))+0.5*tau*norm2(log(x+EPS)-mu)+0.5*n*log(2.*M_PI);
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllLognormal2(const dvar_vector & x, const dvariable & mu, const dvariable & tau)
  { 
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;
    double n=double(x.size());
    nll=-0.5*n*log(tau+EPS)+sum(log(x+EPS))+0.5*tau*norm2(log(x+EPS)-mu)+0.5*n*log(2.*M_PI);
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }

  df1b2variable nllLognormal2(const dvector & x, const df1b2variable & mu, const df1b2variable & tau)
  { 
    df1b2variable nll=.0;
    double n=double(size_count(x)); 
    nll=-0.5*n*log(tau+EPS)+sum(log(x+EPS))+0.5*tau*norm2(log(x+EPS)-mu)+0.5*n*log(2.*M_PI);
    return nll;
  }
  df1b2variable nllLognormal2(const df1b2vector & x, const double mu, const double tau)
  { 
    df1b2variable nll=.0;
    double n=double(size_count(x)); 
    nll=-0.5*n*log(tau+EPS)+sum(log(x+EPS))+0.5*tau*norm2(log(x+EPS)-mu)+0.5*n*log(2.*M_PI);
    return nll;
  }
  df1b2variable nllLognormal2(const df1b2vector & x, const df1b2variable & mu, const double tau)
  { 
    df1b2variable nll=.0;
    double n=double(size_count(x)); 
    nll=-0.5*n*log(tau+EPS)+sum(log(x+EPS))+0.5*tau*norm2(log(x+EPS)-mu)+0.5*n*log(2.*M_PI);
    return nll;
  }
  df1b2variable nllLognormal2(const df1b2vector & x, const double mu, const df1b2variable & tau)
  { 
    df1b2variable nll=.0;
    double n=double(size_count(x)); 
    nll=-0.5*n*log(tau+EPS)+sum(log(x+EPS))+0.5*tau*norm2(log(x+EPS)-mu)+0.5*n*log(2.*M_PI);
    return nll;
  }
  df1b2variable nllLognormal2(const df1b2vector & x, const df1b2variable & mu, const df1b2variable & tau)
  { 
    df1b2variable nll=.0;
    double n=double(size_count(x)); 
    nll=-0.5*n*log(tau+EPS)+sum(log(x+EPS))+0.5*tau*norm2(log(x+EPS)-mu)+0.5*n*log(2.*M_PI);
    return nll;
  }




  dvariable nllLognormal2(const dvar_vector & x, const dvector & mu, const double tau)
  { 
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    double n=double(x.size());
    nll=-0.5*n*log(tau+EPS)+sum(log(x+EPS))+0.5*tau*norm2(log(x+EPS)-mu)+0.5*n*log(2.*M_PI);
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }  
  dvariable nllLognormal2(const dvector & x, const dvar_vector & mu, const dvariable & tau)
  { 
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    double n=double(x.size());
    nll=-0.5*n*log(tau+EPS)+sum(log(x+EPS))+0.5*tau*norm2(log(x+EPS)-mu)+0.5*n*log(2.*M_PI);
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }  
  dvariable nllLognormal2(const dvar_vector & x, const dvar_vector & mu, const dvariable & tau)
  { 
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    double n=double(x.size());
    nll=-0.5*n*log(tau+EPS)+sum(log(x+EPS))+0.5*tau*norm2(log(x+EPS)-mu)+0.5*n*log(2.*M_PI);
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }  

  df1b2variable nllLognormal2(const df1b2vector & x, const dvector & mu, const double tau)
  { 
    df1b2variable nll=.0;   
    double n=double(size_count(x));
    nll=-0.5*n*log(tau+EPS)+sum(log(x+EPS))+0.5*tau*norm2(log(x+EPS)-mu)+0.5*n*log(2.*M_PI);
    return nll;
  }  
  df1b2variable nllLognormal2(const dvector & x, const df1b2vector & mu, const df1b2variable & tau)
  { 
    df1b2variable nll=.0;   
    double n=double(size_count(x));
    nll=-0.5*n*log(tau+EPS)+sum(log(x+EPS))+0.5*tau*norm2(log(x+EPS)-mu)+0.5*n*log(2.*M_PI);
    return nll;
  } 
  df1b2variable nllLognormal2(const df1b2vector & x, const df1b2vector & mu, const df1b2variable & tau)
  { 
    df1b2variable nll=.0;   
    double n=double(size_count(x));
    nll=-0.5*n*log(tau+EPS)+sum(log(x+EPS))+0.5*tau*norm2(log(x+EPS)-mu)+0.5*n*log(2.*M_PI);
    return nll;
  } 







  //
  //                                    Poisson
  //

  double nllPoisson(const double x, const double lambda)
  {
    double nll=0;
    nll = lambda- x*log(lambda+EPS)+gammln(x+1.0);//full likelihood
    return nll;
  }
  dvariable nllPoisson(const dvariable & x, const double lambda)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;
    nll = lambda- x*log(lambda+EPS)+gammln(x+1.0);//full likelihood
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllPoisson(const double x, const dvariable & lambda)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;
    nll = lambda- x*log(lambda+EPS)+gammln(x+1.0);//full likelihood
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllPoisson(const dvariable & x, const dvariable & lambda)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;
    nll = lambda- x*log(lambda+EPS)+gammln(x+1.0);//full likelihood
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }

  df1b2variable nllPoisson(const df1b2variable & x, const double lambda)
  {
    df1b2variable nll=.0;
    nll = lambda- x*log(lambda+EPS)+gammln(x+1.0);//full likelihood
    return nll;
  }
  df1b2variable nllPoisson(const double x, const df1b2variable & lambda)
  {
    df1b2variable nll=.0;
    nll = lambda- x*log(lambda+EPS)+gammln(x+1.0);//full likelihood
    return nll;
  }
  df1b2variable nllPoisson(const df1b2variable & x, const df1b2variable & lambda)
  {
    df1b2variable nll=.0;
    nll = lambda- x*log(lambda+EPS)+gammln(x+1.0);//full likelihood
    return nll;
  }



  dvariable nllPoisson(const dvar_vector & x, const double lambda)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;
    double n=double(x.size());
    nll = n*lambda- sum(x)*log(lambda+EPS)+sum(gammln(x+1.0));//full likelihood, can drop the last term
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }  
  dvariable nllPoisson(const dvector & x, const dvariable & lambda)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;
    double n=double(x.size());
    nll = n*lambda- sum(x)*log(lambda+EPS)+sum(gammln(x+1.0));//full likelihood, can drop the last term
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }  
  dvariable nllPoisson(const dvar_vector & x, const dvariable & lambda)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;
    double n=double(x.size());
    nll = n*lambda- sum(x)*log(lambda+EPS)+sum(gammln(x+1.0));//full likelihood, can drop the last term
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }  

  df1b2variable nllPoisson(const dvector & x, const df1b2variable & lambda)
  {
    df1b2variable nll=.0;
    double n=double(size_count(x)); 
    nll = n*lambda- sum(x)*log(lambda+EPS)+sum(gammln(x+1.0));//full likelihood, can drop the last term
    return nll;
  }  
  df1b2variable nllPoisson(const df1b2vector & x, const double lambda)
  {
    df1b2variable nll=.0;
    double n=double(size_count(x)); 
    nll = n*lambda- sum(x)*log(lambda+EPS)+sum(gammln(x+1.0));//full likelihood, can drop the last term
    return nll;
  }  
  df1b2variable nllPoisson(const df1b2vector & x, const df1b2variable & lambda)
  {
    df1b2variable nll=.0;
    double n=double(size_count(x)); 
    nll = n*lambda- sum(x)*log(lambda+EPS)+sum(gammln(x+1.0));//full likelihood, can drop the last term
    return nll;
  }  






  //
  //                                    Inverse Gamma
  //

  double nllInverseGamma(const double x, const double a, const double b)
  {
    double nll=0;  
    nll= - a*log(b)+gammln(a)  +(a+1.)*log(x+EPS) + b/(x+EPS);
    return nll;
  }
  dvariable nllInverseGamma(const dvariable & x, const double a, const double b)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    nll= - a*log(b)+gammln(a)  +(a+1.)*log(x+EPS) + b/(x+EPS);
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllInverseGamma(const double x, const dvariable & a, const dvariable & b)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    nll= - a*log(b)+gammln(a)  +(a+1.)*log(x+EPS) + b/(x+EPS);
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllInverseGamma(const dvariable & x, const dvariable & a, const dvariable & b)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    nll= - a*log(b)+gammln(a)  +(a+1.)*log(x+EPS) + b/(x+EPS);
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }

  df1b2variable nllInverseGamma(const df1b2variable & x, const double a, const double b)
  {
    df1b2variable nll=.0;  
    nll= - a*log(b)+gammln(a)  +(a+1.)*log(x+EPS) + b/(x+EPS);
    return nll;
  }
  df1b2variable nllInverseGamma(const double x, const df1b2variable & a, const df1b2variable & b)
  {
    df1b2variable nll=.0;  
    nll= - a*log(b)+gammln(a)  +(a+1.)*log(x+EPS) + b/(x+EPS);
    return nll;
  }
  df1b2variable nllInverseGamma(const df1b2variable & x, const df1b2variable & a, const double b)
  {
    df1b2variable nll=.0;  
    nll= - a*log(b)+gammln(a)  +(a+1.)*log(x+EPS) + b/(x+EPS);
    return nll;
  }
  df1b2variable nllInverseGamma(const df1b2variable & x, const double a, const df1b2variable & b)
  {
    df1b2variable nll=.0;  
    nll= - a*log(b)+gammln(a)  +(a+1.)*log(x+EPS) + b/(x+EPS);
    return nll;
  }
  df1b2variable nllInverseGamma(const df1b2variable & x, const df1b2variable & a, const df1b2variable & b)
  {
    df1b2variable nll=.0;  
    nll= - a*log(b)+gammln(a)  +(a+1.)*log(x+EPS) + b/(x+EPS);
    return nll;
  }




  dvariable nllInverseGamma(const dvar_vector & x, const double a, const double b)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    double n=double(x.size());
    nll= - n*a*log(b)+ n*gammln(a)  +(a+1.)*sum(log(x+EPS)) + b/sum(x+EPS);
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllInverseGamma(const dvector & x, const dvariable & a, const dvariable & b)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    double n=double(x.size());
    nll= - n*a*log(b)+ n*gammln(a)  +(a+1.)*sum(log(x+EPS)) + b/sum(x+EPS);
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }
  dvariable nllInverseGamma(const dvar_vector & x, const dvariable & a, const dvariable & b)
  {
    RETURN_ARRAYS_INCREMENT();
    dvariable nll=0;  
    double n=double(x.size());
    nll= - n*a*log(b)+ n*gammln(a)  +(a+1.)*sum(log(x+EPS)) + b/sum(x+EPS);
    RETURN_ARRAYS_DECREMENT();
    return nll;
  }

  df1b2variable nllInverseGamma(const dvector & x, const df1b2variable & a, const df1b2variable & b)
  {
    df1b2variable nll=.0;  
    double n=double(size_count(x)); 
    nll= - n*a*log(b)+ n*gammln(a)  +(a+1.)*sum(log(x+EPS)) + b/sum(x+EPS);
    return nll;
  }
  df1b2variable nllInverseGamma(const df1b2vector & x, const double a, const double b)
  {
    df1b2variable nll=.0;  
    double n=double(size_count(x)); 
    nll= - n*a*log(b)+ n*gammln(a)  +(a+1.)*sum(log(x+EPS)) + b/sum(x+EPS);
    return nll;
  }
  df1b2variable nllInverseGamma(const df1b2vector & x, const df1b2variable & a, const df1b2variable & b)
  {
    df1b2variable nll=.0;  
    double n=double(size_count(x)); 
    nll= - n*a*log(b)+ n*gammln(a)  +(a+1.)*sum(log(x+EPS)) + b/sum(x+EPS);
    return nll;
  }