Notice: On April 23, 2014, Statalist moved from an email list to a forum, based at statalist.org.
From | "Ariel Linden, DrPH" <ariel.linden@gmail.com> |
To | <statalist@hsphsun2.harvard.edu> |
Subject | st: simulated value of a ratio does not compare to manually calculated value |
Date | Fri, 26 Apr 2013 16:28:23 -0400 |
Hi All, Below is a simulation that I am running in which I generate a ratio (mse_ratio) which is simply the ratio of the mse_true / mse_tt When running the simulation, the mse_ratio does not calculate correctly, even though each component does come through correctly... More specifically, you can see from the output below that the mse_true = .0038127, and the mse_tt = .0352618 , which should elicit a ratio of .10812551 However, the mse_ratio value is 91.79034. I cannot figure out what in my code is not working here... Any help is appreciated... See my code below... Thanks in advance! Ariel * this is the results of the simulation code from below . sum mse_true mse_tt mse_ratio Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------- mse_true | 50 .0038127 .0050518 9.49e-07 .0231653 mse_tt | 50 .0352618 .0266384 4.62e-07 .1037486 mse_ratio | 50 91.79034 619.6749 .0000267 4382.196 ********************************* * Start code ******************************* program drop _all program define sim, rclass drop _all set obs 1000 generate x=1 replace x=0 if _n>1000/2 generate ybase=invnorm(uniform()) generate e1=invnorm(uniform()) generate pi=(1-x)/(1+exp(-1+0.1*ybase)) + x/(1+exp(+0.9-0.1*ybase)) generate u1=uniform() gen a=0 replace a=1 if pi>=u1 gen ysqr=ybase^2 gen intxy=x*ybase gen intxy2=x*ysqr generate y=0.5*a-0.45+0.35*ybase+0.1*x+0.3*intxy+0.3*intxy2+0.25*ysqr+e1 * true estimates of Y regress y a x ybase intxy ysqr intxy2 scalar tx = _b[a] return scalar tx = tx scalar mse_true = (_b[a]-0.50)^2 return scalar mse_true = mse_true *ttest estimates - naïve ttest y, by(a) scalar ttest1 = r(mu_2)-r(mu_1) return scalar ttest1 = ttest1 scalar mse_tt = (ttest1-0.50)^2 return scalar mse_tt = mse_tt scalar mse_ratio = mse_true / mse_tt return scalar mse_ratio = mse_ratio end set seed 1234 simulate true=r(tx) ttest1=r(ttest1) mse_true=r(mse_true) /// mse_tt=r(mse_tt) mse_ratio=r(mse_ratio), reps(50): sim sum mse_true mse_tt mse_ratio * * For searches and help try: * http://www.stata.com/help.cgi?search * http://www.stata.com/support/faqs/resources/statalist-faq/ * http://www.ats.ucla.edu/stat/stata/