Re: st: peculiar issue generating truncated normals

[Nick Cox <njcoxstata@gmail.com>]

This does not particularly surprise me given what you are doing, but
that is a visceral reaction. You seem to be going for the tails of the
distribution.

A different issue is that getting results one by one does not seem
essential: you can get a vector of results all at once and -- when
appropriate -- reject unwanted values all at once. Mata supports
elementwise operations with vectors and matrices.

On Fri, Jan 13, 2012 at 4:03 AM, Patrick Roland
<patrick.rolande@gmail.com> 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/