<html><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8"></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space;" class="">dear hassan,<div class=""><br class=""></div><div class="">thank you very much for this detailed description of the code. for each line of code, i will copy and paste your corresponding description in my exercise script.</div><div class=""><br class=""></div><div class="">i also get your point about the perfect not being perfect in keeping the probability within the (0,1) boundaries, but i have not exactly understood your solution. Maybe I will briefly bother you on Monady about that. But my more general question about this point would be, how likely it is to sample a number greater than 2. maybe the reason why McElreath used ‘2' as a threshold in the code is because ‘2' is a safe one (that is big enough)?</div><div class=""><br class=""></div><div class="">if you think something might not be perfect with the code, you can also trying writing directly to the author (<a href="mailto:richard_mcelreath@eva.mpg.de" class="">richard_mcelreath@eva.mpg.de</a>) . I had a small complain about something he has written in chapeter 3 and sent him an email, and he responded. This does not mean that he will always respond, but if your solution might improve his code, it would be in his interest to do so.</div><div class=""><br class=""></div><div class="">thanks again and a good week-end to all</div><div class=""><br class=""></div><div class="">antonella</div><div class=""><br class=""></div><div class=""><div><br class=""><blockquote type="cite" class=""><div class="">On 11. Feb 2022, at 17:19, Hassan Shafiey <<a href="mailto:h.shafiey@gmail.com" class="">h.shafiey@gmail.com</a>> wrote:</div><br class="Apple-interchange-newline"><div class=""><div dir="ltr" class=""><div class="">Dear Statistical rethinkers, </div><div class="">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 class=""><br class=""></div><div class="">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 class=""><br class=""></div><div class=""># You decide on the size of the Markov chain. </div><font color="#0000ff" class="">n_samples <- 1000</font><div class=""><br class=""></div><div class=""># Memory allocation. Array p contains our Markov chain. <br class=""><font color="#0000ff" class="">p <- rep( NA , n_samples )</font></div><div class=""><br class=""></div><div class=""># 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 class=""><font color="#0000ff" class="">p[1] <- 0.5</font></div><div class=""><br class=""></div><div class=""># These are our observations. <br class=""><font color="#0000ff" class="">W <- 6<br class="">L <- 3</font></div><div class=""><br class=""></div><div class=""># 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 class=""><font color="#0000ff" class="">for ( i in 2:n_samples ) {</font></div><div class=""><br class=""></div><div class=""># 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 class=""><font color="#0000ff" class="">p_new <- rnorm( 1 , p[i-1] , 0.1 )</font></div><div class=""><br class=""></div><div class=""># 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 class="">while(p_new <0 | p_new >1) p_new <- rnorm( 1 , p[i-1] , 0.1 ). <br class=""><font color="#0000ff" class="">if ( p_new < 0 ) p_new <- abs( p_new )<br class="">if ( p_new > 1 ) p_new <- 2 - p_new</font></div><div class=""><br class=""></div><div class=""># 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 class=""><font color="#0000ff" class="">q0 <- dbinom( W , W+L , p[i-1] )<br class="">q1 <- dbinom( W , W+L , p_new )</font></div><div class=""><br class=""></div><div class=""># 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 class=""><font color="#0000ff" class="">p[i] <- ifelse( runif(1) < q1/q0 , p_new , p[i-1] )</font><br class="">}<br class=""></div><div class=""><br class=""></div><div class="">PS1: Antonella noticed that upon increasing the size of the Markov chain, you get a better approximation to the analytical solution. </div><div class=""><br class=""></div><div class="">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 class=""><br class=""></div><div class="">This is my understanding. Hope it helps a bit. </div><div class=""><br class=""></div><div class="">All the best </div><div class="">Hassan </div></div>
_______________________________________________<br class="">Mitarbeiter.zoologie mailing list<br class=""><a href="mailto:Mitarbeiter.zoologie@lists.uni-halle.de" class="">Mitarbeiter.zoologie@lists.uni-halle.de</a><br class="">https://lists.uni-halle.de/cgi-bin/mailman/listinfo/mitarbeiter.zoologie<br class=""></div></blockquote></div><br class=""></div></body></html>