distributions_R.hpp

Go to the documentation of this file.

// Copyright (C) 2013-2015 Kasper Kristensen
// License: GPL-2

template<class Type>
Type pnorm(Type q, Type mean = 0., Type sd = 1.){
  CppAD::vector<Type> tx(1);
  tx[0] = (q - mean) / sd;
  return atomic::pnorm1(tx)[0];
}
VECTORIZE3_ttt(pnorm)
VECTORIZE1_t(pnorm)


template<class Type>
Type qnorm(Type p, Type mean = 0., Type sd = 1.){
  CppAD::vector<Type> tx(1);
  tx[0] = p;
  return sd*atomic::qnorm1(tx)[0] + mean;
}
VECTORIZE3_ttt(qnorm)
VECTORIZE1_t(qnorm)


template<class Type>
Type pgamma(Type q, Type shape, Type scale = 1.){
  CppAD::vector<Type> tx(4);
  tx[0] = q/scale;
  tx[1] = shape;
  tx[2] = Type(0);        // 0'order deriv
  tx[3] = -lgamma(shape); // normalize
  return atomic::D_incpl_gamma_shape(tx)[0];
}
VECTORIZE3_ttt(pgamma)


template<class Type>
Type qgamma(Type q, Type shape, Type scale = 1.){
  CppAD::vector<Type> tx(3);
  tx[0] = q;
  tx[1] = shape;
  tx[2] = -lgamma(shape); // normalize
  return atomic::inv_incpl_gamma(tx)[0] * scale;
}
VECTORIZE3_ttt(qgamma)


template<class Type>
Type ppois(Type q, Type lambda){
  CppAD::vector<Type> tx(2);
  tx[0] = q;
  tx[1] = lambda;
  return atomic::ppois(tx)[0];
}
VECTORIZE2_tt(ppois)


template<class Type> 
Type pexp(Type x, Type rate)
{
        return CppAD::CondExpGe(x,Type(0),1-exp(-rate*x),Type(0));
}

// Vectorize pexp
VECTORIZE2_tt(pexp)


template<class Type> 
Type dexp(Type x, Type rate, int give_log=0)
{
        if(!give_log)
                return CppAD::CondExpGe(x,Type(0),rate*exp(-rate*x),Type(0));
        else
                return CppAD::CondExpGe(x,Type(0),log(rate)-rate*x,Type(-INFINITY));
}

// Vectorize dexp
VECTORIZE3_tti(dexp)


template <class Type>
Type qexp(Type p, Type rate)
{
        return -log(1-p)/rate;
}

// Vectorize qexp.
VECTORIZE2_tt(qexp)
template<class Type> 
Type pweibull(Type x, Type shape, Type scale)
{
        return CppAD::CondExpGe(x,Type(0),1-exp(-pow(x/scale,shape)),Type(0));
}

// Vectorize pweibull
VECTORIZE3_ttt(pweibull)


template<class Type> 
Type dweibull(Type x, Type shape, Type scale, int give_log=0)
{
        if(!give_log)
                return CppAD::CondExpGe(x,Type(0),shape/scale * pow(x/scale,shape-1) * exp(-pow(x/scale,shape)),Type(0));
        else
                return CppAD::CondExpGe(x,Type(0),log(shape) - log(scale) + (shape-1)*(log(x)-log(scale)) - pow(x/scale,shape),Type(-INFINITY));
}

// Vectorize dweibull
VECTORIZE4_ttti(dweibull)


template<class Type> 
Type qweibull(Type p, Type shape, Type scale)
{
        Type res = scale * pow( (-log(1-p)) , 1/shape );
        res = CppAD::CondExpLt(p,Type(0),Type(0),res);
        res = CppAD::CondExpGt(p,Type(1),Type(0),res);
        return res;
}

