The density plot I gave when describing different sources of error in the variance-bias decomposition post was made up, with two arbitrary normal densities plotted together. So, it would be nice if I could give a real example.

For the sake of simplicity, say we have one parameter with a simple-normal distribution, and we have one observation that’s simple-normally distributed around the parameter. The true posterior is then also a normal distribution. However, since the ABC posterior is conditional on the observation in a certain ball, we have an expression in terms of the difference between two distribution functions. The resulting curve is almost normal, but not quite.

This isn’t too hard to plot for some specific values, but really what we’d like is a way to change values and see how the plot changes. Ideally I’d like to find a way to get this into a PDF, so I can use it during presentations, but at the moment I’ve got it running in an R script.

Code is below. It’s nothing too fancy, and was mostly a matter of fooling around with the *tcltk* package demo code to see what happened. Eventually I’d like to add:

-Sliders for prior and data variance

-Shifting the plot axes as the variables change

-A way to embed this in a PDF, if that’s possible

-Exception handling for zero tolerance

cdfs <- function(x.star,delta,cmean,cvar) { #cdf of x*+delta-cmean/sqrt(cvar) minus that of x*-delta-cmean/sqrt(cvar) pnorm(x.star, mean=cmean-delta, sd=sqrt(cvar) ) -pnorm(x.star, mean=cmean+delta, sd=sqrt(cvar) ) } dabcpost <- function(x,x.star,pmean,pvar,datvar,delta) { dnorm(x,mean=pmean,sd=sqrt(pvar) ) * cdfs(x.star,delta,x,datvar) /cdfs(x.star,delta,pmean,pvar+datvar) } ABCplot <- function(x.star,delta,xlim,ylim) { curve(dnorm(x,mean=x.star/2,sd=1/sqrt(2) ) , lty="dashed",xlab=expression(theta) , ylab="Density",main="ABC Posterior Error", xlim=xlim,ylim=ylim) curve(dnorm,lty="dotted",add=T) curve(dabcpost(x,x.star,pmean=0,pvar=1,datvar=1, delta=delta) , add=T) legend(x="topright",c("Prior","True posterior","ABC posterior") , lty=c("dotted","dashed","solid") ) } require(tcltk) || stop("tcltk support is absent") require(graphics); require(stats) local({ have_ttk <- as.character(tcl("info", "tclversion")) >= "8.5" if(have_ttk) { tkbutton <- ttkbutton tkframe <- ttkframe tklabel <- ttklabel tkradiobutton <- ttkradiobutton } xlim <- c(-5,5) ylim <- c(0,0.6) x.star <- tclVar(3) x.star.sav <- 3 bw <- tclVar(1) bw.sav <- 1 # in case replot.maybe is called too early replot <- function(...) { bw.sav <<- b <- as.numeric(tclObj(bw)) x.star.sav <<- xs <- as.numeric(tclObj(x.star)) eval(substitute(ABCplot(xs,b,xlim,ylim))) } replot.maybe <- function(...) { if (as.numeric(tclObj(bw)) != bw.sav || as.numeric(tclObj(x.star)) != x.star.sav) replot() } regen <- function(...) { xlim <<- c(min(0,as.numeric(tclObj(x.star) ) /2) -5, max(0,as.numeric(tclObj(x.star) ) /2) +5) replot() } grDevices::devAskNewPage(FALSE) # override setting in demo() tclServiceMode(FALSE) base <- tktoplevel() tkwm.title(base, "Density") spec.frm <- tkframe(base,borderwidth=2) right.frm <- tkframe(spec.frm) frame3 <-tkframe(right.frm,relief="groove",borderwidth=2) tkpack(tklabel(frame3,text="Observation") ) tkpack(tkscale(frame3,command=replot.maybe,from=1,to=10, showvalue=T,variable=x.star, resolution=0.1,orient="horiz") ) frame4 <-tkframe(right.frm, relief="groove", borderwidth=2) tkpack(tklabel (frame4, text="Tolerance")) tkpack(tkscale(frame4, command=replot.maybe, from=0.05, to=16.00, showvalue=T, variable=bw, resolution=0.05, orient="horiz")) tkpack(frame3,frame4, fill="x") tkpack(right.frm,side="left", anchor="n") ## `Bottom frame' (on base): q.but <- tkbutton(base,text="Quit", command=function() tkdestroy(base)) tkpack(spec.frm, q.but) tclServiceMode(TRUE) regen() })