Two R functions for Bayesian updating with different output Basically, I have made two posterior (ps) updating R functions; auto1 and auto2. That is, ps from one set of n and s becomes prior (pr) for the new set of n and s.
But it seems the two functions produce different posteriors, I'm wondering which function is accurate; auto1 or auto2?
auto1 <- function(n, s, top = 5){

  props <- seq(0, 1, 1e-4)
    pr <- dbeta(props, 1, 1)
    pr <- 1e4 * pr / sum(pr)
plot(pr~props, ylim = c(0, top*length(n)), type = "l", yaxt = "n", ylab = NA)

for(i in 1:length(n)) {
ps <- dbinom(s[i], n[i], props) * pr
ps <- 1e4 * ps / sum(ps)
lines((ps+i+1)~props, col = i +1)
pr <- ps
  }
}
auto1(n = c(100, 90, 120), s = c(55, 60, 80), top = 7)

#=============================================================
auto2 <- function(n, s, top = 2.5){

  U <- function(x) dbeta(x, 1, 1)
  props <- seq(0, 1, 1e-4)
  pr <- U(props)

plot(pr~props, ylim = c(0, top*length(n)), type = "l", lty = 2, yaxt = "n", ylab = NA)

for(i in 1:length(n)) {

lk <- function(x) dbinom(s[i], n[i], x)
k <- integrate(function(x) U(x)*lk(x), 0, 1)[[1]]
ps <- function(x) U(x)*lk(x)/k

lines((ps(props)+i+1)~props, col = i +1)
pr <- ps
   }
}
auto2(n = c(100, 90, 120), s = c(55, 60, 80), top = 5)

 A: Note that in auto2 you are never actually using pr inside the for-loop, you are always referring back to U(x), which is the uniform prior.  So... no updating is occurring across observations, only the uniform prior is being updated with the single observation indexed by i.
Examine this modification with care:
auto2 <- function(n, s, top = 2.5){

   prior <- function(new_x, x, y) approx(x, y=y, xout=new_x)$y
   props <- seq(0, 1, 1e-4)
   pr <- dbeta(props,1,1)

   plot(pr~props, ylim = c(0, top*length(n)), type = "l", lty = 2, yaxt = "n", ylab = NA)

   for(i in 1:length(n)) {

      lk <- function(x) dbinom(s[i], n[i], x)
      k <- integrate(function(x) prior(x,props,pr)*lk(x), 0, 1)[[1]]
      ps <- function(x) prior(x,props,pr)*lk(x)/k

      lines((ps(props)+i+1)~props, col = i +1)
      pr <- ps(props)
   }
}
auto2(n = c(100, 90, 120), s = c(55, 60, 80), top = 7)

approx does linear interpolation between the values passed to it.  It provides a functional version of the prior, instead of just the list of values provided in auto1.  The functional version is based on the x-values, in this case props, and the values of the distribution function, stored in pr, but, because it's a function, can be integrated.  
At every iteration / value of i we update the vector pr, which acts to update the values returned from prior with the just-formed posterior.  You will note that the two sets of curves are identical now!
Edit in response to a follow-up question:
You can also get confidence intervals etc. from the two functions.  Consider auto1; edit it so it returns list(pr=pr, props=props).  We can then find a confidence interval using the following function:
res <- auto1(n = c(100, 90, 120), s = c(55, 60, 80), top = 7)

ci_1 <- function(res, alpha) {
       res$cdf <- cumsum(res$pr)
       list(lower=res$props[max(which(res$cdf<alpha))],
            upper=res$props[min(which(res$cdf>(1-alpha)))])   
}

which in this case results in:
> ci_1(res, alpha=0.05)
$lower
[1] 0.504

$upper
[1] 0.5236

Note that this is not the purist's confidence interval, although it is one; the purist's CI would be the shortest of all confidence intervals that contain, in this case, 95% of the posterior mass.  That's not hard to find, but in the case of nearly symmetric posteriors, the two will be very similar.
