x <- seq(0,10)
plot(x, dbinom(x,10, prob=.5))
abline(v=7, lty=2)

pi <- c(.25,.50,.75)
prior <- rep(1/3,3)

plot(pi, prior, main="Prior Belief", pch=16)


dbinom(7,10,.25)

dbinom(7,10,.25)*(1/3)

post <- dbinom(7,10,pi)*prior
plot(pi,post, pch=16, main="Posterior, Un-normalized")

post.norm <- post/sum(post)
plot(pi, post.norm, pch=16, main="Posterior, Normalized")
post.norm

#Set priors
prior <- c(.025,.975,.025)
plot(pi, prior, pch=16, main="Prior")

#Calculate likelihoods*priors
post <- dbinom(7,10,pi)*prior

#Normalise by dividing by the sum
post.norm <- post/sum(post)

#Display posterior
plot(pi, prior, pch=16, main="Posterior", col="gray80")
points(pi, post.norm, pch=16)
post.norm


#Calculate likelihoods*priors (This is where the change in the data is)
post <- dbinom(70,100,pi)*prior

#Normalise by dividing by the sum
post.norm <- post/sum(post)

#Display posterior
plot(pi, prior, pch=16, main="Posterior", col="gray80")
points(pi, post.norm, pch=16)
post.norm

pi <- seq(0,1,.001)

prior <- rep(1/1001,1001)
plot(pi, prior, type="l", main="Prior")

post <- dbinom(7,10,pi)*prior
post.norm <- post/sum(post)

plot(pi, post.norm, type="l", main="Posterior")


sum(post.norm[pi>.5])

prior <- dbeta(pi, 100,100)
prior <- prior/sum(prior)
plot(pi, prior, type="l", main="Strong Prior")


post <- dbinom(7,10,pi)*prior
post.norm <- post/sum(post)
plot(pi, prior, type="l", col="gray80", main="Posterior")
points(pi, post.norm, type="l")

post <- dbinom(70,100,pi)*prior
post.norm <- post/sum(post)
plot(pi, prior, type="l", col="gray80", main="Posterior")
points(pi, post.norm, type="l")

sum(post.norm[pi<.6]) - sum(post.norm[pi<.55])

#Calculate cumulative probability
cumu.post <- cumsum(post.norm)
#Find the boundaries
pi[min(which(cumu.post>=.975))]
pi[max(which(cumu.post<=.025))]


x <- rbinom(1000,1,.52)
mean(x)

t.test(x, mu=.50)

dem.support <- seq(0,1,.01)
prior <- dbeta(dem.support, 50,50)
prior <- prior/sum(prior)

plot(dem.support, prior, type="l", main="Prior")

post <- dbinom(514,1000,dem.support)*prior
post.norm <- post/sum(post)

plot(dem.support, prior, type="l", main="Prior", col="gray80")
points(dem.support, post.norm, type="l")
abline(v=.5, lty=2)

sum(post.norm[dem.support>=.55])

#Calculate cumulative probability
cumu.post <- cumsum(post.norm)
#Find the boundaries
dem.support[min(which(cumu.post>=.975))]
dem.support[max(which(cumu.post<=.025))]

prior <- dbeta(dem.support, 510+550+495, 490+450+505)
prior <- prior/sum(prior)
plot(dem.support, prior, type="l", main="Prior")
abline(v=.5,lty=2)

post <- dbinom(514,1000,dem.support)*prior
post.norm <- post/sum(post)
plot(dem.support, prior, type="l", main="Posterior", col="gray80")
points(dem.support, post.norm, type="l")
abline(v=.5,lty=2)

#Calculate cumulative probability
cumu.post <- cumsum(post.norm)
#Find the boundaries
dem.support[min(which(cumu.post>=.975))]
dem.support[max(which(cumu.post<=.025))]

post.poll <- function(p) {
  if (p <= 0 || p >= 1) return(0)  #p has to be a probability
  dbinom(514, 1000, p) * dbeta(p, 50, 50)
}

set.seed(19104)
n.iter <- 20000
p.draws <- rep(NA, n.iter)

