[R] rgamma function

[Thomas Lumley]

On Mon, Jul 16, 2012 at 12:49 PM, Chandler Zuo <zuo at stat.wisc.edu> wrote:
> Sorry that I posted the wrong syntax. My initial program is very long and I
> tried to post the section where I have been narrowed to locate the problem.
>
> In this sample I am simulating a Gamma with size parameter 5000 and scale
> parameter 1.
>
> I plugged in many breaks in my initial program and what I found was that
> most of the time the program stops when encountering a call to rgamma()
> function. It just freezes without popping out any error.
>
> #include<Rmath.h>
> #include<time.h>
> #include<Rinternals.h>
>
> SEXP generateGamma ()
> {
>     SEXP a;
>     PROTECT(a=allocVector(REALSXP,1));
>     srand(time(NULL));
>     REAL(a)[0]=rgamma(5000,1);
>     UNPROTECT(1);
>     return (a);
> }
>

srand() isn't relevant to the R random number generator, and you
haven't included the header file that defines it (this should lead to
a warning about implicit declaration).

More importantly, you haven't included the R equivalents GetRNGstate()
and PutRNGstate().  You want something like

#include "R.h"
#include<Rmath.h>
#include<Rinternals.h>

SEXP generateGamma ()
{
    SEXP a;
    PROTECT(a=allocVector(REALSXP,1));
    GetRNGstate();
    REAL(a)[0]=rgamma(5000,1);
    PutRNGstate();
    UNPROTECT(1);
    return (a);
}

    - thomas

-- 
Thomas Lumley
Professor of Biostatistics
University of Auckland