Bayesian_chocolates.Rmd

---
title: "Intro to Bayesian Inference"
author: "Prof. Gwen Eadie"
date: "08/03/2021"
output: html_document
---



#  Beta prior distribution


```{r}
# parameters for the prior distribution
a = 3
b = 8

# set the margins to make room for labels (don't have to do this, there are defaults)
# numbers here are bottom, left, top, right
par(mar=c(5,7,4,2))

# with curve(), with fancy labels
curve( dbeta(x, shape1 = a, shape2 = b), xlab = expression(theta), ylab = expression(p(theta)) , cex.lab=2, cex.axis=2, lwd = 2.5, main="Prior Distribution for proportion of blue m&m's")
# ... and a grid:
grid()

```



# Write a function for plotting the posterior distribution

```{r}
# function to plot posterior distribution
posteriorfunc = function(x, a, b, n, y){
  
  dbeta(x, shape1 = (a + y), shape2 = (b + n - y))
  
}

```

# Data

```{r}

# number of blue m&m's 
nblue = 28
# total number of m&m's
ntotal = 120
  
# plot posterior  
par(mar=c(5,5,2,2))
curve( posteriorfunc(x, a=a, b=b, n=ntotal, y=nblue),  xlab = expression(theta), ylab = expression(p(theta|y)) , cex.lab=2, cex.axis=2, lwd = 2.5, main="Posterior Distribution for proportion of blue m&m's")

# add a curve showing the prior distribution
curve( dbeta(x, shape1 = a, shape2 = b), add=TRUE, col="blue", lty=2)

# add a legend
legend("topright", legend=c("prior", "posterior"), col = c("blue", "black"), lwd = 2.5, lty = c(2,1))

nblue = 2
  
# total number of m&m's
ntotal = 20
```

  
# Plot the posterior distribution, calculate some summary statistics

```{r}
```

# Calculate some summary statistics

```{r}
# e.g. function to calculate the expected value of theta (the mean) of the posterior
mean.beta = function(a, b) a/(a+b)

mean.beta(a=(a+nblue), b=(b+ntotal-nblue))  


```