For Bayesian linear regression, we can write the model using the following equation:

\[P(\beta|y,X)=\frac{P(y,X|\beta)P(\beta)}{P(y,X)}\] where, \(P(y,X)=\int_{-\infty}^\infty P(y,X|\beta)P(\beta)d\beta\)

Usually, it is very difficult to calculate the unconditional probability \(P(y,X)\) i.e. the integral \(\int_{-\infty}^\infty P(y,X|\beta)P(\beta)d\beta\). One method to calculate the unconditional probability is to use Markov Chain Monte Carlo method.

Here, I use a very simple weather model from Wikipedia (https://en.wikipedia.org/wiki/Examples_of_Markov_chains) to show how to conduct the MCMC by hand then we will discuss the mathematics behind the MCMC.

Suppose, we want to calculate the unconditional probability of sunny(S) and Rainy(R) of weather, however, we only know the following four conditional probabilities:

\[P(S_{t+1}|R_t)=0.5 \tag{1}\] i.e given today’s weather is rainy, then tomorrow’s probability of sunny is 0.5 \[P(R_{t+1}|R_t)=0.5 \tag{2}\] i.e given today’s weather is rainy, then tomorrow’s probability of rainy is 0.5 \[P(R_{t+1}|S_t)=0.1 \tag{3}\] i.e given today’s weather is sunny, then tomorrow’s probability of rainy is 0.1 \[P(S_{t+1}|S_t)=0.9 \tag{4}\] i.e given today’s weather is sunny, then tomorrow’s probability of sunny is 0.9

And this four conditional probabilities satisfy the Markov property.

Now we design a random experiment for MCMC:

Step 1: Create Markov Chains

First chain:

Randomly select the first day, suppose it is a sunny day.We randomly select a number from 100 numbers i.e from 1 to 100, if this number is less or equal to 90 (note,we have 90% probability to select such numbers) we assign the next day is sunny because of condition \((4)\), for the third day, we continue to select a random number, suppose we picked 95 which is big than 90 (we have less than 10% probability to select such numbers from 1 to 100) therefore we assign the third day as rainy due to condition \((3)\), we continue to select the number for the fourth day such as if the random number greater than 50 we assign the fourth day as sunny and less than 50 we assign the fourth day as rainy because of \((1)\) or \((2)\), we can continue this process for a very long chain we may get the follow 15 days chain:

\[S-S-R-S-S-S-S-R-S-R-S-S-S-R-R\] Second chain:

we randomly select the first day again, suppose it was a rainy day, we continue the above procedures of random experiment we might the following 15 days chain:

\[R-S-S-S-R-S-R-S-S-S-S-S-S-S-S\] We continue above procedures we can create the third chain, fourth chain, and many other chains, the length of chains can be much longer.

step 2:

We sum the number of sunny days and rainy days of these chains, the frequency of these days are the unconditional probability of the sunny and rainy weather.

The following r code creates 2000 Markov chains and each chain includes 1000 days.

MCMC_sim<-function(n_sim,n_chain){
  
chain<-data.frame(mc="")

for(j in 1:n_sim){
  
chain1<-sample(LETTERS[18:19], 1, replace=TRUE)
p<-c()
for(i in 1:(n_chain-1)){
  p[i]=sample(1:100,1)

 
  if (substring(chain1,i,i)=='S' && p[i]<=90){chain1<-gsub(" ","", paste(chain1,"S"))}
  if (substring(chain1,i,i)=='S' && p[i]>90){chain1<-gsub(" ","", paste(chain1,"R"))}
  if (substring(chain1,i,i)=='R'){chain1<-gsub(" ","", paste(chain1,sample(LETTERS[18:19], 1, replace=TRUE)))}
}
chain<-rbind(chain, chain1)
}

chain$sum_R<-lengths(regmatches(chain$mc, gregexpr("R", chain$mc)))
p<-sum(chain$sum_R)/(n_sim*n_chain)
return(p)
}

MCMC_sim(2000,1000)
## [1] 0.167

We see the unconditional probability of rainy calculated by the Markov Chain Monte Carlo method is \(0.167\) which is the same as using the transition matrix method in Wikipedia (https://en.wikipedia.org/wiki/Examples_of_Markov_chains).

Next, I will review the mathematics behind the MCMC methods.