// Vectorize qweibull
VECTORIZE3_ttt(qweibull)
template<class Type> 
Type dbinom(Type k, Type size, Type prob, int give_log=0)
{
  Type logres = lgamma(size + 1) - lgamma(k + 1) - lgamma(size - k + 1);
  // Add 'k * log(prob)' only if k > 0
  logres += CppAD::CondExpGt(k, Type(0), k * log(prob), Type(0) );
  // Add '(size - k) * log(1 - prob)' only if size > k
  logres += CppAD::CondExpGt(size, k, (size - k) * log(1 - prob), Type(0) );
  if (!give_log) return exp(logres);
  else return logres;
}

// Vectorize dbinom
VECTORIZE4_ttti(dbinom)


template<class Type>
Type dbinom_robust(Type k, Type size, Type logit_p, int give_log=0)
{
  CppAD::vector<Type> tx(4);
  tx[0] = k;
  tx[1] = size;
  tx[2] = logit_p;
  tx[3] = 0;
  Type ans = atomic::log_dbinom_robust(tx)[0]; /* without norm. constant */
  if (size > 1) {
    ans += lgamma(size+1.) - lgamma(k+1.) - lgamma(size-k+1.);
  }
  return ( give_log ? ans : exp(ans) );
}
VECTORIZE4_ttti(dbinom_robust)


template <class Type>
Type dbeta(Type x, Type shape1, Type shape2, int give_log)
{
        Type res = exp(lgamma(shape1+shape2) - lgamma(shape1) - lgamma(shape2)) * pow(x,shape1-1) * pow(1-x,shape2-1);
        if(!give_log) 
                return res;
        else 
                return CppAD::CondExpEq(x,Type(0),log(res),lgamma(shape1+shape2) - lgamma(shape1) - lgamma(shape2) + (shape1-1)*log(x) + (shape2-1)*log(1-x));
}

// Vectorize dbeta
VECTORIZE4_ttti(dbeta)


template <class Type>
Type df(Type x, Type df1, Type df2, int give_log)
{
        Type logres = lgamma((df1+df2)/2.) - lgamma(df1/2.) - lgamma(df2/2.) + df1/2.*log(Type(df1)/df2) + (df1/2.-1)*log(x) - (df1+df2)/2.*log(1+Type(df1)/df2*x);
        if(!give_log) return exp(logres);
        else return logres;
}

//Vectorize df
VECTORIZE4_ttti(df)


template <class Type>
Type dlogis(Type x, Type location, Type scale, int give_log)
{
        Type logres = -(x-location)/scale - log(scale) - 2*log(1+exp(-(x-location)/scale));
        if(!give_log) return exp(logres);
        else return logres;
}

// Vectorize dlogis
VECTORIZE4_ttti(dlogis)


template <class Type>
Type dsn(Type x, Type alpha, int give_log=0)
{
        
        if(!give_log) return 2 * dnorm(x,Type(0),Type(1),0) * pnorm(alpha*x);
        else return log(2.0) + log(dnorm(x,Type(0),Type(1),0)) + log(pnorm(alpha*x));
}

// Vectorize dsn
VECTORIZE3_tti(dsn)


template <class Type>
Type dt(Type x, Type df, int give_log)
{
        Type logres = lgamma((df+1)/2) - Type(1)/2*log(df*M_PI) -lgamma(df/2) - (df+1)/2*log(1+x*x/df);
        if(!give_log) return exp(logres);
        else return logres;
}

// Vectorize dt
VECTORIZE3_tti(dt)


template <class Type>
Type dmultinom(vector<Type> x, vector<Type> p, int give_log=0)
{
        vector<Type> xp1 = x+Type(1);
        Type logres = lgamma(x.sum() + Type(1)) - lgamma(xp1).sum() + (x*log(p)).sum();
        if(give_log) return logres;
        else return exp(logres);
}

