InvStablePrior Package

Description

Contains random number generation, plotting, and estimation algorithms for performing Bayesian inference on a parameter of some widely-used exponential data models. The package contains algorithms for the inverse stable-Poisson, inverse stable-exponential, inverse stable-double exponential, inverse stable-inverse gamma, and inverse stable-half-normal models.

Details

References:

Cahoy and Sedransk (2019)<doi:10.1007/s42519-018-0027-2>

Meerschaert and Straka (2013)<doi:10.1051/mmnp/20138201>

Mainardi, Mura, and Pagnini (2010) <doi:10.1155/2010/104505>

Author(s)

Dexter Cahoy cahoyd@uhd.edu

Bayesian inference for the true rate of double exponential (Laplace) distribution.

Description

Generates random numbers from the prior and posterior distributions of the inverse stable-double exponential model. The random variates can be used to simulate prior and posterior predictive distributions as well.

Usage

isdexp(x, B, alp, rho)

Arguments

Value

list consisting of the vectors of random numbers from the prior and posterior distributions, the accepted sample size, and the acceptance probability of the adaptive rejection sampling procedure (Algorithm 2 of the first reference below).

References

Cahoy and Sedransk (2019). Inverse stable prior for exponential models. Journal of Statistical Theory and Practice, 13, Article 29. <doi:10.1007/s42519-018-0027-2>

Meerschaert and Straka (2013). Inverse stable subordinators. Math. Model. Nat. Phenom., 8(2), 1-16. <doi:10.1051/mmnp/20138201>

Mainardi, Mura, and Pagnini (2010). The M-Wright Function in Time-Fractional Diffusion Processes: A Tutorial Survey. Int. J. Differ. Equ., Volume 2010. <doi:10.1155/2010/104505>

Examples

alp=0.95
require(nimble)
dat=rdexp(30, location = 0,  rate = 2)
rho=1/sd(dat)


#b=n
#a=sum(abs(dat) )
#rho=optimize(function(r){exp(-b)*(b/a)^b - (r^b)*exp(-a*r)}, c(0,20),  tol=10^(-50)  )$min

out= isdexp(dat, B=1000000, alp , rho)
#prior samples
thetprior=unlist(out[2])
summary(thetprior)

#posterior samples
thet=unlist(out[1])

#95% Credible intervals
quantile (thet, c(0.025,0.975) )
summary(thet)

#The accepted sample size:
unlist(out[3])

#The acceptance probability:
unlist(out[4])

#Plotting with normalization to have a maximum of 1
#for comparing prior and posterior
out2=density(thet)
ymaxpost=max(out2$y)
out3=density(thetprior)
ymaxprior=max(out3$y)
plot(out2$x,out2$y/ymaxpost, xlim=c(0,5), col="blue", type="l",
 xlab="theta", ylab="density", lwd=2, frame.plot=FALSE)
lines(out3$x,out3$y/ymaxprior,lwd=2,col="red")
#points(dat,rep(0,length(dat)), pch='*')


#Generating 1000 random numbers from the Inverse Stable (alpha=0.4,rho=5) prior
U1 = runif(1000)
U2 = runif(1000)
alp=0.4
rho=5
stab = ( ( sin(alp*pi*U1)*sin((1-alp)*pi*U1)^(1/alp-1) )
/ (  ( sin(pi*U1)^(1/alp) )*abs(log(U2))^(1/alp-1))  )
#Inverse stable random numbers are below:
#rho*stab^(-alp)

Bayesian inference for the true rate of exponential distribution.

Description

Generates random numbers from the prior and posterior distributions of the inverse stable-exponential model. The random variates can be used to simulate prior and posterior predictive distributions as well.

Usage

isexp(x, B, alp, rho)

Arguments

Value

list consisting of the vectors of random numbers from the prior and posterior distributions, the accepted sample size, and the acceptance probability of the adaptive rejection sampling procedure (Algorithm 2 of the first reference below).

References

Cahoy and Sedransk (2019). Inverse stable prior for exponential models. Journal of Statistical Theory and Practice, 13, Article 29. <doi:10.1007/s42519-018-0027-2>

Meerschaert and Straka (2013). Inverse stable subordinators. Math. Model. Nat. Phenom., 8(2), 1-16. <doi:10.1051/mmnp/20138201>

Mainardi, Mura, and Pagnini (2010). The M-Wright Function in Time-Fractional Diffusion Processes: A Tutorial Survey. Int. J. Differ. Equ., Volume 2010. <doi:10.1155/2010/104505>

Examples

alp=0.5
dat=rexp(10,rate=0.5)
rho=1/mean(dat)
#rho=1/mean(dat) + 3*sd(dat)
#rho=1/mean(dat) - 3*sd(dat)

#b=length(dat)
#a=sum(dat)
#rho=optimize(function(r){exp(-b)*(b/a)^b - (r^b)*exp(-a*r)}, c(0,20),  tol=10^(-50)  )$min

