Distribution fitting

I’ve just dumped a load of rpanel animations here (all to do with distribution fitting for a Binomial).   (Need to tidy this up when my toothache gets better).

The first bit of code just draws a picture of a standard Binomial mass function.

require(rpanel)
redraw <- function(panel){
rp.tkrreplot(panel, tkrp)
panel
}

draw <- function(panel){
N <- as.numeric(panel$ns)
Np1 <- N+1
t1 <-  bquote( paste(pi== .(panel$prob)))
st1 <- bquote(paste("Number in sample of size ", .(N)))
plot(c(0:N), rep(0.2, Np1), type = "n", ylim = c(0, panel$ylim), main = t1, ylab = "Prob[X=x]", xlab = st1)
y <-  dbinom(c(0:Np1), size=N, prob=panel$prob)
arrows(x0=c(0:N), y0=rep(0, Np1), x1=c(0:N),  y1=y, length = 0.1, col = "red", lwd = 2)

}

ninit <- c(30)

binomial1 <- rp.control(title = "Fitting a distribution", prob=0.5, ns = ninit, ylim = 0.35)

rp.tkrplot(binomial1, tkrp, draw)
rp.slider(binomial1, prob, action = redraw, from = 0.01, to = 0.99, resolution = 0.01)

rp.textentry(binomial1, ns, redraw, c("Size or trial please"),title = "Trial size", initval = ninit) # Input data

rp.slider(binomial1, ylim, action = redraw, from = 0.01, to = 0.99, resolution = 0.01)

x <- rnorm(1000)

The next bit of code does the same thing, but the x axis has been labelled so that instead of giving x, it gives x/n


redraw <- function(panel){
   rp.tkrreplot(panel, tkrp)
   panel
 }

draw <- function(panel){
   t1 <-  bquote( paste(pi== .(panel$prob)))
   plot(c(0:30)/30, rep(0.2, 31), type = "n", ylim = c(0, 0.25), main = t1, ylab = "Prob[X=x]", xlab = "Proportion in sample of size 30")
 y <-  dbinom(c(0:30), size=30, prob=panel$prob)
 arrows(x0=c(0:30)/30, y0=rep(0, 31), x1=c(0:30)/30,  y1=y, length = 0.1, col = "red", lwd = 2)

}

binomial2 <- rp.control(title = "Fitting a distribution", prob=0.5)

rp.tkrplot(binomial2, tkrp, draw)
 rp.slider(binomial2, prob, action = redraw, from = 0.01, to = 0.99, resolution = 0.01)

Finally, the last bit of code takes some data from the Lancet LQAS article, so we can superimpose the pmf on the data using the sliders.


x.lancet <- c(1,3,7,2,6,2,5,3,2,2,4,0,9,8,4,3,6,7,5,5,4,4,3,4)/9

x.lancet.d <- rbinom(30, size=30, prob=0.12)
 x.lancet <- x.lancet.d/30
 xtabs(~x.lancet.d)

## second binomial, now take proportion (30/30)

require(MASS)

redraw <- function(panel){
   rp.tkrreplot(panel, tkrp)
   panel
 }

draw <- function(panel){
   t1 <-  bquote( paste(pi== .(panel$prob)))
   hist(x.lancet, probability = TRUE, col = "pink", main = t1, ylab = "Prob[X=x]", xlab = "Proportion in sample of size 30", xlim = c(0,1), yaxt = "n")
   par(new = TRUE)
    plot(c(0:30)/30, rep(0.2, 31), type = "n", ylim = c(0, 0.25), main = "", ylab = "", xlab = "")
 y <-  dbinom(c(0:30), size=30, prob=panel$prob)
 arrows(x0=c(0:30)/30, y0=rep(0, 31), x1=c(0:30)/30,  y1=y, length = 0.1, col = "red", lwd = 2)

}

binomial3 <- rp.control(title = "Fitting a distribution", prob=0.5)

rp.tkrplot(binomial3, tkrp, draw)
 rp.slider(binomial3, prob, action = redraw, from = 0.01, to = 0.99, resolution = 0.01)