template <class Type>
Type dSHASHo(Type x, Type mu, Type sigma, Type nu, Type tau, int give_log = 0)
{
        // TODO : Replace log(x+sqrt(x^2+1)) by a better approximation for asinh(x).
                
        Type z = (x-mu)/sigma;
        Type c = cosh(tau*log(z+sqrt(z*z+1))-nu);
        Type r = sinh(tau*log(z+sqrt(z*z+1))-nu);
        Type logres = -log(sigma) + log(tau) -0.5*log(2*M_PI) -0.5*log(1+(z*z)) +log(c) -0.5*(r*r);
                                                
        if(!give_log) return exp(logres);
        else return logres;
}

// Vectorize dSHASHo
VECTORIZE6_ttttti(dSHASHo)


template <class Type>
Type pSHASHo(Type q,Type mu,Type sigma,Type nu,Type tau,int give_log=0)
{
        // TODO : Replace log(x+sqrt(x^2+1)) by a better approximation for asinh(x).

        Type z = (q-mu)/sigma;
        Type r = sinh(tau * log(z+sqrt(z*z+1)) - nu);
        Type p = pnorm(r);
                                        
        if (!give_log) return p;
        else return log(p);
}

// Vectorize pSHASHo
VECTORIZE6_ttttti(pSHASHo)


template <class Type>
Type qSHASHo(Type p, Type mu, Type sigma, Type nu, Type tau, int log_p = 0)
{
        // TODO : Replace log(x+sqrt(x^2+1)) by a better approximation for asinh(x).

        if(!log_p) return mu + sigma*sinh((1/tau)* log(qnorm(p)+sqrt(qnorm(p)*qnorm(p)+1)) + (nu/tau));
        else return mu + sigma*sinh((1/tau)*log(qnorm(exp(p))+sqrt(qnorm(exp(p))*qnorm(exp(p))+1))+(nu/tau));
}

// Vectorize qSHASHo
VECTORIZE6_ttttti(qSHASHo)


template <class Type>
Type norm2SHASHo(Type x, Type mu, Type sigma, Type nu, Type tau, int log_p = 0)
{

        return qSHASHo(pnorm(x),mu,sigma,nu,tau,log_p);
}

// Vectorize norm2SHASHo
VECTORIZE6_ttttti(norm2SHASHo)
template<class Type>
Type pbeta(Type q, Type shape1, Type shape2){
  CppAD::vector<Type> tx(4);
  tx[0] = q;
  tx[1] = shape1;
  tx[2] = shape2;
  tx[3] = 0; // order
  Type ans = atomic::pbeta(tx)[0];
  return ans;
}
VECTORIZE3_ttt(pbeta)


template<class Type>
Type qbeta(Type p, Type shape1, Type shape2){
  CppAD::vector<Type> tx(3);
  tx[0] = p;
  tx[1] = shape1;
  tx[2] = shape2;
  Type ans = atomic::qbeta(tx)[0];
  return ans;
}
VECTORIZE3_ttt(qbeta)


template<class Type>
Type besselK(Type x, Type nu){
  Type ans;
  if(CppAD::Variable(nu)) {
    CppAD::vector<Type> tx(3);
    tx[0] = x;
    tx[1] = nu;
    tx[2] = 0;
    ans = atomic::bessel_k(tx)[0];
  } else {
    CppAD::vector<Type> tx(2);
    tx[0] = x;
    tx[1] = nu;
    ans = atomic::bessel_k_10(tx)[0];
  }
  return ans;
}
VECTORIZE2_tt(besselK)


template<class Type>
Type besselI(Type x, Type nu){
  Type ans;
  if(CppAD::Variable(nu)) {
    CppAD::vector<Type> tx(3);
    tx[0] = x;
    tx[1] = nu;
    tx[2] = 0;
    ans = atomic::bessel_i(tx)[0];
  } else {
    CppAD::vector<Type> tx(2);
    tx[0] = x;
    tx[1] = nu;
    ans = atomic::bessel_i_10(tx)[0];
  }
  return ans;
}
VECTORIZE2_tt(besselI)


