R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.1.0  (2005-04-18), ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for a HTML browser interface to help.
Type 'q()' to quit R.

> ### * <HEADER>
> ###
> attach(NULL, name = "CheckExEnv")
> assign(".CheckExEnv", as.environment(2), pos = length(search())) # base
> ## add some hooks to label plot pages for base and grid graphics
> setHook("plot.new", ".newplot.hook")
> setHook("persp", ".newplot.hook")
> setHook("grid.newpage", ".gridplot.hook")
> 
> assign("cleanEx",
+        function(env = .GlobalEnv) {
+ 	   rm(list = ls(envir = env, all.names = TRUE), envir = env)
+            RNGkind("default", "default")
+ 	   set.seed(1)
+    	   options(warn = 1)
+ 	   delayedAssign("T", stop("T used instead of TRUE"),
+ 		  assign.env = .CheckExEnv)
+ 	   delayedAssign("F", stop("F used instead of FALSE"),
+ 		  assign.env = .CheckExEnv)
+ 	   sch <- search()
+ 	   newitems <- sch[! sch %in% .oldSearch]
+ 	   for(item in rev(newitems))
+                eval(substitute(detach(item), list(item=item)))
+ 	   missitems <- .oldSearch[! .oldSearch %in% sch]
+ 	   if(length(missitems))
+ 	       warning("items ", paste(missitems, collapse=", "),
+ 		       " have been removed from the search path")
+        },
+        env = .CheckExEnv)
> assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now
> assign("ptime", proc.time(), env = .CheckExEnv)
> grDevices::postscript("HI-Examples.ps")
> assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv)
> options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly"))
> options(warn = 1)    
> library('HI')
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "arms"
> 
> ### * arms
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: arms
> ### Title: Function to perform Adaptive Rejection Metropolis Sampling
> ### Aliases: arms
> ### Keywords: distribution multivariate misc
> 
> ### ** Examples
> 
> #### ==> Warning: running the examples may take a few minutes! <== ####    
> ### Univariate densities
> ## Unif(-r,r) 
> y <- arms(runif(1,-1,1), function(x,r) 1, function(x,r) (x>-r)*(x<r), 5000, r=2)
> summary(y); hist(y,prob=TRUE,main="Unif(-r,r); r=2")
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-1.99800 -1.02500 -0.02835 -0.01868  1.02500  2.00000 
> ## Normal(mean,1)
> norldens <- function(x,mean) -(x-mean)^2/2 
> y <- arms(runif(1,3,17), norldens, function(x,mean) ((x-mean)>-7)*((x-mean)<7),
+           5000, mean=10)
> summary(y); hist(y,prob=TRUE,main="Gaussian(m,1); m=10")
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  6.474   9.295   9.981   9.986  10.680  13.710 
> curve(dnorm(x,mean=10),3,17,add=TRUE)
> ## Exponential(1)
> y <- arms(5, function(x) -x, function(x) (x>0)*(x<70), 5000)
> summary(y); hist(y,prob=TRUE,main="Exponential(1)")
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
0.0001741 0.2793000 0.6996000 0.9766000 1.3660000 6.9470000 
> curve(exp(-x),0,8,add=TRUE)
> ## Gamma(4.5,1) 
> y <- arms(runif(1,1e-4,20), function(x) 3.5*log(x)-x,
+           function(x) (x>1e-4)*(x<20), 5000)
> summary(y); hist(y,prob=TRUE,main="Gamma(4.5,1)")
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.2794  2.9440  4.2080  4.5460  5.7500 16.9300 
> curve(dgamma(x,shape=4.5,scale=1),1e-4,20,add=TRUE)
> ## Gamma(0.5,1) (this one is not log-concave)
> y <- arms(runif(1,1e-8,10), function(x) -0.5*log(x)-x,
+           function(x) (x>1e-8)*(x<10), 5000)
> summary(y); hist(y,prob=TRUE,main="Gamma(0.5,1)")
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
9.867e-06 4.265e-02 2.103e-01 4.840e-01 6.471e-01 6.940e+00 
> curve(dgamma(x,shape=0.5,scale=1),1e-8,10,add=TRUE)
> ## Beta(.2,.2) (this one neither)
> y <- arms(runif(1), function(x) (0.2-1)*log(x)+(0.2-1)*log(1-x),
+           function(x) (x>1e-5)*(x<1-1e-5), 5000)
> summary(y); hist(y,prob=TRUE,main="Beta(0.2,0.2)")
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
0.0000215 0.0077390 0.2922000 0.4311000 0.9216000 0.9998000 
> curve(dbeta(x,0.2,0.2),1e-5,1-1e-5,add=TRUE)
> ## Triangular
> y <- arms(runif(1,-1,1), function(x) log(1-abs(x)), function(x) abs(x)<1, 5000)     
> summary(y); hist(y,prob=TRUE,ylim=c(0,1),main="Triangular")
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-0.974600 -0.298400 -0.007794 -0.010020  0.279700  0.985600 
> curve(1-abs(x),-1,1,add=TRUE)
> ## Multimodal examples (Mixture of normals)
> lmixnorm <- function(x,weights,means,sds) {
+     log(crossprod(weights, exp(-0.5*((x-means)/sds)^2 - log(sds))))
+ }
> y <- arms(0, lmixnorm, function(x,...) (x>(-100))*(x<100), 5000, weights=c(1,3,2),
+           means=c(-10,0,10), sds=c(1.5,3,1.5))
> summary(y); hist(y,prob=TRUE,main="Mixture of Normals")
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-14.280  -2.875   1.344   1.675   8.905  14.910 
> curve(colSums(c(1,3,2)/6*dnorm(matrix(x,3,length(x),byrow=TRUE),c(-10,0,10),c(1.5,3,1.5))),
+       par("usr")[1], par("usr")[2], add=TRUE)
> 
> ### Bivariate densities 
> ## Bivariate standard normal
> y <- arms(c(0,2), function(x) -crossprod(x)/2,
+           function(x) (min(x)>-5)*(max(x)<5), 500)
> plot(y, main="Bivariate standard normal", asp=1)
> ## Uniform in the unit square
> y <- arms(c(0.2,.6), function(x) 1,
+           function(x) (min(x)>0)*(max(x)<1), 500)
> plot(y, main="Uniform in the unit square", asp=1)
> polygon(c(0,1,1,0),c(0,0,1,1))
> ## Uniform in the circle of radius r
> y <- arms(c(0.2,0), function(x,...) 1,
+           function(x,r2) sum(x^2)<r2, 500, r2=2^2)
> plot(y, main="Uniform in the circle of radius r; r=2", asp=1)
> curve(-sqrt(4-x^2), -2, 2, add=TRUE)
> curve(sqrt(4-x^2), -2, 2, add=TRUE)
> ## Uniform on the simplex
> simp <- function(x) if ( any(x<0) || (sum(x)>1) ) 0 else 1
> y <- arms(c(0.2,0.2), function(x) 1, simp, 500)
> plot(y, xlim=c(0,1), ylim=c(0,1), main="Uniform in the simplex", asp=1)
> polygon(c(0,1,0), c(0,0,1))
> ## A bimodal distribution (mixture of normals)
> bimodal <- function(x) { log(prod(dnorm(x,mean=3))+prod(dnorm(x,mean=-3))) }
> y <- arms(c(-2,2), bimodal, function(x) all(x>(-10))*all(x<(10)), 500)
> plot(y, main="Mixture of bivariate Normals", asp=1)
> 
> ## A bivariate distribution with non-convex support
> support <- function(x) {
+     return(as.numeric( -1 < x[2] && x[2] < 1 &&
+                       -2 < x[1] &&
+                       ( x[1] < 1 || crossprod(x-c(1,0)) < 1 ) ) )
+ }
> Min.log <- log(.Machine$double.xmin) + 10
> logf <- function(x) {
+     if ( x[1] < 0 ) return(log(1/4))
+     else
+         if (crossprod(x-c(1,0)) < 1 ) return(log(1/pi))
+     return(Min.log)
+ }
> x <- as.matrix(expand.grid(seq(-2.2,2.2,length=40),seq(-1.1,1.1,length=40)))
> y <- sapply(1:nrow(x), function(i) support(x[i,]))
> plot(x,type='n',asp=1)
> points(x[y==1,],pch=1,cex=1,col='green')
> z <- arms(c(0,0), logf, support, 1000)
> points(z,pch=20,cex=0.5,col='blue')
> polygon(c(-2,0,0,-2),c(-1,-1,1,1))
> curve(-sqrt(1-(x-1)^2),0,2,add=TRUE)
> curve(sqrt(1-(x-1)^2),0,2,add=TRUE)
> sum( z[,1] < 0 ) # sampled points in the square
[1] 531
> sum( apply(t(z)-c(1,0),2,crossprod) < 1 ) # sampled points in the circle
[1] 469
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> cleanEx(); ..nameEx <- "convex.bounds"
> 
> ### * convex.bounds
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: convex.bounds
> ### Title: Function to find the boundaries of a convex set
> ### Aliases: convex.bounds
> ### Keywords: misc
> 
> ### ** Examples
> 
> ## boundaries of a unit circle
> convex.bounds(c(0,0), c(1,1), indFunc=function(x) crossprod(x)<1)
[1] -0.7071068  0.7071068
> 
> 
> 
> cleanEx(); ..nameEx <- "trans.dens"
> 
> ### * trans.dens
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: trans.dens
> ### Title: Functions for transdimensional MCMC
> ### Aliases: trans.dens trans.up
> ### Keywords: distribution
> 
> ### ** Examples
> 
> #### ==> Warning: running the examples may take a few minutes! <== ####    
> ### Generate a sample from a mixture of 0,1,2-dim standard normals
> ldens.list <- list(f0 = function(x) sum(dnorm(x,log=TRUE)),
+                    f1 = function(x) dnorm(x,log=TRUE),
+                    f2 = function() 0)
> trans.mix <- function(y) {
+     trans.dens(y, ldens.list=ldens.list, which.models=0:2)
+ }
> 
> trans.rmix <- arms(c(0,0), trans.mix, function(x) crossprod(x)<1e4, 500)
> rmix <- trans.dens(y=trans.rmix, ldens.list=ldens.list,
+                       which.models=0:2, back.transform = TRUE)
> table(rmix[,2])/nrow(rmix) # should be about equally distributed

    0     1     2 
0.338 0.310 0.352 
> plot(trans.rmix,col=rmix[,2]+3,asp=1, xlab="y.1", ylab="y.2",
+      main="A sample from the auxiliary continuous distribution")
> x <- rmix[,-(1:2)]
> plot(x, col=rmix[,2]+3, asp=1,
+      main="The sample transformed back to the original space")
> ### trans.up as a right inverse of trans.dens
> y <- trans.up(x, ldens.list, 0:2)
> stopifnot(all.equal(x, trans.dens(y, ldens.list, 0:2, back.transform=TRUE)[,-(1:2)]))
> 
> ### More trans.up
> z <- trans.up(matrix(0,1000,2), ldens.list, 0:2)
> plot(z,asp=1,col=5) # should look uniform in a circle corresponding to model 2
> z <- trans.up(cbind(runif(1000,-3,3),0), ldens.list, 0:2)
> plot(z,asp=1,col=4) # should look uniform in a region corresponding to model 1
> 
> 
> 
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  18.74 0.06 18.8 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')