If you try to understand how Metropolis-Hastings MCMC works in producing a posterior for a parameter. You can read various website like this

But unless you understand R very well, it would make your head spin. But you can understand it easily with this simple julia code. You would need of course to understand some simple Bayesian theory beforehand.

```
using Distributions
using StatsPlots
# We perform 9 coin tosses and obtain 6 heads and 3 tails.
# What is the probability density function for obtaining a Head?
heads = 6
tails = 3
# define the number of points used (Not strictly required for MCMC)
numofpoints = 1001
# define probability grid from 0.0 to 1.0 with 1001 points
# (Not strictly required for MCMC)
p_grid = [ k for k in range(0,stop=1,length=numofpoints) ]
p_grid_lb = p_grid[1] # lower bound ( required for MCMC )
p_grid_ub = p_grid[end] # upper bound ( required for MCMC )
# define our prior which is a uniform prior ( required for MCMC )
prior = [ 1.0 for _ in 1:numofpoints ]
# define the binomial distribution
function dbinom(n::Integer,size::Integer,p::Real)
binomial(size, n) * p^n * (1.0-p)^(size-n)
end
# define the number of samples used in MCMC
# And create array p with that many elements (filled with 0.0)
n_samples = 1_000_000
p = zeros(n_samples)
# Set the starting value to a random number between 0.0 and 1.0
# strictly speaking it should be between p_grid_lb and p_grid_ub
p[1] = rand()
JumpingWidth = 0.1
for i in 2:n_samples
# p_new has a value around the vicinity of p[i-1]
p_new = rand( Normal( p[i-1] , JumpingWidth ) )
# make sure p_new is between 0.0 and 1.0
if p_new < p_grid_lb # if p_new < 0.0
p_new = p_grid_lb + abs( p_new - p_grid_lb ) # abs( p_new )
elseif p_new > p_grid_ub # if p_new > 1.0
p_new = p_grid_ub - abs( p_new - p_grid_ub ) # 2.0 - p_new
end
# Calc the probability of obtaining Heads and Tails
# assuming that probability = p[i-1] and probability = p_new
" q0 is posterior0 = likelihood0 * prior0 "
" q1 is posterior1 = likelihood1 * prior1 "
" In our cases, the prior is always 1.0 "
" or calc prior as prior[1+ceil(p*(numofpoints-1))] "
" or calc prior as Distributions.pdf(Distributions.Beta(1, 1),p) "
" where p is the probability "
" and the likelihood is dbinom( heads , heads+tails , probability ) "
q0 = dbinom( heads , heads+tails , p[i-1] ) * 1.0
q1 = dbinom( heads , heads+tails , p_new ) * 1.0
# The value of p[i] depends on whether the
# random number is less than q1/q0
p[i] = rand() < q1/q0 ? p_new : p[i-1]
end
# Next we need to display the histogram of the array p
# this is the shape of the posterior produced by MCMC
histogram(p,legend=false) |> display
```