R Script: Interactive Plot for ABC Posterior Against True Posterior

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.

Interactive1

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()
})
Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s