|
|
|
|
Subject: Random numbers .... suite
From: Jourdain, Damien (D.JOURDAIN@CGIAR.ORG)
Date: ven nov 30 2001 - 15:42:13 CET
Hi there,
I am still struggling with my correlated random numbers. It progressed
fairly well since my last message.
This is an update on where I am at the moment.
Part I : Normal random numbers.
First step is to create random numbers following a normal distribution from
random numbers following a random numbers following a uniform (0 to 1)
distributiuon as the Cormas random is doing. I found a rather elegant and
robust method over the internet. This is described at
www.taygeta.com/random/gaussian.html . I wrote a method randomNormal that I
inserted as a Cormas class method:
randomNormal
"genere des nombres aleatoires suivant une loi de distribution
normale N(0,1)"
"Ecrit par Damien Jourdain .... attention ne fait pas partie du core
certifie par Cormas !!"
" Source: www.taygeta.com/random/gaussian.html "
| x1 x2 w u y1 |
w := 2. "initialise w avec une valeur plus grande que 1 "
[w >= 1] whileTrue:
[x1 := (2.0 * Cormas random asFloat - 1.0) asFloat.
x2 := (2.0 * Cormas random asFloat - 1.0) asFloat.
w := x1 * x1 + (x2 * x2)].
u := (w ln * -2 / w) sqrt.
y1 := x1 * u.
^y1
Testing several times this method, it seems rather reasonable, and is not
consuming to much time .. so I kept it for the moment !
If anybody has something better to propose ... most welcome !
Part II: Correlated random numbers (theory)
Some of you have followed the out-of-forum discussion. I'll post here the
last message so you can keep track of what we are talking about.
Let say I want four variables to be correlated with the following
correlation matrix:
R = 1 0.4 0.8 0.3
0.4 1 0.3 0.2
0.8 0.3 1 0.2
0.3 0.2 0.2 1
Let see what Rudolf K. Bock say about correlated random numbers.
http://rkb.home.cern.ch/rkb/AN16pp/node232.html
"The notion of correlation is linked with that of variance and elements in
an error matrix. Correlated random numbers arise from uncorrelated random
numbers by error propagation. If correlated random numbers have to be
generated according to a known error matrix, the inverse operation (of error
propagation) is required: what is wanted is the matrix A which transforms
the unit matrix I into the (given) error matrix E when propagating errors.
This is exactly the problem of Cholesky decomposition: A will be a
triangular matrix, and it can be found from E by that comparatively simple
algorithm. "
For those who remember matrix algebra (!)
For a matrix A, the LU decomposition factorizes the matrix A = LU. If A
satisfies the above criteria, one can decompose more efficiently
into A = U'.U, where U (which can be seen as the "square root" of A) is a
upper triangular matrix with positive diagonal elements.
Our upper triangular matrix U (I did it with SAS, but I guess you could
easily do it in Excel) corresponding to our initial R matrix is:
U = 1 0.4 0.8 0.3
0 0.916 -0.021 0.087
0 0 0.599 -0.063
0 0 0 0.947
Now, life will almost be easy, we just have to draw four random numbers X Y
Z T corresponding to non correlated normal distributions N(0,1).
Each sample drawn is a (1,4) matrix, that we will postmultiply by the U
(4,4) matrix. We will therefore obtain a (1,4) matrix of the transformed
variables.
If we make repeated draws, you will find out after a while that those
variables are correlated as expected !
Once you have those four variables, you can adjust them to get another mean
and variance. For example
X2 = 100 + (X * 10) give a variable N(100,10)... Since linear
transformations do not affect correlation between variables, we have the
desired variables in terms of mean, variance and correlation between
variables!
Part III: Correlated random numbers : a proposition of method
I have written a method (again as a Cormas class method) called
randomCorrelated: . This method receives as argument the transpose of the
cholesky factorisation matrix (lower triangular matrix), so I can reproduce
the matrix product somehow easilier.
Here is the method. I welcome comments on how to improve it:
correlatedRandom: laMatrice
"this method receives as argument an OrderedCollection of OrderedCollection"
"that is the root of the correlation matrix (Cholesky) To ease calculations
I took the transpose of this matrix, and used a premultiplication"
| dimension randInit rand2 r xx randFin |
dimension := laMatrice size. "square matrix, so we need only one
dimension "
"genere les nombres aleatoires initiaux tires de distributions normales
reduites centrees"
randInit := OrderedCollection new.
1 to: dimension
do: [:a | randInit add: Cormas randomNormal ].
"create a new empty matrix called rand2, this matrix will receive the inner
product the matrix randInit and laMatrice"
rand2 := OrderedCollection new.
1 to: dimension
do: [:a | r := OrderedCollection new: dimension.
rand2 add: r].
1 to: dimension
do: [:i |
(laMatrice at: i) with: randInit do: [:x :y | (rand2
at: i) addLast: x * y]]. "inner product"
randFin := OrderedCollection new. "initialise la collection de sortie"
rand2 do:
[:col |
xx := col inject: 0 into: [:s :suivant | s + suivant]. "sum of the
lines of the rand2 "
randFin add: xx].
^randFin
This method gives as an output an OrderedCollection with the desired
correlated random numbers, and should normally (hopefully) work whatever the
size of the correlation matrix, i.e. it can generate more (or less) than 4
correlated series of random numbers.
If you want to try it, you may try it in a workspace as followed:
| c1 c2 c3 c4 uneMatrice |
c1 := #(1 0 0 0) asOrderedCollection.
"correspond to the matrix first column "
c2 := #(0.4 0.916 0 0 ) asOrderedCollection.
c3 := #(0.8 -0.021 0.599 0 ) asOrderedCollection.
c4 := #(0.3 0.087 -0.063 0.947) asOrderedCollection.
uneMatrice := OrderedCollection new.
uneMatrice add: c1 ; add: c2 ; add: c3 ; add: c4.
^ (Cormas correlatedRandom: uneMatrice).
I tried that several times with large samples (200), and it seems to
reproduce fairly well the correlations between variables.
Sorry for this long message. But since we started the discussion here, I
wanted to keep you informed on the subject.
Looking forward to get some feedback on the methods developed, and ways to
improve them. For example, is it OK to use OrderedCollection to reproduce
matrix product, or should I sue other kind of collection?
Also a problem is that the Cholesky factorisation of the correlation matrix
has to be calculated outside the program (I used SAS, but any statistical
package would do). I guess it is not too difficult to calculate it outside
the model, and introduce it as an argument. And I do not think any one wants
to enter into matrix algebra (at least not me!) to propose a method for
Cholesky factorization?
Cheers,
Damien