p.cur    <- 0.5                    #start somewhere
post.cur <- post.poll(p.cur)
step     <- 0.02                   #how big a step we propose

for (i in 1:n.iter) {
  p.prop    <- p.cur + rnorm(1, 0, step)   #propose a nearby value
  post.prop <- post.poll(p.prop)
  #move there with probability equal to the ratio of the two heights
  if (runif(1) < post.prop / post.cur) {
    p.cur    <- p.prop
    post.cur <- post.prop
  }
  p.draws[i] <- p.cur                       #write down where we are
}

#throw away the early "burn-in" steps that depend on where we started
p.draws <- p.draws[2001:n.iter]

#Fine-grid posterior as the benchmark
p.grid <- seq(0, 1, 0.001)
prior  <- dbeta(p.grid, 50, 50); prior <- prior / sum(prior)
grid.post <- dbinom(514, 1000, p.grid) * prior
grid.post <- grid.post / sum(grid.post)

hist(p.draws, breaks = 40, freq = FALSE, xlim = c(.45, .58),
     main = "MCMC draws vs. the grid posterior",
     xlab = "p (Democratic support)")
#overlay the grid posterior (rescaled to a density for comparison)
points(p.grid, grid.post / mean(diff(p.grid)), type = "l", col = "firebrick", lwd = 2)

mean(p.draws)                     #posterior mean
quantile(p.draws, c(.025, .975))  #95% credible interval
mean(p.draws > 0.5)               #P(Democrat ahead)

set.seed(19104)
n <- 50
x <- runif(n, 0, 10)
y <- 2 + 0.5 * x + rnorm(n, mean = 0, sd = 1)
plot(x, y, pch = 16, main = "Our data")

#Log-score of a candidate line: bigger (less negative) = better fit
log.post <- function(alpha, beta) {
  sum(dnorm(y, mean = alpha + beta * x, sd = 1, log = TRUE)) +
  dnorm(alpha, 0, 10, log = TRUE) + dnorm(beta, 0, 10, log = TRUE)
}

plot(x, y, pch = 16, main = "Three candidate lines")
abline(a = 0, b = 1,   col = "firebrick")   #(0, 1)
abline(a = 4, b = 0,   col = "dodgerblue")  #(4, 0)
abline(a = 2, b = 0.5, col = "forestgreen") #(2, 0.5)

log.post(0, 1)    #red
log.post(4, 0)    #blue
log.post(2, 0.5)  #green

set.seed(19104)
n.iter <- 20000
draws <- matrix(NA, nrow = n.iter, ncol = 2)
colnames(draws) <- c("alpha", "beta")

alpha.cur <- 0; beta.cur <- 0
lp.cur <- log.post(alpha.cur, beta.cur)
step <- 0.15

for (i in 1:n.iter) {
  alpha.prop <- alpha.cur + rnorm(1, 0, step)
  beta.prop  <- beta.cur  + rnorm(1, 0, step)
  lp.prop    <- log.post(alpha.prop, beta.prop)
  #move with probability equal to the ratio of the two heights
  if (runif(1) < exp(lp.prop - lp.cur)) {
    alpha.cur <- alpha.prop; beta.cur <- beta.prop; lp.cur <- lp.prop
  }
  draws[i, ] <- c(alpha.cur, beta.cur)
}
draws <- draws[2001:n.iter, ]

plot(draws[, "alpha"], draws[, "beta"], pch = 16, cex = .3,
     col = rgb(0, 0, 0, .1),
     xlab = "Intercept (alpha)", ylab = "Slope (beta)",
     main = "Where the walker went")

plot(x, y, pch = 16, main = "The posterior is a cloud of lines")
for (i in seq(1, nrow(draws), length.out = 200)) {
  abline(a = draws[i, "alpha"], b = draws[i, "beta"], col = rgb(.2, .2, .8, .05))
}
points(x, y, pch = 16)

colMeans(draws)                            #posterior means
quantile(draws[, "beta"], c(.025, .975))   #95% credible interval for the slope
mean(draws[, "beta"] > 0)                  #P(slope is positive)

