# apply.c "LIGHTHOUSE" NESTED SAMPLING APPLICATION # (GNU General Public License software, (C) Sivia and Skilling 2006) # u=0 u=1 # ------------------------------------- # y=2 |:::::::::::::::::::::::::::::::::::::| v=1 # |::::::::::::::::::::::LIGHT::::::::::| # north|::::::::::::::::::::::HOUSE::::::::::| # |:::::::::::::::::::::::::::::::::::::| # |:::::::::::::::::::::::::::::::::::::| # y=0 |:::::::::::::::::::::::::::::::::::::| v=0 # --*--------------*----*--------*-**--**--*-*-------------*-------- # x=-2 coastline -->east x=2 # Problem: # Lighthouse at (x,y) emitted n flashes observed at D[.] on coast. # Inputs: # Prior(u) is uniform (=1) over (0,1), mapped to x = 4*u - 2; and # Prior(v) is uniform (=1) over (0,1), mapped to y = 2*v; so that # Position is 2-dimensional -2 < x < 2, 0 < y < 2 with flat prior # Likelihood is L(x,y) = PRODUCT[k] (y/pi) / ((D[k] - x)^2 + y^2) # Outputs: # Evidence is Z = INTEGRAL L(x,y) Prior(x,y) dxdy # Posterior is P(x,y) = L(x,y) / Z estimating lighthouse position # Information is H = INTEGRAL P(x,y) log(P(x,y)/Prior(x,y)) dxdy ###################################################################### source("mininest.r") ###################################################################### ## object layout ## ------------- ## u Uniform-prior controlling parameter for x ## v Uniform-prior controlling parameter for y ## x Geographical easterly position of lighthouse ## y Geographical northerly position of lighthouse ## logL logLikelihood = ln Prob(data | position) ## logWt log(Weight), adding to SUM(Wt) = Evidence Z # Data set: arrival positions D <- c( 4.73, 0.45, -1.73, 1.09, 2.19, 0.12, 1.31, 1.00, 1.32, 1.07, 0.86, -0.49, -2.59, 1.73, 2.11, 1.61, 4.98, 1.71, 2.23,-57.20, 0.96, 1.25, -1.56, 2.45, 1.19, 2.17,-10.66, 1.91, -4.16, 1.92, 0.10, 1.98, -2.51, 5.55, -0.47, 1.91, 0.95, -0.78, -0.84, 1.72, -0.01, 1.48, 2.70, 1.21, 4.41, -4.79, 1.33, 0.81, 0.20, 1.58, 1.29, 16.19, 2.75, -2.38, -1.79, 6.50, -18.53, 0.72, 0.94, 3.64, 1.94, -0.11, 1.57, 0.57); logLhood <- function( # logLikelihood function x, # Easterly position y) # Northerly position { sum(log((y/pi) / ((D-x)*(D-x) + y*y))) } sample.from.prior <- function() { Obj <- list() Obj$u <- runif(1,0.0,1.0) # uniform in (0,1) Obj$v <- runif(1,0.0,1.0) # uniform in (0,1) Obj$x <- 4.0 * Obj$u - 2.0 # map to x Obj$y <- 2.0 * Obj$v # map to y Obj$logL <- logLhood(Obj$x, Obj$y) Obj } ## Evolve object within likelihood constraint. ## This is different from the C version in that the object passed in ## is not changed, but instead a new version is returned. explore <- function(Obj, # Object being evolved logLstar) # lhood constraint L > Lstar { Ret <- Obj step <- 0.1; # Initial guess suitable step-size in (0,1) m <- 20; accept <- 0; # # MCMC acceptances reject <- 0; # # MCMC rejections Try <- list(); # Trial object for(m in 1:20) { # MCMC counter (pre-judged # steps) ## Create trial object Try$u <- Obj$u + step * (2.*runif(1,0,1) - 1.); # |move| < step Try$v <- Obj$v + step * (2.*runif(1,0,1) - 1.); # |move| < step Try$u <- Try$u - floor(Try$u); # wraparound to stay within (0,1) Try$v <- Try$v - floor(Try$v); # wraparound to stay within (0,1) Try$x <- 4.0 * Try$u - 2.0; # map to x Try$y <- 2.0 * Try$v; # map to y Try$logL <- logLhood(Try$x, Try$y); # trial likelihood value ## Accept if and only if within hard likelihood constraint if( Try$logL > logLstar ) { Ret <- Try; accept <- accept+1; } else reject <- reject+1 ## Refine step-size to let acceptance ratio converge around 50% if( accept > reject ) step <- step * exp(1.0 / accept); if( accept < reject ) step <- step / exp(1.0 / reject); } Ret } ## Output posterior properties, here mean and stddev of x,y show.results <- function(results) { x <- 0.0; xx <- 0.0; # 1st and 2nd moments of x y <- 0.0; yy <- 0.0; # 1st and 2nd moments of y for(s in results$samples) { w <- exp(s$logWt - results$logZ); # proportional weight x <- x + w * s$x; xx <- xx + w * s$x * s$x; y <- y + w * s$y; yy <- yy + w * s$y * s$y; } cat(sprintf("Evidence: ln(Z) = %g +- %g\n",results$logZ, results$logZ.sd)) cat(sprintf("Information: H = %g nats = %g bits\n", results$Hnats, results$Hbits)) cat(sprintf("mean(x) = %g, stddev(x) = %g\n", x, sqrt(xx-x*x))) cat(sprintf("mean(y) = %g, stddev(y) = %g\n", y, sqrt(yy-y*y))) } t0 <- proc.time() results <- nested.sampling(n=100, max.iter=1000, rprior=sample.from.prior, explore=explore) show.results(results) cat('Time elapsed:',(proc.time()-t0)[1],'seconds\n')