rCCARS: Concave-Convex Adaptive Rejection Sampling Algorithm rCCARS...

In teacher7738/AdapSamp: Adaptive Sampling Algorithms

Description

Concave-Convex Adaptive Rejection Sampling Algorithm rCCARS generates a sequence by Concave- Convex Adaptive Rejection Sampling algorithm from target densities with bounded domain. For unbounded domain's rCCARS can also be used for sampling approximately, especially for densities whose probability is too small in tails.

Usage

rCCARS(n, cvformula, ccformula, min, max, sp)

Arguments

n	Desired sample size
cvformula, ccformula	Convex and concave decompositions for -ln(p(x)) where p(x) is the kernal of target density
min, max	Domain except positive and negative infinity
sp	Supporting set

Details

Strictly speaking,Concave-Convex Adaptive Rejection Sampling Algorithm can generate samples from target densities who have bounded domains. However, if the target distribution has a tiny probability in the tails, we can also use rCCARS for sampling approximately. For example, if we want to get a sequence drawn from N(0,1), and we know that

References

Teh Y W. Concave-Convex Adaptive Rejection Sampling[J]. Journal of Computational & Graphical Statistics, 2011, 20(3):670-691.

Examples

#Generalized inverse bounded gaussian distribution with lambda=-1 and a=b=2
x<-rCCARS(100,"x+x^-1","2*log(x)",0.001,100,1)
hist(x,breaks=20,probability = TRUE);lines(density(x,bw=0.1),col="red",lwd=2,lty=2)
f <- function(x) {x^(-2)*exp(-x-x^(-1))/0.2797318}
lines(seq(0,5,0.01),f(seq(0,5,0.01)),lwd=2,lty=3,col="blue")

#Expontional bounded distribution
x<-rCCARS(100,"x^4","-8*x^2+16",-3,4,c(-2,1))
hist(x,breaks=30,probability=TRUE);lines(density(x,bw=0.05),col="blue",lwd=2,lty=2)
f <- function(x) exp(-(x^2-4)^2)/ 0.8974381
lines(seq(-3,4,0.01),f(seq(-3,4,0.01)),col="red",lty=3,lwd=2)

#Makeham bounded distribution
x<-rCCARS(100,"x+1/log(2)*(2^x-1)","-log(1+2^x)",0,5,c(1,2,3))
hist(x,breaks=30,probability=TRUE);lines(density(x,bw=0.05),col="blue",lwd=2,lty=2)
f <- function(x){(1+2^x)*exp(-x-1/log(2)*(2^x-1))}
plot(seq(0,5,0.01),f(seq(0,5,0.01)),col="red",lty=3,lwd=2,type="l")

teacher7738/AdapSamp documentation built on May 4, 2019, 4:22 p.m.