Stata The Stata listserver
[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

RE: st: Stata vs. SAS (random number generation) stata plugin help


From   "JP Azevedo" <jazevedo@provide.com.br>
To   <statalist@hsphsun2.harvard.edu>
Subject   RE: st: Stata vs. SAS (random number generation) stata plugin help
Date   Tue, 12 Apr 2005 00:27:45 -0300

Hello,

I would like to thank the list once again for their help.

Given that the procedure that I'm trying to replicate in Stata relies in
several sequences of random number (over 120,000), I believe that the
suggestion of generating the random number in one package and them exporting
to the other would not be very practical.

The path suggested by Joseph and Naji seem to be the most promising, that is
relying on a C or C++ pseudo-random number generator on both packages. Stata
has its plug-in capabilities and SAS has a toolkit called "user written
functions", which should work in a similar manner.

A friend has passed me a reference of the following pseudo-random number
algorithm in C which I'm trying to turn into a Stata plug-in:


long int kiss (long *i, long *j, long *k) { *j = *j ^ (*j << 17); *k = ( *k
^ (*k << 18) ) & 0X7FFFFFFF;

return ( (*i = 69069 * (*i) + 23606797) + (*j ^ = (*j >> 15)) +
         (*k ^ = (*k >> 13));

}

Christian P. Robert & George Casella (1999) Monte Carlo Statistical Methods.

Springer

Below you will find my first attempt to turn the code algorithm in a Stata
plug-in format. I'm a bit unsure how I should make Stata retrieve the
pseudo-random sequence generated by the code.

Any help will be much appreciated.

Cheers,

JP




*********************  start   ************************
*********************  kiss.c  ************************

/* KISS random number generator */

#include "stplugin.h"

STDLL stata_call(int argc, char *argv[])

    {

    ST_long         i, j, k, q ;
    ST_double       rnd ;
    ST_retcode      rc ;


    if(SF_nvars() < 3) {
         return(102) ;                            /* not enough variables
specified */
    }

    SF_scal_use(argv[0],&i) ; /* seed i */
    SF_scal_use(argv[1],&j) ; /* seed j */
    SF_scal_use(argv[2],&k) ; /* seed k */


    *j = *j ^ (*j << 17); *k = ( *k ^ (*k << 18) ) & 0X7FFFFFFF;

    return ( (*i = 69069 * (*i) + 23606797) + (*j ^ = (*j >> 15)) +
             (*k ^ = (*k >> 13)) );

    SF_vstore( ST_int(t), q, sum) ;

    return(0) ;

    }
*********************  end  ************************


-----Original Message-----
From: owner-statalist@hsphsun2.harvard.edu
[mailto:owner-statalist@hsphsun2.harvard.edu] On Behalf Of Naji
Sent: 08 April 2005 11:19
To: statalist@hsphsun2.harvard.edu
Subject: Re : st: Stata vs. SAS (random number generation)

Hi


Here is a link to MARSAGLIA Generator & Diehard battery test
I wonder if 
1) One can transpose C code in order to be used either in SAS/Stata (I don't
know how)
2) if the the same random seed will produce the same Pseudo Random Numbers


http://stat.fsu.edu/pub/diehard/

Hope it helps
Regards
Naji


Le 8/04/05 11:16, « Joseph Coveney » <jcoveney@bigplanet.com> a écrit :

> JP Azevedo wrote:
> 
> Many thanks you all for the valuable comments and suggestions.
> 
> I did a little bit of more digging, and found out that the references on
the
> "pseudo-random number generator" on the Stata 8 manuals and the SAS 5
> manuals are quite different.
> 
> While Stata talks about using George Marsalis (1994) 32-bits pseudo-random
> number generator, SAS actually has two different functions RANUNI [Fishman
> and Moore, 1982] and UNIFORM [Lewis, Goodman and Miller (1969) and Kennedy
> and Gentle (1980)].
> 
> Unfortunately the SAS manuals did not present a section with the full
> references, and the Stata one only said:
> 
> Marsaglia, G (1994). Personal communication
> 
> I would like to know if anyone familiar with this literature could advise
me
> how feasible it is to create a compatible pseudo-random number generator
> which would be compatible across these two packages?
> 
>
----------------------------------------------------------------------------
> 
> The Stata user manual does describe the algorithm in detail.  I guess that
> Stata's manual cites the personal communication, because George Marsaglia
> hadn't published his KISS generator at the time.  And apparently it still
> hasn't been formally published:  a recent article* cites his 1999 posting
to
> sci.stat.math (which provides C source code and can be found archived at
> www.cs.yorku.ca/~oz/marsaglia-rng.html among other places) to reference
the
> KISS generator.
> 
> SAS got slammed several years ago for its random number generator
> ( www.amstat.org/publications/tas/mccull-1.pdf and
> www.amstat.org/publications/tas/mccull.pdf ) primarily because of the
> inadequacy of its period length.  SAS Institute put out a rejoinder of
sorts
> ( support.sas.com/rnd/app/papers/statisticalaccuracy.pdf ).  (SAS's JMP
> offering has, however, switched to the Mersenne Twister.)  The Stata
manual
> states a period of ca. 2^126 for its generator.  Based upon this
> consideration, if you wanted to use the straightforward method that
Timothy
> McKenna suggests, then Stata -> SAS would be the better direction, unless
> SAS has since switched to a modern algorithm in its flagship product as
> well.
> 
> If you wanted a statistical-package-neutral source of numbers to use in
this
> way, Professor Marsaglia's DIEHARD site at Florida State University held
> downloadable sets of pseudorandom numbers.  That's gone, but there are
> numerous other Web sites that offer pseudorandom number sets for
> downloading, and a few with natural-origin random number sets.
> 
> Joseph Coveney
> 
> * Philip H.W. Leong, Guanglie Zhang, Dong-U Lee, Wayne Luk and John D.
> Villasenor, A comment on the implementation of the Ziggurat method.
_Journal
> of Statistical Software_ Volume 12, Issue 7, February 2005.
> www.jstatsoft.org/counter.php?id=114&url=v12/i07/v12i07.pdf&ct=1
> 
> 
> 
> 
> 
> 
> 
> 
> 
> *
> *   For searches and help try:
> *   http://www.stata.com/support/faqs/res/findit.html
> *   http://www.stata.com/support/statalist/faq
> *   http://www.ats.ucla.edu/stat/stata/
> 



*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/


*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/



© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index