Skip to content

rgamma function

5 messages · Chandler Zuo, Peter Dalgaard, Thomas Lumley

#
Hi,

Has anyone encountered the problem of rgamma function in C? The 
following simplified program always dies for me, and I wonder if anyone 
can tell me the reason.

#include <Rmath.h>
#include <time.h>
#include <Rinternals.h>

SEXP generateGamma ()
{
     srand(time(NULL));
     return (rgamma(5000,1));
}

Has anyone encountered a similar problem before? Is there another way of 
generating Gamma random variable in C?

P.S. I have no problem compiling and loading this function in R.

Thanks for suggestions in advance!

--Chandler
1 day later
#
On Jul 14, 2012, at 04:55 , Chandler Zuo wrote:

            
It doesn't even give off a warning??

The prototype in Rmath.h is

double  rgamma(double, double);

and you should be returning an SEXP. As soon as something tries to interpret the double value as a pointer -- Poof!

Notice that rgamma in C is not the same function as the R counterpart, in particular it isn't vectorized, so only generates one random number at a time. The long and the short of it is that you need to read up on sections 5.9 and 5.10 of Writing R Extensions.

  
    
#
On Fri, Jul 13, 2012 at 7:55 PM, Chandler Zuo <zuo at stat.wisc.edu> wrote:
rgamma doesn't return an SEXP, it returns a double.  Also, the srand()
call is pointless.
Strange. You should get compiler warnings that the return type is
incompatible.  I get

foo.c: In function ?generateGamma?:
foo.c:7: warning: implicit declaration of function ?srand?
foo.c:8: error: incompatible types in return

I thought the ANSI standard actually *required* a diagnostic for the
incompatible return types.

   -thomas
#
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);
}


Thanks for your suggestions!
On 07/15/12 07:07, peter dalgaard wrote:
#
On Mon, Jul 16, 2012 at 12:49 PM, Chandler Zuo <zuo at stat.wisc.edu> wrote:
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