```
library(gamlss.dist)
PIG <- function(y, mu, tau) {
a <- sqrt(1 / tau^2 + 2*mu / tau)
(2*a/pi)^0.5 * mu^y*exp(1/tau)*besselK(a,y-0.5)/((a*tau)^y*factorial(y))
}
PIG.pgf <- function(y, mu, tau) {
if(y==0)
return(exp(1/tau*(1-(1+2*tau*mu)^0.5)));
if(y==1)
return(mu*(1+2*tau*mu)^-0.5*PIG.pgf(0,mu,tau));
return (2*tau*mu/(1+2*tau*mu) *
(1-3/(2*y)) * PIG.pgf(y-1,mu,tau) +
mu^2 / (1+2*tau*mu) / (y*(y-1)) * PIG.pgf(y-2,mu,tau));
}
PIG.pgf2 <- function(y, mu, tau) {
exp(tau^-1*(1-(1-2*tau*mu*(y-1))^0.5))
}
gamlss.dist::dPIG(r, exp(mup[i]), exp(mup[i])*beta+tau)
real PIG(int y, real mu, real tau);
real PIG(int y, real mu, real tau) {
if(y==0)
return exp(1/tau*(1-(1+2*tau*mu)^0.5));
if(y==1)
return mu*(1+2*tau*mu)^-0.5*PIG(0,mu,tau);
return 2.*tau*mu/(1+2*tau*mu) *
(1. - 3. /(2.*y)) * PIG(y-1,mu,tau) +
mu^2 / (1+2.*tau*mu) / (y*(y-1)) * PIG(y-2,mu,tau);
}
beta ~ student_t(4, 0, 1);
tau ~ student_t(4, 0, 1);
target += log(PIG(y, exp(mup[i]), exp(mup[i])*beta+tau));
LinseyPoisson<- function(x,theta,a) gamma(x+a)*theta^(a+1)*(a+(x+a)/(1+theta))/(factorial(x)*gamma(a+1)*(1+theta)^(1+x+a))
LinseyPoisson_Log <- function(x,theta,a) (1 + a)*log(theta) - (1 + a + x)*log(1 + theta) + log(a + (a + x)/(1 + theta)) -
lgamma(x+1) - lgamma(1 + a) + lgamma(a + x)
la = exp(mup);
alpha_inv ~ cauchy(0, 1);
real<lower=0> alpha = 1.0 / alpha_inv;
real LinseyPoisson_log(int x, real theta, real a) {
return (1 + a)*log(theta) - (1 + a + x)*log(1 + theta) + log(a + (a + x)/(1 + theta)) -
lgamma(x+1) - lgamma(1 + a) + lgamma(a + x);
}
target += LinseyPoisson_log(y, 0.5/la[i]*(alpha+la[i]*(-1+sqrt((4*la[i]+square(alpha+la[i]))/square(la[i])))), alpha);
// Poisson Transmuted Exponential
//#https://arxiv.org/pdf/1504.01097.pdf
pow<-function(x,k) ((x)^(k))
square<-function(x) ((x)**2)
#PTED <- function(x,theta,alpha) theta*((1-alpha)/pow(1+theta,x+1)+2*alpha/pow(1+2*theta,x+1))
PTED <- function(x,theta,alpha) exp( log(theta) + log( (1-alpha) / pow(1 + theta, x+1) + 2 * alpha / pow(1 + 2*theta,x+1)))
vector<lower=-1,upper=1>[2] alpha;
real PTED_log(int x, real theta, real alpha) {
return log(theta) + log( (1-alpha) / pow(1 + theta, x+1) + 2 * alpha / pow(1 + 2*theta,x+1));
}
target += PTED_log(y[1], (2.0 - alpha[1]) / (2.0 * exp(mup[i, 1])), alpha[1]);
# shanker
https://arxiv.org/pdf/1504.01097.pdf
shanker <-function(x,theta)
square(theta)*pow(1+theta,-2-x)*(1+x+theta+square(theta))/(1+square(theta))
real shanker_log(int x, real theta) {
return 2*log(theta)-(2+x)*log(1+theta)-log(1+square(theta))+log(1+x+theta+square(theta));
}
real la1 = exp(mup[i, 1]);
real theta1 =(2 + (2*(1 - 3*pow(la1,2)))/pow(1 + (45*pow(la1,2))/2. + (3*sqrt(3)*sqrt(pow(la1,2)*(8 + 71*pow(la1,2) + 4*pow(la1,4))))/2.,0.3333333333333333) +
pow(2,0.6666666666666666)*pow(2 + 45*pow(la1,2) + 3*sqrt(3)*sqrt(pow(la1,2)*(8 + 71*pow(la1,2) + 4*pow(la1,4))),0.3333333333333333))/(6. *la1) ;
target += shanker_log(y, theta1);
```

Hereâ€™s how I verified the Shanker-Poisson:

Sorry I not have enough time to write a lot of documentation to it.