confint(lm(y ~ x))

set.seed(19104)

#Two-tailed p-value for 7 heads in 10 flips, fair coin
p.value <- pbinom(2, 10, 0.5) + (1 - pbinom(6, 10, 0.5))
p.value

#Define a grid of possible pi values
pi.vals <- seq(0, 1, 0.001)

#Uniform prior: every value equally likely
prior <- rep(1, length(pi.vals))
#Normalize so it sums to 1
prior <- prior / sum(prior)

plot(pi.vals, prior, type="l",
     xlab="Pi (Probability of Heads)", ylab="Prior Belief",
     main="Uniform Prior", ylim=c(0, max(prior)*3))

#Likelihood: P(7 heads in 10 flips | pi) for each possible pi
likelihood <- dbinom(7, 10, pi.vals)

plot(pi.vals, likelihood, type="l",
     xlab="Pi (Probability of Heads)", ylab="Likelihood",
     main="Likelihood of 7 Heads in 10 Flips")

#Posterior = Prior x Likelihood (then normalize)
posterior <- prior * likelihood
posterior <- posterior / sum(posterior)

plot(pi.vals, posterior, type="l",
     xlab="Pi (Probability of Heads)", ylab="Posterior Belief",
     main="Posterior Distribution After 7/10 Heads")
abline(v=0.5, lty=2, col="gray50")
abline(v=0.7, lty=2, col="firebrick")
legend("topleft", c("Pi = 0.5", "Pi = 0.7"),
       lty=2, col=c("gray50","firebrick"))

#What is the probability that the coin is biased towards heads?
prob.biased <- sum(posterior[pi.vals > 0.5])
prob.biased

#Strong prior: we think the coin is probably fair
strong.prior <- dbeta(pi.vals, 20, 20)
strong.prior <- strong.prior / sum(strong.prior)

#Weak prior (uniform)
weak.prior <- rep(1, length(pi.vals))
weak.prior <- weak.prior / sum(weak.prior)

plot(pi.vals, strong.prior, type="l", col="dodgerblue", lwd=2,
     xlab="Pi", ylab="Prior Belief", main="Two Different Priors")
points(pi.vals, weak.prior, type="l", col="firebrick", lwd=2)
legend("topright", c("Strong prior (probably fair)", "Weak prior (uniform)"),
       col=c("dodgerblue","firebrick"), lwd=2)

#Posterior with strong prior
posterior.strong <- strong.prior * likelihood
posterior.strong <- posterior.strong / sum(posterior.strong)

#Posterior with weak prior
posterior.weak <- weak.prior * likelihood
posterior.weak <- posterior.weak / sum(posterior.weak)

plot(pi.vals, posterior.weak, type="l", col="firebrick", lwd=2,
     xlab="Pi", ylab="Posterior Belief",
     main="Posteriors Under Different Priors",
     ylim=c(0, max(c(posterior.weak, posterior.strong))))
points(pi.vals, posterior.strong, type="l", col="dodgerblue", lwd=2)
abline(v=0.7, lty=3)
abline(v=0.5, lty=3)
legend("topleft", c("Weak prior (uniform)", "Strong prior (probably fair)"),
       col=c("firebrick","dodgerblue"), lwd=2)

#Likelihood with more data
likelihood.big <- dbinom(700, 1000, pi.vals)

#Posterior with strong prior and lots of data
posterior.strong.big <- strong.prior * likelihood.big
posterior.strong.big <- posterior.strong.big / sum(posterior.strong.big)

#Posterior with weak prior and lots of data
posterior.weak.big <- weak.prior * likelihood.big
posterior.weak.big <- posterior.weak.big / sum(posterior.weak.big)

plot(pi.vals, posterior.weak.big, type="l", col="firebrick", lwd=2,
     xlab="Pi", ylab="Posterior Belief",
     main="Posteriors Converge With More Data (n=1000)",
     ylim=c(0, max(c(posterior.weak.big, posterior.strong.big))))
