A nice aspect of the book is that he uses many examples in R.

**Grid approximation for estimating a posterior distribution**

In chapter 2, he shows how to estimate a posterior distribution using a 'grid approximation' approach. I've wrapped up his lines of code into a little function that lets you try out different number of grid points.

In this case, we are interested in calculating the posterior distribution for

*p*, the probability of success for a binomial variable.

We can specifiy different prior distributions for

*p*.

Let's assume that we have carried out 9 Bernoulli trials and observed 6 successes so far.

Let's try a uniform prior:

plot_posterior_using_uniform_prior <- function(num_points_in_grid)

{

` p_grid <- seq(from = 0, to=1, length.out=num_points_in_grid)`

` prior <- rep(1,num_points_in_grid)`

` likelihood <- dbinom(6, size=9, prob=p_grid)`

` ustd.posterior <- likelihood * prior`

` posterior <- ustd.posterior/sum(ustd.posterior)`

` plot(p_grid, posterior, type="b", xlab="probability of success", ylab="posterior probability for p")`

`}`

```
```

`Do a plot with a grid with 30 values for p:`

`plot_posterior_using_uniform_prior(num_points_in_grid=30) `

```
```

Let's try 10 values for p: plot_posterior_using_uniform_prior(num_points_in_grid=10)

`This doesn't give us such a good estimate of the posterior distribution for p:`

```
```

`What if we had a different prior, that was a step function, assigning 0 probability to all values of p less than 0.5:`

```
```

plot_posterior_using_step_prior <- function(num_points_in_grid){

p_grid <- seq(from = 0, to=1, length.out=num_points_in_grid)

prior <- ifelse ( p_grid < 0.5, 0, 1)

` likelihood <- dbinom(6, size=9, prob=p_grid) `

` ustd.posterior <- likelihood * prior `

` posterior <- ustd.posterior/sum(ustd.posterior) `

` plot(p_grid, posterior, type="b", xlab="probability of success", ylab="posterior probability for p")`

}`plot_posterior_using_step_prior(num_points_in_grid=30)`

` `

` `

**Quadratic approximation for estimating a posterior distribution**

Chapter 2 of Statistical Rethinking also explains that the grid approximation approach doesn't scale well and is slow, and that a faster, more scaleable approach is quadratic approximation. McElreath shows how it can be used to estimate the posterior distribution for the same parameter

*p*, as above: (see his book for how to install the 'rethinking' R package):

`library(rethinking)`

globe.qa <- quap(alist(W ~ dbinom(W+L,p), p ~ dunif(0,1)), data=list(W=6,L=3) + )

`precis(globe.qa)`

```
mean sd 5.5% 94.5%
p 0.67 0.16 0.42 0.92
```

` `

`This tells us that, using a uniform prior and assuming the posterior for `*p* is Gaussian (the key assumption of the quadratic approximation approach), its maximum (peak) is at 0.67 and its standard deviation is 0.16.

` `

`This is similar to the peak of the posterior distribution for `*p* that we estimated using the grid approximation approach above (when using a large number of grid points such as 30).

```
```

```
```

```
```

## No comments:

Post a Comment