------------------------------------------------------------------------------- name: log: Q:\C-modelling\runmlwin\website\logfiles\2020-03-27\16\1_Introduct > ion_to_MCMC_Estimation_and_Bayesian_Modelling.smcl log type: smcl opened on: 27 Mar 2020, 17:43:30 . **************************************************************************** . * MLwiN MCMC Manual . * . * 1 Introduction to MCMC Estimation and Bayesian Modelling . . . . . . . 1 . * . * Browne, W. J. (2009). MCMC Estimation in MLwiN, v2.26. Centre for . * Multilevel Modelling, University of Bristol. . **************************************************************************** . * Stata do-file to replicate all analyses using runmlwin . * . * Chris Charlton and George Leckie, . * Centre for Multilevel Modelling, 2012 . * http://www.bristol.ac.uk/cmm/software/runmlwin/ . **************************************************************************** . . * 1.1 Bayesian modelling using Markov Chain Monte Carlo methods . . . . . .1 . . * 1.2 MCMC methods and Bayesian modelling . . . . . . . . . . . . . . . . .2 . . * 1.3 Default prior distributions . . . . . . . . . . . . . . . . . . . . .4 . . * 1.4 MCMC estimation . . . . . . . . . . . . . . . . . . . . . . . . . . .5 . . * 1.5 Gibbs sampling . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 . . * 1.6 Metropolis Hastings sampling . . . . . . . . . . . . . . . . . . . . 8 . . * 1.7 Running macros to perform Gibbs sampling and Metropolis . * Hastings sampling on the simple linear regression model . . . . . . 10 . . use "http://www.bristol.ac.uk/cmm/media/runmlwin/tutorial.dta", clear . . set seed 12345 . . generate y = normexam . . generate x = standlrt . . generate xsq = x^2 . . generate xy = x*y . . foreach var of varlist y x xsq xy { 2. . quietly summarize `var' 3. . local sum`var' = r(sum) 4. . } . . quietly count . . local N = r(N) . . local beta0 = 0 . . local beta1 = 0 . . local sigma2e = 1 . . local epsilon = 0.001 . . local burnin = 500 . . local chain = 5000 . . local totaliterations = `burnin' + `chain' . . local refresh = 50 . . if (`chain'>`N') set obs `chain' number of observations (_N) was 4,059, now 5,000 . . quietly gen iteration = . . . quietly gen beta0 = . . . quietly gen beta1 = . . . quietly gen sigma2e = . . . quietly gen e2i = . . . forvalues i = 1/`totaliterations' { 2. . local beta0 = (`sumy' - `beta1'*`sumx')/`N' + sqrt(`sigma2e'/`N')*inv > normal(uniform()) 3. . local beta1 = (`sumxy' - `beta0'*`sumx')/`sumxsq' + sqrt(`sigma2e'/`s > umxsq')*invnormal(uniform()) 4. . quietly replace e2i = (y - (`beta0' + (`beta1'*x)))^2 5. . quietly summarize e2i 6. . local sume2i = r(sum) 7. . local sigma2e = 1/((1/(`epsilon' + (`sume2i'/2)))*(invgammap(`epsilon > ' + (`N'/2),uniform()))) 8. . if (`i' > `burnin') { 9. . local c = `i' - `burnin' 10. . quietly replace iteration = `i' if _n==`c' 11. . quietly replace beta0 = `beta0' if _n==`c' 12. . quietly replace beta1 = `beta1' if _n==`c' 13. . quietly replace sigma2e = `sigma2e' if _n==`c' 14. . } 15. . if mod(`i',`refresh')==0 { 16. . display as text "Iteration `i'" 17. . } 18. . } Iteration 50 Iteration 100 Iteration 150 Iteration 200 Iteration 250 Iteration 300 Iteration 350 Iteration 400 Iteration 450 Iteration 500 Iteration 550 Iteration 600 Iteration 650 Iteration 700 Iteration 750 Iteration 800 Iteration 850 Iteration 900 Iteration 950 Iteration 1000 Iteration 1050 Iteration 1100 Iteration 1150 Iteration 1200 Iteration 1250 Iteration 1300 Iteration 1350 Iteration 1400 Iteration 1450 Iteration 1500 Iteration 1550 Iteration 1600 Iteration 1650 Iteration 1700 Iteration 1750 Iteration 1800 Iteration 1850 Iteration 1900 Iteration 1950 Iteration 2000 Iteration 2050 Iteration 2100 Iteration 2150 Iteration 2200 Iteration 2250 Iteration 2300 Iteration 2350 Iteration 2400 Iteration 2450 Iteration 2500 Iteration 2550 Iteration 2600 Iteration 2650 Iteration 2700 Iteration 2750 Iteration 2800 Iteration 2850 Iteration 2900 Iteration 2950 Iteration 3000 Iteration 3050 Iteration 3100 Iteration 3150 Iteration 3200 Iteration 3250 Iteration 3300 Iteration 3350 Iteration 3400 Iteration 3450 Iteration 3500 Iteration 3550 Iteration 3600 Iteration 3650 Iteration 3700 Iteration 3750 Iteration 3800 Iteration 3850 Iteration 3900 Iteration 3950 Iteration 4000 Iteration 4050 Iteration 4100 Iteration 4150 Iteration 4200 Iteration 4250 Iteration 4300 Iteration 4350 Iteration 4400 Iteration 4450 Iteration 4500 Iteration 4550 Iteration 4600 Iteration 4650 Iteration 4700 Iteration 4750 Iteration 4800 Iteration 4850 Iteration 4900 Iteration 4950 Iteration 5000 Iteration 5050 Iteration 5100 Iteration 5150 Iteration 5200 Iteration 5250 Iteration 5300 Iteration 5350 Iteration 5400 Iteration 5450 Iteration 5500 . . tabstat beta0 beta1 sigma2e, stats(mean sd) format(%4.3f) stats | beta0 beta1 sigma2e ---------+------------------------------ mean | -0.001 0.595 0.649 sd | 0.013 0.013 0.014 ---------------------------------------- . . . . * 1.8 Dynamic traces for MCMC . . . . . . . . . . . . . . . . . . . . . . 12 . . line beta0 iteration . . line beta1 iteration . . line sigma2e iteration . . . . * 1.9 Macro to run a hybrid Metropolis and Gibbs sampling method . * for a linear regression example . . . . . . . . . . . . . . . . . . 15 . . use "http://www.bristol.ac.uk/cmm/media/runmlwin/tutorial.dta", clear . . set seed 12345 . . generate y = normexam . . generate x = standlrt . . generate xsq = x^2 . . generate xy = x*y . . foreach var of varlist y x xsq xy { 2. . quietly summarize `var' 3. . local sum`var' = r(sum) 4. . } . . quietly count . . local N = r(N) . . local beta0 = 0 . . local beta1 = 0 . . local sigma2e = 1 . . local epsilon = 0.001 . . local burnin = 500 . . local chain = 5000 . . local beta0sd = 0.01 . . local beta1sd = 0.01 . . local beta0accept = 0 . . local beta1accept = 0 . . local totaliterations = `burnin' + `chain' . . local refresh = 50 . . if (`chain'>`N') set obs `chain' number of observations (_N) was 4,059, now 5,000 . . quietly gen iteration = . . . quietly gen beta0 = . . . quietly gen beta1 = . . . quietly gen sigma2e = . . . quietly gen e2i = . . . forvalues i = 1/`totaliterations' { 2. . * Update beta0 . * Calculate new proposed beta0 . local beta0prop = `beta0' + `beta0sd'*invnormal(uniform()) 3. . local beta0logpostdiff = -1*(2*(`beta0' - `beta0prop')*(`sumy' - `bet > a1'*`sumx') + `N'*((`beta0prop')^2 - (`beta0')^2))/(2*`sigma2e') 4. . if (`beta0logpostdiff' > 0) { 5. . * Definitely accept as higher posterior . local beta0 = `beta0prop' 6. . local beta0accept = `beta0accept' + 1 7. . } 8. . else { 9. . * Only sometimes accept . if (uniform() < exp(`beta0logpostdiff')) { 10. . local beta0 = `beta0prop' 11. . local beta0accept = `beta0accept' + 1 12. . } 13. . } 14. . * Update beta1 . local beta1prop = `beta1' + `beta1sd'*invnormal(uniform()) 15. . local beta1logpostdiff = -1*(2*(`beta1' - `beta1prop')*(`sumxy' - `be > ta0'*`sumx') + `sumxsq'*((`beta1prop')^2 - (`beta1')^2))/(2*`sigma2e') 16. . if (`beta1logpostdiff' > 0) { 17. . * Definitely accept as higher posterior . local beta1 = `beta1prop' 18. . local beta1accept = `beta1accept' + 1 19. . } 20. . else { 21. . * Only sometimes accept . if (uniform() < exp(`beta1logpostdiff')) { 22. . local beta1 = `beta1prop' 23. . local beta1accept = `beta1accept' + 1 24. . } 25. . } 26. . * Update sigma2e . quietly replace e2i = (y - (`beta0' + (`beta1'*x)))^2 27. . quietly summarize e2i 28. . local sume2i = r(sum) 29. . local sigma2e = 1/((1/(`epsilon' + (`sume2i'/2)))*(invgammap(`epsilon > ' + (`N'/2),uniform()))) 30. . if (`i' > `burnin') { 31. . local c = `i' - `burnin' 32. . quietly replace iteration = `i' if _n==`c' 33. . quietly replace beta0 = `beta0' if _n==`c' 34. . quietly replace beta1 = `beta1' if _n==`c' 35. . quietly replace sigma2e = `sigma2e' if _n==`c' 36. . } 37. . if mod(`i',`refresh')==0 { 38. . display as text "Iteration `i'" 39. . } 40. . } Iteration 50 Iteration 100 Iteration 150 Iteration 200 Iteration 250 Iteration 300 Iteration 350 Iteration 400 Iteration 450 Iteration 500 Iteration 550 Iteration 600 Iteration 650 Iteration 700 Iteration 750 Iteration 800 Iteration 850 Iteration 900 Iteration 950 Iteration 1000 Iteration 1050 Iteration 1100 Iteration 1150 Iteration 1200 Iteration 1250 Iteration 1300 Iteration 1350 Iteration 1400 Iteration 1450 Iteration 1500 Iteration 1550 Iteration 1600 Iteration 1650 Iteration 1700 Iteration 1750 Iteration 1800 Iteration 1850 Iteration 1900 Iteration 1950 Iteration 2000 Iteration 2050 Iteration 2100 Iteration 2150 Iteration 2200 Iteration 2250 Iteration 2300 Iteration 2350 Iteration 2400 Iteration 2450 Iteration 2500 Iteration 2550 Iteration 2600 Iteration 2650 Iteration 2700 Iteration 2750 Iteration 2800 Iteration 2850 Iteration 2900 Iteration 2950 Iteration 3000 Iteration 3050 Iteration 3100 Iteration 3150 Iteration 3200 Iteration 3250 Iteration 3300 Iteration 3350 Iteration 3400 Iteration 3450 Iteration 3500 Iteration 3550 Iteration 3600 Iteration 3650 Iteration 3700 Iteration 3750 Iteration 3800 Iteration 3850 Iteration 3900 Iteration 3950 Iteration 4000 Iteration 4050 Iteration 4100 Iteration 4150 Iteration 4200 Iteration 4250 Iteration 4300 Iteration 4350 Iteration 4400 Iteration 4450 Iteration 4500 Iteration 4550 Iteration 4600 Iteration 4650 Iteration 4700 Iteration 4750 Iteration 4800 Iteration 4850 Iteration 4900 Iteration 4950 Iteration 5000 Iteration 5050 Iteration 5100 Iteration 5150 Iteration 5200 Iteration 5250 Iteration 5300 Iteration 5350 Iteration 5400 Iteration 5450 Iteration 5500 . . tabstat beta0 beta1 sigma2e, stats(mean sd) format(%4.3f) stats | beta0 beta1 sigma2e ---------+------------------------------ mean | -0.001 0.596 0.649 sd | 0.012 0.013 0.014 ---------------------------------------- . . display as text "beta0 proposal accepted `=`beta0accept'/5500' times" beta0 proposal accepted .7625454545454545 times . . display as text "beta1 proposal accepted `=`beta1accept'/5500' times" beta1 proposal accepted .7501818181818182 times . . . . * 1.10 MCMC estimation of multilevel models in MLwiN . . . . . . . . . . .18 . . * Chapter learning outcomes . . . . . . . . . . . . . . . . . . . . . . . 19 . . . . . . **************************************************************************** . exit end of do-file