Re: [Rd] RNG stuck via Fortran call

From: Duncan Murdoch <murdoch_at_stats.uwo.ca>
Date: Wed 30 Nov 2005 - 15:50:53 GMT

On 11/30/2005 10:06 AM, Gilles GUILLOT wrote:
> Having not much success with my previous question I try to reformulate it:
>
> I'm simulating a Markow chain in Fortran interfaced with R.
> Each loop of my Fortran calls various functions of the R RNG through
> the wrapper given below.
>
> In a run of 100 iterations of the Markov kernel,
> after 20 iterations, the RNG seems like frozen.
>
> For example, the first call to the RNG in my loop is:
>
> rpostlamb = ggrgam(dble(m),1.d0)
> write(*,*) 'rpostlamb=',rpostlamb
>
>
> after the 20th iteration, it will return the same value
> rpostlamb= 1.24634557
>
> for all the following interations.
>
> Is there any suggestion of explanation for this strange fact ?

You need to put together code that someone else could compile and run, or show us your real code if it's short enough. It sounds as though you're not handling the RNG state properly.

Please don't send this code to me, post it to the list. (Actually, in putting together a short reproducible example, you're very likely to discover the bug yourself, and won't need to post anything.)

Duncan Murdoch
>
>
> Thanks in advance
>
> Gilles
>
> R Version 2.2.0 compiled under Mandrake 10.1
>
>
> #include <R.h>
> #include <Rmath.h>
>
> /******************/
> /* random numbers */
> /******************/
>
> void F77_SUB(rndstart)(void)
> { GetRNGstate(); }
>
> void F77_SUB(rndend)(void)
> { PutRNGstate(); }
>
> double F77_SUB(ggrnorm)(double *mu, double *sigma)
> { return rnorm(*mu, *sigma); }
>
> double F77_SUB(ggrexp)(double *scale)
> {return rexp(*scale);}
>
> double F77_SUB(ggrgam)(double *a, double *scale)
> {return rgamma(*a, *scale);}
>
> double F77_SUB(ggrunif)(double *a, double *b)
> {return runif(*a, *b);}
>
> double F77_SUB(ggrbinom)(double *n, double *p)
> {return rbinom(*n, *p);}
>
> double F77_SUB(ggrpois)(double *lambda)
> {return rpois(*lambda);}
>
> ______________________________________________
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel



R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Thu Dec 01 03:42:04 2005

This archive was generated by hypermail 2.1.8 : Wed 30 Nov 2005 - 21:25:24 GMT