Bookmark and Share

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


[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: st: peculiar issue generating truncated normals


From   Patrick Roland <[email protected]>
To   [email protected]
Subject   Re: st: peculiar issue generating truncated normals
Date   Fri, 13 Jan 2012 10:27:46 -0800

For the accept-reject sampler I am drawing a vector of results all at
once. The other procedure is inherently sequential.

This result surprises me because the random number generators in Stata
are meant to pass all kinds of tests for randomness. Yet it seems that
this is a test of randomness they don't pass...

On Fri, Jan 13, 2012 at 12:38 AM, Nick Cox <[email protected]> wrote:
> 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
> <[email protected]> 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/


© Copyright 1996–2018 StataCorp LLC   |   Terms of use   |   Privacy   |   Contact us   |   Site index