-------------------------------------------------------------------------------
      name:  <unnamed>
       log:  Q:\C-modelling\runmlwin\website\logfiles\2024-10-11\18\16_An_Intro
> duction_to_Simulation_Methods_of_Estimation.smcl
  log type:  smcl
 opened on:  11 Oct 2024, 16:45:15

. ****************************************************************************
. *     MLwiN User Manual 
. *
. * 16  An Introduction to Simulation Methods of Estimation                241
. *
. *     Rasbash, J., Steele, F., Browne, W. J. and Goldstein, H. (2012).
. *     A User’s Guide to MLwiN, v2.26. Centre for Multilevel Modelling,
. *     University of Bristol.
. ****************************************************************************
. *     Stata do-file to replicate all analyses using runmlwin
. *
. *     George Leckie and Chris Charlton, 
. *     Centre for Multilevel Modelling, 2012
. *     https://www.bristol.ac.uk/cmm/software/runmlwin/
. ****************************************************************************
. 
. * 16.1 An illustration of parameter estimation with Normally distributed . . 
. * . . .data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .242
. 
. use "https://www.bristol.ac.uk/cmm/media/runmlwin/height.dta", clear

. 
. summarize height

    Variable |        Obs        Mean    Std. dev.       Min        Max
-------------+---------------------------------------------------------
      height |        100      175.35     10.0024        154        204

. 
. histogram height
(bin=10, start=154, width=5)

. 
. quietly summarize height

. 
. display 1 - normal((200 - r(mean))/r(sd))
.00686197

. 
. capture program drop heightsim1

. 
. program define heightsim1, rclass
  1. 
.         drop _all
  2. 
.         set obs 100
  3. 
.         generate heightsim = 175.35 + 10.002*invnorm(uniform())
  4. 
.         summarize heightsim
  5. 
.         return scalar pmean = r(mean)
  6. 
.         return scalar pvar  = r(Var)
  7. end

. 
. simulate, reps(1000) seed(12345): heightsim1

      Command: heightsim1
         pvar: r(pvar)
        pmean: r(pmean)

Simulations (1,000): .........10.........20.........30.........40.........50...
> ......60.........70.........80.........90.........100.........110.........120
> .........130.........140.........150.........160.........170.........180.....
> ....190.........200.........210.........220.........230.........240.........2
> 50.........260.........270.........280.........290.........300.........310...
> ......320.........330.........340.........350.........360.........370........
> .380.........390.........400.........410.........420.........430.........440.
> ........450.........460.........470.........480.........490.........500......
> ...510.........520.........530.........540.........550.........560.........57
> 0.........580.........590.........600.........610.........620.........630....
> .....640.........650.........660.........670.........680.........690.........
> 700.........710.........720.........730.........740.........750.........760..
> .......770.........780.........790.........800.........810.........820.......
> ..830.........840.........850.........860.........870.........880.........890
> .........900.........910.........920.........930.........940.........950.....
> ....960.........970.........980.........990.........1,000 done

. // Note: To obtain estimates as close as possible to the manual, ncrease the
. // number of reps to 10000.
. 
. generate iteration = _n

. 
. line pmean iteration

. 
. kdensity pmean

. 
. centile pmean, centile(2.5 97.5)

                                                          Binom. interp.   
    Variable |       Obs  Percentile    Centile        [95% conf. interval]
-------------+-------------------------------------------------------------
       pmean |     1,000        2.5    173.2063        173.0833    173.4998
             |                 97.5    177.2738        177.1252    177.4803

. 
. line pvar iteration

. 
. kdensity pvar

. 
. centile pvar, centile(2.5 97.5)

                                                          Binom. interp.   
    Variable |       Obs  Percentile    Centile        [95% conf. interval]
-------------+-------------------------------------------------------------
        pvar |     1,000        2.5     72.8153        70.54017    74.12077
             |                 97.5    128.0812        126.7303    132.4153

. 
. 
. 
. use "https://www.bristol.ac.uk/cmm/media/runmlwin/height.dta", clear

. 
. save "height.dta", replace
(file height.dta not found)
file height.dta saved

. 
. capture program drop heightsim2

. 
. program define heightsim2, rclass
  1. 
.         use "height.dta", clear
  2. 
.         bsample
  3. 
.         summarize height
  4. 
.         return scalar npmean = r(mean)
  5. 
.         return scalar npvar  = r(Var)
  6. end

. 
. simulate, reps(1000) seed(12345): heightsim2

      Command: heightsim2
        npvar: r(npvar)
       npmean: r(npmean)

Simulations (1,000): .........10.........20.........30.........40.........50...
> ......60.........70.........80.........90.........100.........110.........120
> .........130.........140.........150.........160.........170.........180.....
> ....190.........200.........210.........220.........230.........240.........2
> 50.........260.........270.........280.........290.........300.........310...
> ......320.........330.........340.........350.........360.........370........
> .380.........390.........400.........410.........420.........430.........440.
> ........450.........460.........470.........480.........490.........500......
> ...510.........520.........530.........540.........550.........560.........57
> 0.........580.........590.........600.........610.........620.........630....
> .....640.........650.........660.........670.........680.........690.........
> 700.........710.........720.........730.........740.........750.........760..
> .......770.........780.........790.........800.........810.........820.......
> ..830.........840.........850.........860.........870.........880.........890
> .........900.........910.........920.........930.........940.........950.....
> ....960.........970.........980.........990.........1,000 done

. // Note: To obtain estimates as close as possible to the manual, increase the
. // number of reps to 10000.
. 
. erase "height.dta"

. 
. generate iteration = _n

. 
. line npmean iteration

. 
. kdensity npmean

. 
. centile npmean, centile(2.5 97.5)

                                                          Binom. interp.   
    Variable |       Obs  Percentile    Centile        [95% conf. interval]
-------------+-------------------------------------------------------------
      npmean |     1,000        2.5     173.391        173.2132      173.59
             |                 97.5      177.29        177.1536    177.4074

. 
. line npvar iteration

. 
. kdensity npvar

. 
. centile npvar, centile(2.5 97.5)

                                                          Binom. interp.   
    Variable |       Obs  Percentile    Centile        [95% conf. interval]
-------------+-------------------------------------------------------------
       npvar |     1,000        2.5    73.62636        70.95973    76.35469
             |                 97.5    127.5858         124.643    129.8879

. 
. 
. 
. * 16.2 Generating random numbers in MLwiN . . . . . . . . . . . . . . . .249
. 
. clear

. 
. set obs 100
Number of observations (_N) was 0, now 100.

. 
. generate female = (uniform()<=0.6)

. 
. generate height2 = (1 - female)*(175 + 10*invnorm(uniform())) + female*(160 +
>  8*invnorm(uniform()))

. 
. histogram height2, by(female)

. 
. 
. 
. * Chapter learning outcomes . . . . . . . . . . . . . . . . . . . . . . .253
. 
. 
. 
. ****************************************************************************
. exit

end of do-file