out= isexp(dat, B=1000000, alp , rho)
#prior samples
thetprior=unlist(out[2])
summary(thetprior)

#posterior samples
thet=unlist(out[1])

#95% Credible intervals
quantile (thet, c(0.025,0.975) )
summary(thet)

#The accepted sample size:
unlist(out[3])

#The acceptance probability:
unlist(out[4])

#Plotting with normalization to have a maximum of 1
#for comparing prior and posterior
out2=density(thet)
ymaxpost=max(out2$y)
out3=density(thetprior)
ymaxprior=max(out3$y)
plot(out2$x,out2$y/ymaxpost, xlim=c(0,1), col="blue", type="l",
 xlab="theta", ylab="density",lwd=2,  frame.plot=FALSE)
lines(out3$x,out3$y/ymaxprior,lwd=2, col="red")
#points(dat,rep(0,length(dat)), pch='*')


#Generating 1000 random numbers from the Inverse Stable (alpha=0.4,rho=5) prior
U1 = runif(1000)
U2 = runif(1000)
alp=0.4
rho=5
stab = ( ( sin(alp*pi*U1)*sin((1-alp)*pi*U1)^(1/alp-1) )
/ (  ( sin(pi*U1)^(1/alp) )*abs(log(U2))^(1/alp-1))  )
#Inverse stable random numbers are below:
#rho*stab^(-alp)

Bayesian inference for the true inverse mean/rate of half-normal/half-Gaussian distribution.

Description

Generates random numbers from the prior and posterior distributions of the inverse stable-half-normal model. The random variates can be used to simulate prior and posterior predictive distributions as well.

Usage

ishalfn(x, B, alp, rho)

Arguments

Value

list consisting of the vectors of random numbers from the prior and posterior distributions, the accepted sample size, and the acceptance probability of the adaptive rejection sampling procedure (Algorithm 2 of the first reference below).

References

Cahoy and Sedransk (2019). Inverse stable prior for exponential models. Journal of Statistical Theory and Practice, 13, Article 29. <doi:10.1007/s42519-018-0027-2>

Meerschaert and Straka (2013). Inverse stable subordinators. Math. Model. Nat. Phenom., 8(2), 1-16. <doi:10.1051/mmnp/20138201>

Mainardi, Mura, and Pagnini (2010). The M-Wright Function in Time-Fractional Diffusion Processes: A Tutorial Survey. Int. J. Differ. Equ., Volume 2010. <doi:10.1155/2010/104505>

Examples

alp=0.95
require(fdrtool)
dat=rhalfnorm(100, theta=sqrt(pi/2) )
rho=1/mean(dat)


#b=length(dat)/2
#a=sum(dat^2)/pi
#rho=optimize(function(r){exp(-b)*(b/a)^b - (r^b)*exp(-a*r)}, c(0,20),  tol=10^(-50)  )$min

out= ishalfn(dat, B=1000000, alp , rho)
#prior samples
thetprior=unlist(out[2])
summary(thetprior)

#posterior samples
thet=unlist(out[1])

#95% Credible intervals
quantile (thet, c(0.025,0.975) )
summary(thet)

#The accepted sample size:
unlist(out[3])

#The acceptance probability:
unlist(out[4])

#Plotting with normalization to have a maximum of 1
#for comparing prior and posterior
out2=density(thet)
ymaxpost=max(out2$y)
out3=density(thetprior)
ymaxprior=max(out3$y)
plot(out2$x,out2$y/ymaxpost, xlim=c(0,2), col="blue", type="l",
 xlab="theta", ylab="density", lwd=2,  frame.plot=FALSE)
lines(out3$x,out3$y/ymaxprior, lwd=2,col="red")
#points(dat,rep(0,length(dat)), pch='*')


#Generating 1000 random numbers from the Inverse Stable (alpha=0.4,rho=5) prior
U1 = runif(1000)
U2 = runif(1000)
alp=0.4
rho=5
stab = ( ( sin(alp*pi*U1)*sin((1-alp)*pi*U1)^(1/alp-1) )
/ (  ( sin(pi*U1)^(1/alp) )*abs(log(U2))^(1/alp-1))  )
#Inverse stable random numbers are below:
#rho*stab^(-alp)

Bayesian inference for the true scale of inverse gamma distribution.

Description

Generates random numbers from the prior and posterior distributions of the inverse stable-inverse gamma model. The random variates can be used to simulate prior and posterior predictive distributions as well.

Usage

isinvgam(x, B, alp, rho, sh)

Arguments

Value

list consisting of the vectors of random numbers from the prior and posterior distributions, the accepted sample size, and the acceptance probability of the adaptive rejection sampling procedure (Algorithm 2 of the first reference below).

References

Cahoy and Sedransk (2019). Inverse stable prior for exponential models. Journal of Statistical Theory and Practice, 13, Article 29. <doi:10.1007/s42519-018-0027-2>

