<div dir="ltr"><div>Dear Statistical rethinkers, </div><div>I've thought a bit more on code 2.8 and I want to share it with you. Correct me please if I am wrong.  </div><div><br></div><div>the aim of the code: reconstruction of posterior probability distribution of proportion of the Earth surface covered with water given 6 waters out of 9 observations. This will be done by repeatedly generating random numbers under the light of observation (6 W out of 9). Each random number depends on the previous random number hence it is a Markov chain. </div><div><br></div><div># You decide on the size of the Markov chain. </div><font color="#0000ff">n_samples <- 1000</font><div><br></div><div># Memory allocation. Array p contains our Markov chain. <br><font color="#0000ff">p <- rep( NA , n_samples )</font></div><div><br></div><div># Setting the first number. It could be any valid value as long as you have enough random numbers. But let's be unbiased and start with 0.5. <br><font color="#0000ff">p[1] <- 0.5</font></div><div><br></div><div># These are our observations. <br><font color="#0000ff">W <- 6<br>L <- 3</font></div><div><br></div><div># Through this loop, in each step we generate a random number and if it is verified it will be stored in array p. You start with 2 because you already set the first one.  <br><font color="#0000ff">for ( i in 2:n_samples ) {</font></div><div><br></div><div># You generate a random number from normal dist with the mean equals to your previous (old) random number and call it p_new. You can use any other distribution (I think) for generating random numbers as long as it is dependent on the previous random number. <br><font color="#0000ff">p_new <- rnorm( 1 , p[i-1] , 0.1 )</font></div><div><br></div><div># Since the nature of p_new is probability, it must be between 0 and 1. Here you check for this and if it is outside this interval you bring it back inside. This code however is not perfect. Suppose you got 8, then based on the second checkpoint, you subtract it from 2 which you get -6 which is still outside the valid interval. You could use a command line like this instead:</div><div>while(p_new <0 | p_new >1)  p_new <- rnorm( 1 , p[i-1] , 0.1 ). <br><font color="#0000ff">if ( p_new < 0 ) p_new <- abs( p_new )<br>if ( p_new > 1 ) p_new <- 2 - p_new</font></div><div><br></div><div># Now you have a valid new random number and you need to decide whether to accept it or stay with the old random number. This will be done based on likelihoods of the old random number (q0) and new one (q1). The likelihoods are calculated under the light of observations. Here is the meaning of q0 for example (suppose the old value is still 0.5): If 50% of the Earth surface is covered by water, what is the probability that out of 9 sampling you get 6 waters. <br><font color="#0000ff">q0 <- dbinom( W , W+L , p[i-1] )<br>q1 <- dbinom( W , W+L , p_new )</font></div><div><br></div><div># Once you have q0 and q1, you accept or reject the new random number proportional to q1/q0. If for example q1/q0 happens to be 3, it means that the new random number is 3 times more plausible than the old one but this doesn't mean that you definitely accept the new one. You must always make some room for less plausible value as well, otherwise you don't get a distribution. <br><font color="#0000ff">p[i] <- ifelse( runif(1) < q1/q0 , p_new , p[i-1] )</font><br>}<br></div><div><br></div><div>PS1: Antonella noticed that upon increasing the size of the Markov chain, you get a better approximation to the analytical solution. </div><div><br></div><div>PS2: If you are not sure that you started the Markov chain with an unbiased value, it's better to discard some of the first random numbers and then convert it to the plot. BEAST discards the first 10000 numbers by default I think. Belinda! Is that correct? </div><div><br></div><div>This is my understanding. Hope it helps a bit. </div><div><br></div><div>All the best </div><div>Hassan </div></div>