points(pi.vals, posterior.strong.big, type="l", col="dodgerblue", lwd=2)
legend("topleft", c("Weak prior", "Strong prior"),
       col=c("firebrick","dodgerblue"), lwd=2)

set.seed(19104)

#True probability of heads
true.pi <- 0.6

#Generate 100 flips
flips <- rbinom(100, 1, true.pi)

#Start with a uniform prior (Beta(1,1))
a <- 1
b <- 1

#Track the posterior mean after each flip
posterior.means <- rep(NA, 100)

par(mfrow=c(2,3))
for(i in 1:100){
  if(flips[i] == 1){
    a <- a + 1  #Heads: increase a
  } else {
    b <- b + 1  #Tails: increase b
  }

  posterior.means[i] <- a / (a + b)

  #Plot the posterior at selected points
  if(i %in% c(1, 5, 10, 25, 50, 100)){
    curve(dbeta(x, a, b), from=0, to=1,
          xlab="Pi", ylab="Density",
          main=paste("After", i, "flips"))
    abline(v=0.6, lty=2, col="firebrick")
  }
}
par(mfrow=c(1,1))

plot(1:100, posterior.means, type="l",
     xlab="Number of Flips", ylab="Posterior Mean (Best Guess for Pi)",
     main="Bayesian Updating of Coin Bias Estimate")
abline(h=0.6, lty=2, col="firebrick")
legend("topright", "True Pi = 0.6", lty=2, col="firebrick")

set.seed(19104)

pi.vals <- seq(0, 1, 0.001)

#Prior: centered on 0.50, most mass between .45 and .55
#Beta(50, 50) has mean 0.5 and is concentrated
prior <- dbeta(pi.vals, 50, 50)
prior <- prior / sum(prior)

plot(pi.vals, prior, type="l", lwd=2,
     xlab="Democratic Two-Party Vote Share",
     ylab="Prior Belief",
     main="Prior Belief About Election Outcome",
     xlim=c(.3,.7))

#Frequentist approach
n.poll <- 1000
p.hat <- 0.53
se <- sqrt(p.hat * (1 - p.hat) / n.poll)

#95% confidence interval
ci.lower <- p.hat - 1.96 * se
ci.upper <- p.hat + 1.96 * se
cat("Frequentist 95% CI:", round(ci.lower, 3), "to", round(ci.upper, 3), "\n")

#Hypothesis test: H0: pi = 0.5
z <- (p.hat - 0.5) / se
p.val <- 2 * pnorm(-abs(z))
cat("P-value for H0: pi = 0.5:", round(p.val, 4))

#Likelihood: 530 out of 1000 support Democrat
likelihood <- dbinom(530, 1000, pi.vals)

#Posterior
posterior <- prior * likelihood
posterior <- posterior / sum(posterior)

plot(pi.vals, prior, type="l", lwd=2, col="gray60",
     xlab="Democratic Two-Party Vote Share",
     ylab="Density",
     main="Prior, Likelihood, and Posterior",
     xlim=c(.4, .6),
     ylim=c(0, max(c(prior, posterior))*1.1))
points(pi.vals, likelihood / max(likelihood) * max(posterior),
       type="l", lwd=2, col="darkorange", lty=2)
points(pi.vals, posterior, type="l", lwd=2, col="dodgerblue")
legend("topleft", c("Prior", "Likelihood (scaled)", "Posterior"),
       col=c("gray60","darkorange","dodgerblue"), lwd=2, lty=c(1,2,1))

#Probability that the Democrat is ahead
prob.dem.ahead <- sum(posterior[pi.vals > 0.5])
cat("Probability Democrat is ahead:", round(prob.dem.ahead, 3), "\n")

#Bayesian "credible interval" -- central 95% of the posterior
#Find the 2.5th and 97.5th percentiles of the posterior
cumulative <- cumsum(posterior)
ci.bayes.lower <- pi.vals[min(which(cumulative >= 0.025))]
ci.bayes.upper <- pi.vals[min(which(cumulative >= 0.975))]
cat("Bayesian 95% Credible Interval:", ci.bayes.lower, "to", ci.bayes.upper)

