Notice: On April 23, 2014, Statalist moved from an email list to a forum, based at statalist.org.

# Re: st: peculiar issue generating truncated normals

 From Steve Samuels To statalist@hsphsun2.harvard.edu Subject Re: st: peculiar issue generating truncated normals Date Sun, 15 Jan 2012 14:34:24 -0500

```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

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/
```