## 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), ' ', s==.(plotit)))
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+panel\$A), shape2=((plotit-plotit)+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