template<class Type>
Type besselJ(Type x, Type nu){
  CppAD::vector<Type> tx(3);
  tx[0] = x;
  tx[1] = nu;
  tx[2] = 0;
  Type ans = atomic::bessel_j(tx)[0];
  return ans;
}
VECTORIZE2_tt(besselJ)


template<class Type>
Type besselY(Type x, Type nu){
  CppAD::vector<Type> tx(3);
  tx[0] = x;
  tx[1] = nu;
  tx[2] = 0;
  Type ans = atomic::bessel_y(tx)[0];
  return ans;
}
VECTORIZE2_tt(besselY)


template<class Type>
Type dtweedie(Type y, Type mu, Type phi, Type p, int give_log = 0) {
  Type p1 = p - 1.0, p2 = 2.0 - p;
  Type ans = -pow(mu, p2) / (phi * p2); // log(prob(y=0))
  if (y > 0 || CppAD::Variable(y)) {
    CppAD::vector<Type> tx(4);
    tx[0] = y;
    tx[1] = phi;
    tx[2] = p;
    tx[3] = 0;
    Type ans2 = atomic::tweedie_logW(tx)[0];
    ans2 += -y / (phi * p1 * pow(mu, p1)) - log(y);
    if (CppAD::Variable(y)) {
      ans += CppAD::CondExpGt(y, Type(0), ans2, Type(0));
    } else {
      ans += ans2;
    }
  }
  return ( give_log ? ans : exp(ans) );
}
VECTORIZE5_tttti(dtweedie)



template<class Type>
Type compois_calc_logZ(Type loglambda, Type nu) {
  CppAD::vector<Type> tx(3);
  tx[0] = loglambda;
  tx[1] = nu;
  tx[2] = 0;
  return atomic::compois_calc_logZ(tx)[0];
}
VECTORIZE2_tt(compois_calc_logZ)


template<class Type>
Type compois_calc_loglambda(Type logmean, Type nu) {
  CppAD::vector<Type> tx(3);
  tx[0] = logmean;
  tx[1] = nu;
  tx[2] = 0;
  return atomic::compois_calc_loglambda(tx)[0];
}
VECTORIZE2_tt(compois_calc_loglambda)


template<class T1, class T2, class T3>
T1 dcompois(T1 x, T2 mode, T3 nu, int give_log = 0) {
  T2 loglambda = nu * log(mode);
  T1 ans = x * loglambda - nu * lfactorial(x);
  ans -= compois_calc_logZ(loglambda, nu);
  return ( give_log ? ans : exp(ans) );
}

template<class T1, class T2, class T3>
T1 dcompois2(T1 x, T2 mean, T3 nu, int give_log = 0) {
  T2 logmean = log(mean);
  T2 loglambda = compois_calc_loglambda(logmean, nu);
  T1 ans = x * loglambda - nu * lfactorial(x);
  ans -= compois_calc_logZ(loglambda, nu);
  return ( give_log ? ans : exp(ans) );
}

/********************************************************************/
/* SIMULATON CODE                                                   */
/********************************************************************/

extern "C" {
  double Rf_rnorm(double mu, double sigma);
}
template<class Type>
Type rnorm(Type mu, Type sigma)
{
  return Rf_rnorm(asDouble(mu), asDouble(sigma));
}
VECTORIZE2_tt(rnorm)
VECTORIZE2_n(rnorm)

extern "C" {
  double Rf_rpois(double mu);
}
template<class Type>
Type rpois(Type mu)
{
  return Rf_rpois(asDouble(mu));
}
VECTORIZE1_t(rpois)
VECTORIZE1_n(rpois)

extern "C" {
  double Rf_runif(double a, double b);
}
template<class Type>
Type runif(Type a, Type b)
{
  return Rf_runif(asDouble(a), asDouble(b));
}
VECTORIZE2_tt(runif)
VECTORIZE2_n(runif)