And finally, if all the above has made sense, we can plot the likelihood on a panel to the right of the data fitting.


lf <- function(prob){

llike <- sum(dbinom(x.lancet.d, size = 30, prob=prob, log = TRUE))
   return(llike)
 }

grid <- seq(from=0.01, to = 0.99, by = 0.01)

llike <- vector("numeric", length(grid))
 for (i in 1:length(grid)){
   llike[i] <- lf(grid[i])
 }

#min(llike)
 #max(llike)

#plot(grid, llike, type = "l")

#require(MASS)

redraw <- function(panel){
   rp.tkrreplot(panel, tkrp)
   panel
 }

draw <- function(panel){
   par(mfrow = c(1,2))
   t1 <-  bquote( paste(pi== .(panel$prob)))
   hist(x.lancet, probability = TRUE, col = "pink", main = t1, ylab = "Prob[X=x]", xlab = "Proportion in sample of size 30", xlim = c(0,1), yaxt = "n")
   par(new = TRUE)
    plot(c(0:30)/30, rep(0.2, 31), type = "n", ylim = c(0, 0.25), main = "", ylab = "", xlab = "")
 y <-  dbinom(c(0:30), size=30, prob=panel$prob)
 arrows(x0=c(0:30)/30, y0=rep(0, 31), x1=c(0:30)/30,  y1=y, length = 0.1, col = "red", lwd = 2)
   plot(grid, llike, type = "n")
   points(panel$prob, lf(panel$prob), pch = 16, col = "blue")

}

binomial4<- rp.control(title = "Fitting a distribution", prob=0.5)

rp.tkrplot(binomial4, tkrp, draw)
 rp.slider(binomial4, prob, action = redraw, from = 0.01, to = 0.99, resolution = 0.01)

At some point, for some reason, I have pulled out a bit of code that plots the likelihood function alone.


draw <- function(panel){

#n <- 30##n=300
   #s <- 5##s=112
   ##R <- N-S
   par(mfrow = c(1,2))
   t1 <-  bquote( paste(alpha== .(panel$A), ' ', beta==.(panel$B)))
   t2 <-  bquote( paste(n== .(plotit[1]), ' ', s==.(plotit[2])))
 curve(dbeta(x, shape1=panel$A, shape2=panel$B, log = TRUE), from = 0.01, to = 0.99, n=99, ylab = "", xlab = "Prior", lwd = 2, col = "blue", main = t1)
 curve(dbeta(x, shape1=(plotit[2]+panel$A), shape2=((plotit[1]-plotit[2])+panel$B), log = TRUE), from = 0.01, to = 0.99, n=99, ylim = c(-100, 0), xlab = "Posterior", ylab = "", lwd = 2, col = "red", main = t2)

}

redraw <- function(panel){
   rp.tkrreplot(panel, tkrp)
   panel
 }

# n <- 30##n=300
  # s <- 5##s=112

N <- 6
 S <- 1

N <- 900
 S <- 112
 plotit <- c(N,S)

binomial5 <- rp.control(title = "Fitting a distribution", A=1, B=1)

rp.tkrplot(binomial5, tkrp, draw)
 rp.slider(binomial5, A, action = redraw, from = 0.01, to = 10, resolution = 0.01)
 rp.slider(binomial5, B, action = redraw, from = 0.01, to = 10, resolution = 0.01)

I guess the logic was next to add the ability to have a prior distribution for prob[X=x].

The main todo here is to stitch all of these into one user friendly panel.

I do wonder what happened to my rejection sampler code.

Advertisements
This entry was posted in R, rpanel and tagged , , , . Bookmark the permalink.

2 Responses to Distribution fitting

  1. subramanian says:

    Hi

    Thank you for the examples. i have a query. How do i make rpanel print the value of a calculation in the rpanel itself – say, just below a button? Assume that the panel has two rp.textentry boxes for variables x and y, and an rp.button and i want the sum of the 2 variables and y to be printed below the submit button as say Answer: value

  2. subramanian says:

    it should be “..2 variables x and y to be printed.. ” in my post above

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