Meerschaert and Straka (2013). Inverse stable subordinators. Math. Model. Nat. Phenom., 8(2), 1-16. <doi:10.1051/mmnp/20138201>

Mainardi, Mura, and Pagnini (2010). The M-Wright Function in Time-Fractional Diffusion Processes: A Tutorial Survey. Int. J. Differ. Equ., Volume 2010. <doi:10.1155/2010/104505>

Examples

alp=0.95
require(nimble)
sh=2.1 # a>2 so variance exists, known
dat=rinvgamma(50, shape=sh,  scale = 4)
rho= (sh-1)*mean(dat)


#b=n
#a=sum(1/dat )
#rho=optimize(function(r){exp(-b)*(b/a)^b - (r^b)*exp(-a*r)}, c(0,20),  tol=10^(-50)  )$min

out= isinvgam(dat, B=1000000, alp , rho,sh)
#prior samples
thetprior=unlist(out[2])
summary(thetprior)

#posterior samples
thet=unlist(out[1])
summary(thet)

#95% Credible intervals
quantile (thet, c(0.025,0.975) )
summary(thet)

#The accepted sample size:
unlist(out[3])

#The acceptance probability:
unlist(out[4])

#Plotting with normalization to have a maximum of 1
#for comparing prior and posterior
out2=density(thet)
ymaxpost=max(out2$y)
out3=density(thetprior)
ymaxprior=max(out3$y)
plot(out2$x,out2$y/ymaxpost, xlim=c(0,5), col="blue", type="l",
 xlab="theta", ylab="density",lwd=2,  frame.plot=FALSE)
lines(out3$x,out3$y/ymaxprior,lwd=2,col="red")
#points(dat,rep(0,length(dat)), pch='*')


#Generating 1000 random numbers from the Inverse Stable (alpha=0.4,rho=5) prior
U1 = runif(1000)
U2 = runif(1000)
alp=0.4
rho=5
stab = ( ( sin(alp*pi*U1)*sin((1-alp)*pi*U1)^(1/alp-1) )
/ (  ( sin(pi*U1)^(1/alp) )*abs(log(U2))^(1/alp-1))  )
#Inverse stable random numbers are below:
#rho*stab^(-alp)

Bayesian inference for the true mean of Poisson distribution.

Description

Generates random numbers from the prior and posterior distributions of the inverse stable-Poisson model. The random variates can be used to simulate prior and posterior predictive distributions as well.

Usage

ispoi(x, B, alp, rho)

Arguments

Value

list consisting of the vectors of random numbers from the prior and posterior distributions, the accepted sample size, and the acceptance probability of the adaptive rejection sampling procedure (Algorithm 2 of the first reference below).

References

Cahoy and Sedransk (2019). Inverse stable prior for exponential models. Journal of Statistical Theory and Practice, 13, Article 29. <doi:10.1007/s42519-018-0027-2>

Meerschaert and Straka (2013). Inverse stable subordinators. Math. Model. Nat. Phenom., 8(2), 1-16. <doi:10.1051/mmnp/20138201>

Mainardi, Mura, and Pagnini (2010). The M-Wright Function in Time-Fractional Diffusion Processes: A Tutorial Survey. Int. J. Differ. Equ., Volume 2010. <doi:10.1155/2010/104505>

Examples

alp=0.9
dat=rpois(50,lambda=10)
rho=mean(dat)
#rho=mean(dat) + 3*sd(dat)
#rho=mean(dat) - 3*sd(dat)

#a=length(dat)
#b=sum(dat)
#rho=optimize(function(r){exp(-b)*(b/a)^b - (r^b)*exp(-a*r)}, c(0,20),  tol=10^(-50)  )$min

out= ispoi(dat, B=1000000, alp , rho)
#prior samples
thetprior=unlist(out[2])
summary(thetprior)

#posterior samples
thet=unlist(out[1])

#95% Credible intervals
quantile (thet, c(0.025,0.975) )
summary(thet)

#The accepted sample size:
unlist(out[3])

#The acceptance probability:
unlist(out[4])

#Plotting with normalization to have a maximum of 1
#for comparing prior and posterior
out2=density(thet)
ymaxpost=max(out2$y)
out3=density(thetprior)
ymaxprior=max(out3$y)
plot(out2$x,out2$y/ymaxpost, xlim=c(0,15), col="blue", type="l",
 xlab="theta", ylab="density", lwd=2, frame.plot=FALSE)
lines(out3$x,out3$y/ymaxprior,lwd=2,col="red")
#points(dat,rep(0,length(dat)), pch='*')


#Generating 1000 random numbers from the Inverse Stable (alpha=0.4,rho=5) prior
U1 = runif(1000)
U2 = runif(1000)
alp=0.4
rho=5
stab = ( ( sin(alp*pi*U1)*sin((1-alp)*pi*U1)^(1/alp-1) )
/ (  ( sin(pi*U1)^(1/alp) )*abs(log(U2))^(1/alp-1))  )
#Inverse stable random numbers are below:
#rho*stab^(-alp)