Re: st: peculiar issue generating truncated normals

[Patrick Roland <patrick.rolande@gmail.com>]

Yes, I'm still puzzled here. Any other ideas would be appreciated. If
this indeed related to the random number generators, it would seems
quite troublesome that they should be revealed to be inadequate in
such a simple problem.

On Sun, Jan 15, 2012 at 11:34 AM, Steve Samuels <sjsamuels@gmail.com> wrote:
>
> Patrick:
>
> I see that I misread the GHK reference.  It looks like your equations do match those shown there. Therefore I have no explanation for the discrepancy you've observed.
>
> Steve
>
> On Jan 13, 2012, at 8:32 PM, Steve Samuels wrote:
>
>
> In fact, for full generality, all the c[] terms should be terms of v = c*c'.  Since v = (1 1 /1 5),  in your example, c[1,1] = v[1,1] =1  and c[2,1] = v[2,1] = 1.
>
> Steve
>
> Or rather, the square root of the [2,2] element of c*c'.
>
>
>
>
> I can't be sure, but shouldn't c[2,2] in the r2 = GHK equation be sqrt(c[2,2])?
>
>
> Steve
> sjsamuels@gmail.com
>
> On Jan 12, 2012, at 11:03 PM, Patrick Roland wrote:
>
> I've been trying two ways of generating truncated multivariate normals
> in mata, which give different results. I'd be interested if anyone
> could shed light on why this might be. In the code I generate
> bivariate normal draws with variance c*c', where c = (1,0\1,2) ,
> truncated above at (-2,-2).
>
> In the first case, I use simple accept reject sampling. In the second,
> I sequentially draw truncated normals (this is explained here, for
> example: http://www.hss.caltech.edu/~mshum/gradio/ghk_desc.pdf).
>
> The means are consistently different, across many different seeds. The
> first variable has a mean of around -2.41 with accept reject and -2.37
> with sequential sampling. I'm running Stata SE 11.2. This seems
> peculiar to me - any ideas?
>
> mata
> u = (-2,-2)
> trials = 100000
> c = (1,0\1,2)
> GHK = J(trials,2,0)
> for(i=1;i<=trials;i++){
> r1 = invnormal(runiform(1,1)*normal(u[1]/c[1,1]))
> r2 = invnormal(runiform(1,1)*normal((u[2]-c[2,1]*r1)/c[2,2]))
> GHK[i,.] = (c*(r1\r2))'
> }
>
> i=0
> AR = J(trials,2,0)
> while(i<trials){
>        t = c*rnormal(2,1,0,1)
>        if(t'<u){
>                i=i+1
>                AR[i,.] = t'
>        }
> }
>
> mean(GHK)
> variance(GHK)
> mean(AR)
> variance(AR)
> end
> *
> *   For searches and help try:
> *   http://www.stata.com/help.cgi?search
> *   http://www.stata.com/support/statalist/faq
> *   http://www.ats.ucla.edu/stat/stata/
>
>
>
> *
> *   For searches and help try:
> *   http://www.stata.com/help.cgi?search
> *   http://www.stata.com/support/statalist/faq
> *   http://www.ats.ucla.edu/stat/stata/

*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/