set.seed(19104)

#God says: true support is 52%
true.pi <- 0.52
n <- 500

#Take a sample
sample.data <- rbinom(n, 1, true.pi)
p.hat <- mean(sample.data)

cat("Sample proportion:", round(p.hat, 3), "\n")
cat("Sample size:", n, "\n")

#----- Frequentist -----
se.freq <- sqrt(p.hat * (1 - p.hat) / n)

#95% CI
freq.lower <- p.hat - 1.96 * se.freq
freq.upper <- p.hat + 1.96 * se.freq

#Test H0: pi = 0.5
z.stat <- (p.hat - 0.5) / se.freq
p.val <- 2 * pnorm(-abs(z.stat))

cat("--- Frequentist Results ---\n")
cat("95% CI:", round(freq.lower, 3), "to", round(freq.upper, 3), "\n")
cat("P-value (H0: pi = 0.5):", round(p.val, 4), "\n")
cat("Reject null at alpha = 0.05?", ifelse(p.val < 0.05, "Yes", "No"), "\n")

#----- Bayesian -----
#Prior: Beta(50, 50) -- elections are usually close
a.prior <- 50
b.prior <- 50

#Posterior: Beta(a + successes, b + failures)
successes <- sum(sample.data)
failures <- n - successes
a.post <- a.prior + successes
b.post <- b.prior + failures

#Posterior mean
post.mean <- a.post / (a.post + b.post)

#95% credible interval
cred.lower <- qbeta(0.025, a.post, b.post)
cred.upper <- qbeta(0.975, a.post, b.post)

#P(Democrat ahead)
prob.ahead <- 1 - pbeta(0.5, a.post, b.post)

cat("--- Bayesian Results ---\n")
cat("Posterior mean:", round(post.mean, 3), "\n")
cat("95% Credible Interval:", round(cred.lower, 3), "to", round(cred.upper, 3), "\n")
cat("P(Democrat ahead):", round(prob.ahead, 3), "\n")

#Visualize
curve(dbeta(x, a.post, b.post), from=0.4, to=0.65,
      xlab="Democratic Vote Share", ylab="Density",
      main="Bayesian Posterior with Frequentist CI",
      lwd=2, col="dodgerblue")
abline(v=true.pi, lty=2, col="firebrick", lwd=2)
abline(v=cred.lower, lty=3, col="dodgerblue")
abline(v=cred.upper, lty=3, col="dodgerblue")
abline(v=freq.lower, lty=3, col="darkorange")
abline(v=freq.upper, lty=3, col="darkorange")
legend("topright",
       c("Posterior", "True value", "Bayesian 95% CI", "Frequentist 95% CI"),
       col=c("dodgerblue","firebrick","dodgerblue","darkorange"),
       lwd=c(2,2,1,1), lty=c(1,2,3,3))

n.small <- 50
y.small <- 28
p.hat.small <- y.small / n.small

#Frequentist CI
se.small <- sqrt(p.hat.small * (1 - p.hat.small) / n.small)
cat("Frequentist 95% CI:",
    round(p.hat.small - 1.96*se.small, 3), "to",
    round(p.hat.small + 1.96*se.small, 3), "\n")

#Prior: Beta(24, 26) has mean ~0.48
a.prior.small <- 24
b.prior.small <- 26

a.post.small <- a.prior.small + y.small
b.post.small <- b.prior.small + (n.small - y.small)

post.mean.small <- a.post.small / (a.post.small + b.post.small)
cred.lower.small <- qbeta(0.025, a.post.small, b.post.small)
cred.upper.small <- qbeta(0.975, a.post.small, b.post.small)
prob.ahead.small <- 1 - pbeta(0.5, a.post.small, b.post.small)

cat("Bayesian posterior mean:", round(post.mean.small, 3), "\n")
cat("Bayesian 95% Credible Interval:",
    round(cred.lower.small, 3), "to", round(cred.upper.small, 3), "\n")
cat("P(Your candidate ahead):", round(prob.ahead.small, 3), "\n")