Generate Random Draw from Gamma Distribution Re: Monte Carlo Simulation in R...
May be you are interested in the first `n' for which the sum of iid gamma rvs exceeds 2000, subject to the min-max constraints on each rv. If so, the following one-liner will give it to you: which(cumsum(pmax(1, pmin(rgamma(500, shape=0.067, rate=0.008), 85))) > 2000)[1] Note that I have used a slightly inefficient code (a somewhat generous N of 500), but this permits economy of code (on-liner!). Ravi. ---------------------------------------------------------------------------- ------- Ravi Varadhan, Ph.D. Assistant Professor, The Center on Aging and Health Division of Geriatric Medicine and Gerontology Johns Hopkins University Ph: (410) 502-2619 Fax: (410) 614-9625 Email: rvaradhan at jhmi.edu Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h tml ---------------------------------------------------------------------------- -------- -----Original Message----- From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf Of Hongwei Dong Sent: Tuesday, November 10, 2009 2:57 PM To: David Winsemius Cc: R-help Forum Subject: Re: [R] Generate Random Draw from Gamma Distribution Re: Monte Carlo Simulation in R... Thanks. I tried the rgamma function too. But I'm still wondering how I can set the min, max, and sum of the variates created by the random draws. Anyone has a clue? Thanks. Garry On Tue, Nov 10, 2009 at 11:47 AM, David Winsemius
<dwinsemius at comcast.net>wrote:
On Nov 10, 2009, at 2:26 PM, Hongwei Dong wrote: Exactly! Thanks, Duncan.
Let me re-phrase me question like this: 1) X_i values are independent Gammas, with the shape 0.067 and scale
0.008
2) Min(X)=1 and Max(X)=85
You might want to check that your parameterization in in agreement with that used by the rgamma function. Simply using those numbers yields a distribution that does not look as though it would get many qualifying samples. Here are 20 draws without any exclusions outside a range:
rgamma(20, shape=0.067, scale = 0.008)
[1] 2.213459e-03 2.815705e-05 2.381306e-04 2.264602e-07 1.293713e-07 7.680773e-38 6.441082e-15 6.168961e-13 [9] 5.089033e-06 1.571858e-16 9.869878e-12 1.813121e-13 1.253287e-11 1.852885e-04 4.212802e-07 1.774495e-25 [17] 1.892984e-07 5.927422e-17 1.322638e-12 4.327472e-05 http://finzi.psych.upenn.edu/R/Rhelp02/archive/31459.html 3) SUM(X)=2000
4) Do I also have to define the number of draws? if yes, it could be 250. Based on these restrictions, I want to generate random draw. I'm
wondering
how I can do this in R. Thanks. Garry On Tue, Nov 10, 2009 at 11:17 AM, Duncan Murdoch <murdoch at stats.uwo.ca
wrote:
On 11/10/2009 1:25 PM, Hongwei Dong wrote:
Hi, Dear R users,
I'm wondering if I can do Monte Carlo Simulation in R. My problem is like this: I know variable X follows Gamma distribution with shape parameter 0.067 and scale parameter 0.008. The sum of the X is 2000. I need R
help
me to simulate a vector of X that satisfies both the probability distribution and the sum. Anyone has a clue to this? Much appreciated.
Your requirements are slightly contradictory or incomplete. Here's one way to fully specify the problem: The X_i values are independent Gammas, with the given shape and scale. You want to simulate from the joint distribution conditional on the event sum(X) == 2000. Is that your problem? I don't know how to do the simulation, but maybe someone else does. Duncan Murdoch
[[alternative HTML version deleted]]
______________________________________________ R-help at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
David Winsemius, MD Heritage Laboratories West Hartford, CT
______________________________________________ R-help at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.