-----------------------------------------------------------------------------------------------
      name:  <unnamed>
       log:  U:\projects\misc\twonormals\my2normals.txt
  log type:  text
 opened on:  29 Aug 2013, 15:46:45

. 
. clear all 

. set seed 1234567

. set obs 2000
obs was 0, now 2000

. 
. // Simulate data
. gen y1 = rnormal(5,2)

. gen y2 = rnormal(10,1)

. gen u = uniform() 

. gen y = cond(u<=0.3, y1, y2)

. label var y "N(5,2) with p1=0.3 & N(10,1) with p2=1-p1"

. hist y
(bin=33, start=-2.8955159, width=.47841464)

. 
. 
. // Define likelihood evaluator
. capture program drop twonormals_lf 

. program twonormals_lf
  1.    version 12
  2.         args todo b lnf 
  3.         tempvar xb1 xb2 lpr1 lnsigma1 lnsigma2 sigma1 sigma2 prob fxb1 fxb2 pr1 pr2
  4.         
.         // Define theta
.         mleval `xb1' = `b', eq(1)
  5.         mleval `xb2' = `b', eq(2)
  6.         mleval `lpr1' = `b', scalar eq(3) 
  7.         mleval `lnsigma1' = `b', scalar eq(4) 
  8.         mleval `lnsigma2' = `b', scalar eq(5)  
  9.         
.         // Transform so st deviations are > 0 and pri in [0,1]
.         qui gen double `sigma1' = exp(`lnsigma1')
 10.         qui gen double `sigma2' = exp(`lnsigma2')
 11.         qui gen double `pr1' = exp(`lpr1') / (1+exp(`lpr1'))
 12.         qui gen double `pr2' = 1 - `pr1'
 13.         
.         // Define likelihood
.         qui gen double `fxb1' = normalden($ML_y1,`xb1',`sigma1')
 14.         qui gen double `fxb2' = normalden($ML_y2,`xb2',`sigma2')
 15.         qui gen double `prob' = `pr1'*`fxb1' + `pr2'*`fxb2'
 16.         * Take log
.         mlsum `lnf' = ln(`prob')
 17.         
.         if (`todo'==0 | `lnf' >=.) exit         
 18. end

. 
. // Display all parameters
. capture program drop Replay

. program Replay, eclass
  1.         ml display
  2.         // pr1 and pr2
.         _diparm __lab__, label(    Derived   )
  3.         _diparm logitp1, invlogit label(p1)
  4.         ereturn scalar pi1_est = r(est)
  5.         ereturn scalar pi1_se = r(se)
  6.         _diparm logitp1, func(1-exp(@)/(1+exp(@))) ///
>                                         der(-exp(@)/((1+exp(@))^2)) label(p2)
  7.         ereturn scalar pi2_est = r(est)
  8.         ereturn scalar pi2_se = r(se)
  9.         // Sigma1 and sigma2
.         _diparm lnsigma1, exp label(sigma11)
 10.         ereturn scalar sigma1_est = r(est)
 11.         ereturn scalar sigma1_se = r(se)
 12.         _diparm lnsigma2, exp label(sigma2)
 13.         ereturn scalar sigma2_est = r(est)
 14.         ereturn scalar sigma2_se = r(se)        
 15.         _diparm __bot__
 16. end

. 
. // Fit null model
. *ml model d0 twonormals_lf (eq1: y=) (eq2: y=) /logitp1 /lnsigma1 /lnsigma2
. * can't evaluate without relatively good starting values
. * note that starting values are in order and on the esimation scale (logit for probabilities
. * and log(signa) for standard deviations
. 
. * give starting values
. ml model d0 twonormals_lf (eq1: y=) (eq2: y=) /logitp1 /lnsigma1 /lnsigma2, ///
>          init(5 5 0.2 0.5 0.6, copy)

. ml search
initial:       log likelihood = -9207.1602
rescale:       log likelihood = -5132.0744
rescale eq:    log likelihood = -4292.3953

. qui ml maximize, diff

. Replay

                                                  Number of obs   =       2000
                                                  Wald chi2(0)    =          .
Log likelihood = -4289.0173                       Prob > chi2     =          .

------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
eq1          |
       _cons |    10.0074   .0316402   316.29   0.000     9.945384    10.06941
-------------+----------------------------------------------------------------
eq2          |
       _cons |   4.978248   .1224625    40.65   0.000     4.738226     5.21827
-------------+----------------------------------------------------------------
logitp1      |
       _cons |   .7518246   .0609578    12.33   0.000     .6323496    .8712997
-------------+----------------------------------------------------------------
lnsigma1     |
       _cons |  -.0363745   .0250318    -1.45   0.146     -.085436    .0126869
-------------+----------------------------------------------------------------
lnsigma2     |
       _cons |   .6862525   .0455523    15.07   0.000     .5969716    .7755334
------------------------------------------------------------------------------
     Derived |            
          p1 |   .6795761   .0132737                       .653022    .7050161
          p2 |   .3204239   .0132737                      .2949839     .346978
     sigma11 |   .9642791   .0241376                      .9181119    1.012768
      sigma2 |   1.986258   .0904787                      1.816609     2.17175
------------------------------------------------------------------------------

. * matches true parameters
. 
. log close
      name:  <unnamed>
       log:  U:\projects\misc\twonormals\my2normals.txt
  log type:  text
 closed on:  29 Aug 2013, 15:46:48
-----------------------------------------------------------------------------------------------