<- file stat 97random.html ->
code for Normal simulation (1997)
<![CDATA[
This file has - code for Normal simulation, and
comments on testing of random number generators.
]]>
Code for Normal simulation
<![CDATA[
=======================Hans-Peter Piepho, 17 Apr 1997==========ssc
Message-ID: <9704171450.AA15263@fserv.wiz.uni-kassel.de>
From: Hans-Peter Piepho <piepho@WIZ.UNI-KASSEL.DE>
Subject: Re: code for creating simulated data
>
> I need to simulate a set of normally distributed variables so
>that 2 to 3 of them have a 'known' correlation that I am able to
>vary.. I trust this has been done many times before and I would
>appreciate any shared code or suggestions about where to go look for
>help as I start (I am a SAS user currently....though I have used
>BMDP and SPSS as well and can often figure out how to translate those
>codes).
> I am planning on using simulated data to better understand
>some estimate changes we are seeing in the Proc Mixed models we are
>running. Admiitedly, our real data have a 'correlated error
>structure' due to 2 repeated measures on our subjects. Though a
>large N and focus on quite large changes in estimates has made me believe
>I may not need to simulate the correlated errors to play with this
>initially.... (these simulations are simply to help our group learn and
>'argue' .... ). However, folks are welcome to argue otherwise and
>provide more complicated code for error structure simulation as well
>.....
>
You can do a Cholesky decomposition of your variance-covariance matrix and
use that to generate correlated data.
1. Decompose your variance-covariance matrix V as V=A'A
2. Generate 3 standard normal deviates and store them in vector x
3. y=mu + A'x has mean mu and variance-covariance matrix V.
The IML function root() yields the Cholesky decomposition of a matrix.
Hans-Peter
/*Example*/
proc iml;
n=100; /*number of samples*/
p=3; /*size of sample*/
seed=-1;
v={7 1 3,
1 8 2,
3 2 9}; /*your covariance matrix*/
mu={10, 20, 30}; /*your mean vector*/
a=root(v);
do i=1 to n;
x=normal(repeat(seed,p));
y=mu + a`*x;
print y;
yy=yy//y`;
end;
sum_yy=yy[+,];
cov_est=(yy`*yy-sum_yy`*sum_yy/n)/(n-1);
print cov_est;
quit;
_______________________________________________________________________
Hans-Peter Piepho
*--------
]]>
testing random generators
<![CDATA[
=======================Philip Staite, 03 Feb 1997==========ssm,...
From: Philip Staite <pstaite@vnet.ibm.com>
Subject: Re: Excellent pseudo random number generator
Message-ID: <32F5F657.167E@vnet.ibm.com>
Agner Fog wrote:
<edited>
> 6. Request for feedback
> -----------------------
> If there is any math genius out there who can refine my analysis then
> I would be grateful to hear from you. My heuristic arguments about
> cycle lengths rely on the assumption that the SHUFFLADD algorithm has
> no regularities which can reduce cycle lengths, such as the ADDGEN has.
> Proving this theoretically is beyond my powers. And an experimental
> proof is only possible for systems too small to be useful for the
> intended purposes.
>
> Also, if anybody can prove, theoretically or statistically, that the
> random numbers produced by my algorithm are less than perfectly random,
> then please let me know.
Try generating 1024 numbers from your algorithm then putting them into a
1024 pt FFT. Another test is to "oversample" or take every 4096th or
8192nd number generated and pass that into the 1024pt FFT. (or would
that be "undersample?")
Depending on seed value input you can get some tremendous frequency
correlations. ie. the output groups into a relatively small number of
"bins" whose magnitude dwarfs the rest. This seems to be a "feature" of
addative generators in general. (???) The addative 11//55 generator
from "Numerical Recipies in C" exhibits the same kind of behavior.
So far the best (with regard to lack of frequency correlations) method
I've seen is to use multiple linear-congruential generators, (with
differing constants) and hash together a subset of the mid-bits from
their output.
* * * * * * * * * * * * * * * * * * * * * * * * * * * * *
]]> Document by Rich Ulrich. E-mail to wpilib+@pitt.edu
FAQ top.
Ulrich home page.
Ulrich FAQ.
http://www.pitt.edu/~wpilib/stats99.html