| Title: | Adaptive Sampling Algorithms |
|---|---|
| Description: | For distributions whose probability density functions are log-concave, the adaptive rejection sampling algorithm can be used to build envelope functions for sampling. For others, we can use the modified adaptive rejection sampling algorithm, the concave-convex adaptive rejection sampling algorithm and the adaptive slice sampling algorithm. So we designed an R package mainly including 4 functions: rARS(), rMARS(), rCCARS() and rASS(). These functions can realize sampling based on the algorithms above. |
| Authors: | Dong Zhang |
| Maintainer: | Dong Zhang <[email protected]> |
| License: | GPL-2 |
| Version: | 1.1.1 |
| Built: | 2025-01-20 04:25:04 UTC |
| Source: | https://github.com/cran/AdapSamp |
rARS generates a sequence of random numbers using the adaptive rejection sampling algorithm.
rARS(n, formula, min = -Inf, max = Inf, sp)rARS(n, formula, min = -Inf, max = Inf, sp)
n |
Desired sample size; |
formula |
Kernal of the target density; |
min, max
|
Domain including positive and negative infinity of the target distribution; |
sp |
Supporting set. |
Dong Zhang <[email protected]>
# Example 1: Standard normal distribution x1 <- rARS(100,"exp(-x^2/2)",-Inf,Inf,c(-2,2)) # Example 2: Truncated normal distribution x2 <- rARS(100,"exp(-x^2/2)",-2.1,2.1,c(-2,2)) # Example 3: Normal distribution with mean=2 and sd=2 x3 <- rARS(100,"exp(-(x-2)^2/(2*4))",-Inf,Inf,c(-3,3)) # Example 4: Exponential distribution with rate=3 x4 <- rARS(100,"exp(-3*x)",0,Inf,c(2,3,100)) # Example 5: Beta distribution with alpha=3 and beta=4 x5 <- rARS(100,"x^2*(1-x)^3",0,1,c(0.4,0.6)) # Example 6: Gamma distribution with alpha=5 and lambda=2 x6 <- rARS(100,"x^(5-1)*exp(-2*x)",0,Inf,c(1,10)) # Example 7: Student distribution with df=10 x7 <- rARS(100,"(1+x^2/10)^(-(10+1)/2)",-Inf,Inf,c(-10,2)) # Example 8: F distribution with m=10 and n=5 x8 <- rARS(100,"x^(10/2-1)/(1+10/5*x)^(15/2)",0,Inf,c(3,10)) # Example 9:Cauchy distribution x9 <- rARS(100,"1/(1+(x-1)^2)",-Inf,Inf,c(-2,2,10)) # Example 10:Rayleigh distribution with lambda=1 x10 <- rARS(100,"2*x*exp(-x^2)",0,Inf,c(0.01,10))# Example 1: Standard normal distribution x1 <- rARS(100,"exp(-x^2/2)",-Inf,Inf,c(-2,2)) # Example 2: Truncated normal distribution x2 <- rARS(100,"exp(-x^2/2)",-2.1,2.1,c(-2,2)) # Example 3: Normal distribution with mean=2 and sd=2 x3 <- rARS(100,"exp(-(x-2)^2/(2*4))",-Inf,Inf,c(-3,3)) # Example 4: Exponential distribution with rate=3 x4 <- rARS(100,"exp(-3*x)",0,Inf,c(2,3,100)) # Example 5: Beta distribution with alpha=3 and beta=4 x5 <- rARS(100,"x^2*(1-x)^3",0,1,c(0.4,0.6)) # Example 6: Gamma distribution with alpha=5 and lambda=2 x6 <- rARS(100,"x^(5-1)*exp(-2*x)",0,Inf,c(1,10)) # Example 7: Student distribution with df=10 x7 <- rARS(100,"(1+x^2/10)^(-(10+1)/2)",-Inf,Inf,c(-10,2)) # Example 8: F distribution with m=10 and n=5 x8 <- rARS(100,"x^(10/2-1)/(1+10/5*x)^(15/2)",0,Inf,c(3,10)) # Example 9:Cauchy distribution x9 <- rARS(100,"1/(1+(x-1)^2)",-Inf,Inf,c(-2,2,10)) # Example 10:Rayleigh distribution with lambda=1 x10 <- rARS(100,"2*x*exp(-x^2)",0,Inf,c(0.01,10))
rASS generates a sequence of random numbers by the adaptive slice sampling algorithm with stepping-out procedures.
rASS(n, x0 = 0, formula, w = 3)rASS(n, x0 = 0, formula, w = 3)
n |
Desired sample size; |
x0 |
Initial value; |
formula |
Target density function p(x); |
w |
Length of the coverage interval. |
Dong Zhang <[email protected]>
Neal R M. Slice sampling - Rejoinder[J]. Annals of Statistics, 2003, 31(3):758-767.
# Example 1: Sampling from exponential distribution with bounded domain x<-rASS(100,-1,"1.114283*exp(-(4-x^2)^2)",3) plot(density(x))# Example 1: Sampling from exponential distribution with bounded domain x<-rASS(100,-1,"1.114283*exp(-(4-x^2)^2)",3) plot(density(x))
rCCARS generates a sequence of random numbers by the concave-convex adaptive rejection sampling algorithm from target distributions with bounded domain.
rCCARS(n, cvformula, ccformula, min, max, sp)rCCARS(n, cvformula, ccformula, min, max, sp)
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 |
Strictly speaking, the concave-convex adaptive rejection sampling algorithm can generate samples from target distributions who have bounded domains. For distributions with unbounded domain, rCCARS can also be used for sampling approximately. For example, if we want draw a sequence from N(0,1) by the concave-convex adaptive rejection sampling algorithm. We know that X~N(0,1) has a so small probability in two tails that we can ingore the parts at both ends. Pr(X>20)=P(X<-20)=2.753624e-89, therefore we can get random numbers approximately from N(0,1) with the bound [-20,20]. Also, you can make this bound large enough to reduce sampling error.
Dong Zhang <[email protected]>
Teh Y W. Concave-Convex Adaptive Rejection Sampling[J]. Journal of Computational & Graphical Statistics, 2011, 20(3):670-691.
# Example 1: 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") #The following examples are also available; #But it may take a few minutes to run them. # Example 2: Expontional bounded distribution # x<-rCCARS(1000,"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) # Example 3: Makeham bounded distribution # x<-rCCARS(1000,"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))} # lines(seq(0,5,0.01),f(seq(0,5,0.01)),col="red",lty=3,lwd=2,type="l")# Example 1: 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") #The following examples are also available; #But it may take a few minutes to run them. # Example 2: Expontional bounded distribution # x<-rCCARS(1000,"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) # Example 3: Makeham bounded distribution # x<-rCCARS(1000,"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))} # lines(seq(0,5,0.01),f(seq(0,5,0.01)),col="red",lty=3,lwd=2,type="l")
rMARS generates a sequence of random numbers using the modified adaptive rejection sampling algorithm.
rMARS(n, formula, min = -Inf, max = Inf, sp, infp, m = 10^(-4))rMARS(n, formula, min = -Inf, max = Inf, sp, infp, m = 10^(-4))
n |
Desired sample size; |
formula |
Kernel of the target distribution; |
min, max
|
Domain including positive and negative infinity of the target distribution; |
sp |
Supporting set; |
infp |
Inflexion set; |
m |
A parameter for judging concavity and convexity in a certain interval. |
Dong Zhang <[email protected]>
Martino L, Miguez J. A generalization of the adaptive rejection sampling algorithm[J]. Statistics & Computing, 2011, 21(4):633-647.
# Example 1: Exponential distribution x <- rMARS(100,"exp(-(4-x^2)^2)",-Inf,Inf, c(-2.5,0,2.5),c(-2/sqrt(3),2/sqrt(3))) hist(x,probability=TRUE,xlim=c(-3,3),ylim=c(0,1.2),breaks=20) lines(density(x,bw=0.05),col="blue") f <- function(x)(exp(-(4-x^2)^2)) lines(seq(-3,3,0.01),f(seq(-3,3,0.01))/integrate(f,-3,3)[[1]],lwd=2,lty=2,col="red") #The following examples are also available; #But it may take a few minutes to run them. # Example 2: Distribution with bounded domain # x <- rMARS(1000,"exp(-(x^2-x^3))",-3,2,c(-1,1),1/3) # hist(x,probability=TRUE,xlim=c(-3,2.5),ylim=c(0,1.2),breaks=20) # lines(density(x,bw=0.2),col="blue") # f <- function(x) exp(-(x^2-x^3)) # lines(seq(-3,2,0.01),f(seq(-3,2,0.01))/integrate(f,-3,2)[[1]],lwd=2,lty=2,col="red",type="l") # Example 3: Weibull distribution with k=3 and lambda=1 # x <- rMARS(100,"3*x^2*exp(-x^3)",10^-15,Inf,c(0.01,1),(1/3)^(1/3),m=10^-4) # hist(x,probability=TRUE,breaks=20,xlim=c(0,2)) # lines(density(x,bw=0.15),col="blue") # f <- function(x) 3*x^2*exp(-x^3) # lines(seq(0,2,0.01),f(seq(0,2,0.01)),lwd=2,lty=2,col="red",type="l") # Example 4: Mixed normal distribution with p=0.3,m1=2,m2=8,sigma1=1,sigma2=2 # x <- rMARS(100,"0.3/sqrt(2*pi)*exp(-(x-2)^2/2)+(1-0.3)/sqrt(2*pi)/2*exp(-(x-8)^2/8)",-Inf,Inf, # c(-6,-4,0,3,6,15),c(-5.120801,-3.357761,3.357761,5.120801),m=10^-8) # hist(x,breaks=20,probability=TRUE);lines(density(x,bw=0.45),col="blue",lwd=2) # f <- function(x)0.3/sqrt(2*pi)*exp(-(x-2)^2/2)+(1-0.3)/sqrt(2*pi)/2*exp(-(x-8)^2/8) # lines(seq(0,14,0.01),f(seq(0,14,0.01)),lty=3,col="red",lwd=2 )# Example 1: Exponential distribution x <- rMARS(100,"exp(-(4-x^2)^2)",-Inf,Inf, c(-2.5,0,2.5),c(-2/sqrt(3),2/sqrt(3))) hist(x,probability=TRUE,xlim=c(-3,3),ylim=c(0,1.2),breaks=20) lines(density(x,bw=0.05),col="blue") f <- function(x)(exp(-(4-x^2)^2)) lines(seq(-3,3,0.01),f(seq(-3,3,0.01))/integrate(f,-3,3)[[1]],lwd=2,lty=2,col="red") #The following examples are also available; #But it may take a few minutes to run them. # Example 2: Distribution with bounded domain # x <- rMARS(1000,"exp(-(x^2-x^3))",-3,2,c(-1,1),1/3) # hist(x,probability=TRUE,xlim=c(-3,2.5),ylim=c(0,1.2),breaks=20) # lines(density(x,bw=0.2),col="blue") # f <- function(x) exp(-(x^2-x^3)) # lines(seq(-3,2,0.01),f(seq(-3,2,0.01))/integrate(f,-3,2)[[1]],lwd=2,lty=2,col="red",type="l") # Example 3: Weibull distribution with k=3 and lambda=1 # x <- rMARS(100,"3*x^2*exp(-x^3)",10^-15,Inf,c(0.01,1),(1/3)^(1/3),m=10^-4) # hist(x,probability=TRUE,breaks=20,xlim=c(0,2)) # lines(density(x,bw=0.15),col="blue") # f <- function(x) 3*x^2*exp(-x^3) # lines(seq(0,2,0.01),f(seq(0,2,0.01)),lwd=2,lty=2,col="red",type="l") # Example 4: Mixed normal distribution with p=0.3,m1=2,m2=8,sigma1=1,sigma2=2 # x <- rMARS(100,"0.3/sqrt(2*pi)*exp(-(x-2)^2/2)+(1-0.3)/sqrt(2*pi)/2*exp(-(x-8)^2/8)",-Inf,Inf, # c(-6,-4,0,3,6,15),c(-5.120801,-3.357761,3.357761,5.120801),m=10^-8) # hist(x,breaks=20,probability=TRUE);lines(density(x,bw=0.45),col="blue",lwd=2) # f <- function(x)0.3/sqrt(2*pi)*exp(-(x-2)^2/2)+(1-0.3)/sqrt(2*pi)/2*exp(-(x-8)^2/8) # lines(seq(0,14,0.01),f(seq(0,14,0.01)),lty=3,col="red",lwd=2 )