Showing Bayesian updating in R I'm very new to Bayesian and am very interested in understanding the concept that each posterior from a previous test, can be used as a prior for a current test or in Lindley's (2000) words 'Yesterday's posterior is today's prior'.
To see this in action, however, I was wondering if I could explore this concept in R as a function for the simple problem of estimating the probability of success from binomially distributed data.
I tried something but my R skills were not advanced enough to show this updating concept. Since Stack Overflow was not of much help, I am wondering if any of our colleagues here might know the trick? 
(In R code below, pr = prior; lk = likelihood; ps = posterior)
update <- function(n, dat){

for(i in 1:length(n)){           
pr = function(x) dbeta(x, 1, 1)
lk = function(x) dbinom(dat[i], n[i], x)  # So, here first `n = 100` and 
ps = function(x) pr(x)*lk(x)              # `dat = 55` will go thru first 
 }                                        # round of the loop and produce
curve(ps)                                 # a `ps`. But in the second run of
}                                         # the loop, the `ps` just
# Example of use:                         # produced will be used as `pr`
update(n = c(100, 50), dat = c(55, 60) )

 A: You aren't updating the "new" prior with the posterior!  And it's much harder to do with functions, unless you're good with R functional programming, than with discretized values for $p$.
I'm constructing an example similar to what you were doing, but not exactly the same.  Here I assume the data follows a Binomial($n=10$, $p=0.5$) distribution,that $n=10$ is known, and that the initial prior on $p$ is uniform.  At each of 20 iterations I generate one new data point from the Binomial(10, 0.5) distribution and update the posterior with it.
Here's some R code that will do the job:
# True distribution
p <- 0.5
n <- 10

# Prior on p (n assumed known), discretized
p_values <- seq(0.01,0.99,0.001)
pr <- dbeta(p_values,1,1)
pr <- 1000 * pr / sum(pr)  # Have to normalize given discreteness
plot(pr~p_values, col=1, ylim=c(0, 14), type="l")

# Run for 20 samples
for (i in 1:20) {
   x <- rbinom(1, n, p) 
   ps <- dbinom(x, n, p_values) * pr
   ps <- 1000 * ps / sum(ps)
   lines(ps~p_values, col=(i+1))

   pr = ps
}

And here's the graphical result:

As you can see, the posterior distributions do concentrate rather quickly and, despite a poor initial sample, are shifting their peak closer and closer to the true value of 0.5, albeit with some randomness.
