# Example R code for inverse cdf and rejection sampling.
#
# M.V.Wickerhauser, 2011/03/02
# M.V.Wickerhauser, 2024/03/22



/***************************************************************
 * Maxima code to find the inverse cdf for exp(-2x) on [2,inf) *
 ***************************************************************/
/* pdf without the constant */
f:exp(-2*x);
/* integrate to find the constant */
c:integrate(f,x,2,inf);
/* normalize the pdf */
f:(1/c)*f;
/* ...which gives f(x) = 2*exp(4-2*x)  */
/* integrate to find the cdf F */
F:integrate(f,x,2,x);
/* solve y=F(x) for x to find the inverse of cdf F */
x:solve(y=f,x);
/* This gives x=2+log(2/y)/2 */

##############################################################
# R code for inverse cdf sampler for exp(-2x) on [2,inf)
#
isamp<-function(n=5000, seed=NULL) {
 if( !is.null(seed) ) set.seed(seed);  # in case a seed is provided
 y<- runif(n);     # step 1: sample Unif([0,1])
 samples<- 2+log(2/y)/2 # step 2: apply inverse cdf
 hist(samples, breaks=20, prob=TRUE);
 x<- seq(2,10,by=0.01);  # some points in [2,infinity)
 fx<- 2*exp(4-2*x);  # target pdf       
 lines(x,fx);
}
# Run experiments
isamp(n=1000)
isamp(n=2000,seed=63130)



/*************************************************** 
 * Maxima code to find the max value of beta(2,2): *
 ***************************************************/
/* pdf without the constant */
f:x*(1-x);
/* integrate to find the constant */
c:integrate(f,x,0,1);
/* normalize the pdf */
f:(1/c)*f;
/* differentiate to get f' */
diff(f,x);
/* solve f'(x0)=0 to find max */
x0:solve(0=diff(f,x),x);
/* evaluate f(x0) to get max value */
subst(x0,f);




##############################################################
# From Seefeld and Linder, p.189  --- Rejection sampling, beta(2,2)
rsamp<-function(n=5000, seed=NULL) {
 if( !is.null(seed) ) set.seed(seed);  # in case a seed is provided
 x<-runif(n);  # step 1 
 u<-runif(n);  # step 2
 fx<- dbeta(x,2,2);  # step 3
 acceptx<-x[fx>1.5*u]; # step 4   # scale by 1.5 = max(dbeta)
 hist(acceptx, breaks=20, prob=TRUE);
 points(x,fx);
 length(acceptx);   # announce how many samples were accepted
}
#
# Run experiments
rsamp(n=1000)
rsamp(n=2000,seed=63130)