**long time ago…**, and kindly contributed to R-bloggers)

*n*,

*p*)

*n*is the number of tosses and the parameter

*p*is the probability to obtain a head. Given the number of tosses and the number of heads (

*h*), we can use Bayesian Inference to compute the probability to obtain a head (

*p*).

*p*given the number of heads (

*h*) is:

*h*|

*p*) is the likelihood of the data given a value of

*p*, Pr(

*p*) is the prior probability for

*p*and Pr(

*h*) is the marginal likelihood for the data. The prior probability of our parameter represent our believe about what values can take the parameter. If we think our parameter can take whatever value with the same probability, we use an uniform (flat) prior.

*p*and get the posterior probability of these values. To do this, we can use the Metropolis-Hastings MCMC algorithm, with the next steps (simplified):

*p.*

` p `

*p*, called

*p*‘.

` p_prime `

*p’*|

*h*) / Pr(

*p*|

*h*).

*h*|

*p*) can be defined in R language as:

` # help(dbinom) `

likelihood lh lh

}

*p*) is our prior probability, that is, our believe about what values can take

*p*. The only thing that we know is that it must be a value between 0 and 1, since it is a probability. If we think that all values have the same probability, we can define a flat prior using the

*beta*distribution. A

*beta*distribution with parameters β(1,1) is a flat distribution between 0 and 1 (you can learn more about

*beta*distribution here). Then, in R we can obtain Pr(

*p*) under a β(1,1) as:

` # help(dbeta) `

dbeta(p, 1, 1)

*p*. We express this equation in R language as follows:

` R `

*p*‘) and we update the value of

*p*=

*p*‘. If not, the change is rejected.

` if (R > 1) {R random if (random p } `

*p*.

` posterior[i,1] posterior[i,2] `

*p*. This can be easily done in R using a

*For*loop (check the full code below).

` # Set the numer of tosses. `

n # Set the number of heads obtained.

h # Define our likelihood function.

# Since our model is a binomial model, we can use:

likelihood lh lh

}

# Set the starting value of p

p # Create an empty data.frame to store the accepted p values for each iteration.

# Remember: "the posterior probability is just an updated version of the prior"

posterior # Set the lenght of the loop (Marcov Chain, number of iterations).

nrep # Start the loop (MCMC)

for (i in 1:nrep) {

# Obtain a new proposal value for p

p_prime # Avoid values out of the range 0 - 1

if (p_prime if (p_prime > 1) {p_prime # Compute the acceptance proability using our likelihood function and the

# beta(1,1) distribution as our prior probability.

R # Accept or reject the new value of p

if (R > 1) {R random if (random p }

# Store the likelihood of the accepted p and its value

posterior[i,1] posterior[i,2] print(i)

}

` par(mfrow= c(1,2)) `

prior plot(1:5000 ,posterior$V2, cex=0, xlab = "generations", ylab = "p",

main = "trace of MCMCn accepted values of parameter pn prior = beta(1,1) generations = 5000")

lines(1:5000, posterior$V2, cex=0)

abline(h=mean(posterior$V2), col="red")

plot(density(posterior$V2), xlim = c(min(min(prior),min((posterior$V2))), max(max(prior),max((posterior$V2)))),

ylim = c(0, max(max(density(prior)$y),max((density(posterior$V2)$y)))), main= "prior VS posteriorn prior= beta(1,1)",

lwd=3, col="red")

lines(density(prior), lwd=3, lty=2, col="blue")

legend("topleft", legend=c("prior density","posterior density"),

col=c("blue","red"), lty=c(3,1), lwd=c(3,3), cex = 1)

*p*…

Enjoy!

**leave a comment**for the author, please follow the link and comment on their blog:

**long time ago…**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more…