extern "C" {
  double Rf_rbinom(double size, double prob);
}
template<class Type>
Type rbinom(Type size, Type prob)
{
  return Rf_rbinom(asDouble(size), asDouble(prob));
}
VECTORIZE2_tt(rbinom)
VECTORIZE2_n(rbinom)

extern "C" {
  double Rf_rgamma(double shape, double scale);
}
template<class Type>
Type rgamma(Type shape, Type scale)
{
  return Rf_rgamma(asDouble(shape), asDouble(scale));
}
VECTORIZE2_tt(rgamma)
VECTORIZE2_n(rgamma)

extern "C" {
  double Rf_rexp(double rate);
}
template<class Type>
Type rexp(Type rate)
{
  return Rf_rexp(asDouble(rate));
}

VECTORIZE1_t(rexp)
VECTORIZE1_n(rexp)

extern "C" {
        double Rf_rbeta(double shape1, double shape2);
}
template<class Type>
Type rbeta(Type shape1, Type shape2)
{
        return Rf_rbeta(asDouble(shape1), asDouble(shape2));
}

VECTORIZE2_tt(rbeta)
VECTORIZE2_n(rbeta)

extern "C" {
        double Rf_rf(double df1, double df2);
}
template<class Type>
Type rf(Type df1, Type df2)
{
        return Rf_rf(asDouble(df1), asDouble(df2));
}

VECTORIZE2_tt(rf)
VECTORIZE2_n(rf)

extern "C" {
        double Rf_rlogis(double location, double scale);
}
template<class Type>
Type rlogis(Type location, Type scale)
{
        return Rf_rlogis(asDouble(location), asDouble(scale));
}

VECTORIZE2_tt(rlogis)
VECTORIZE2_n(rlogis)

extern "C" {
        double Rf_rt(double df);
}
template<class Type>
Type rt(Type df)
{
        return Rf_rt(asDouble(df));
}

VECTORIZE1_t(rt)
VECTORIZE1_n(rt)

extern "C" {
        double Rf_rweibull(double shape, double scale);
}
template<class Type>
Type rweibull(Type shape, Type scale)
{
        return Rf_rweibull(asDouble(shape), asDouble(scale));
}

VECTORIZE2_tt(rweibull)
VECTORIZE2_n(rweibull)


template<class Type>
Type rcompois(Type mode, Type nu)
{
  Type loglambda = nu * log(mode);
  return atomic::compois_utils::simulate(asDouble(loglambda), asDouble(nu));
}
VECTORIZE2_tt(rcompois)
VECTORIZE2_n(rcompois)


template<class Type>
Type rcompois2(Type mean, Type nu)
{
  Type logmean = log(mean);
  Type loglambda = compois_calc_loglambda(logmean, nu);
  return atomic::compois_utils::simulate(asDouble(loglambda), asDouble(nu));
}
VECTORIZE2_tt(rcompois2)

// Note: Vectorize manually to avoid many identical calls to
// 'calc_loglambda'.
template<class Type>
vector<Type> rcompois2(int n, Type mean, Type nu)
{
  Type logmean = log(mean);
  Type loglambda = compois_calc_loglambda(logmean, nu);
  Type mode = exp(loglambda / nu);
  return rcompois(n, mode, nu);
}

template<class Type>
Type rtweedie(Type mu, Type phi, Type p) {
  // Copied from R function tweedie::rtweedie
  Type lambda = pow(mu, 2. - p) / (phi * (2. - p));
  Type alpha  = (2. - p) / (1. - p);
  Type gam = phi * (p - 1.) * pow(mu, p - 1.);
  Type N = rpois(lambda);
  Type ans = rgamma( -alpha * N /* shape */, gam /* scale */);
  return ans;
}
VECTORIZE3_ttt(rtweedie)
VECTORIZE3_n(rtweedie)