Path: news.mathworks.com!not-for-mail
From: Peter Perkins <Peter.Perkins@MathRemoveThisWorks.com>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Precision issues with random numbers
Date: Thu, 09 Dec 2010 17:42:25 -0500
Organization: The MathWorks, Inc.
Lines: 40
Message-ID: <idrm0h$3ng$1@fred.mathworks.com>
References: <idr8ih$o1u$1@fred.mathworks.com> <idrdv5$5ej$1@fred.mathworks.com>
NNTP-Posting-Host: perkinsp.dhcp.mathworks.com
Mime-Version: 1.0
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 7bit
X-Trace: fred.mathworks.com 1291934545 3824 172.31.57.40 (9 Dec 2010 22:42:25 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Thu, 9 Dec 2010 22:42:25 +0000 (UTC)
User-Agent: Mozilla/5.0 (Windows; U; Windows NT 5.1; en-US; rv:1.9.2.8) Gecko/20100802 Thunderbird/3.1.2
In-Reply-To: <idrdv5$5ej$1@fred.mathworks.com>
Xref: news.mathworks.com comp.soft-sys.matlab:694062

On 12/9/2010 3:25 PM, Roger Stafford wrote:
> "Mark Salomon" <answer_through@thenewsgroup.thank.you.com> wrote in
> message <idr8ih$o1u$1@fred.mathworks.com>...
>> .......
>> Specifically, the function gamrnd returns zeros very often when the
>> shape parameter is low (e.g., less than 0.001). .......
> - - - - - - - -
> With a shape parameter that small you can expect a tight clustering of
> 'gamrnd' near zero. However the odds are very heavily against ever
> getting a zero.

Roger, for once I have to disagree with you, at least practical terms. 
With that shape parameter, there's something like a 50% chance of 
getting a value smaller than the smallest (nondenormal) double precision 
number.

 >> gamcdf(realmin,.001,1)
ans =
       0.49272

Mark, the problem is that gamrnd(.001,1,n,m) asks for numbers that can't 
be represented in d.p., and so the values are rounded to either realmin 
or zero.

 >> x = gamrnd(.001,1,1000000,1);
 >>
 >> sum(x <= realmin)
ans =
       492535

My question would be, is this really a physically reasonable 
probablility distribution?  If the answer is really yes, I think you'll 
need to generate log-gamma random values.  I don't know how to do that, 
but Luc Devroye's very fine book

    http://cg.scs.carleton.ca/~luc/rnbookindex.html

might be a place to start.

Hope this helps.