Generating random variables in a joint normal distribution?
Hi all, Since I am fresh to use SciPy for simulation, I am not sure about how to get a sequence of pairs of random variables in a joint normal distribution. I have read the info docs in the module scipy.stats. It seems that only the normal distribution rvs for single variable is provided. Thus one sequence of random variables could be get by calling stats.norm.rvs() for times (right?). But if I want pairs of random variables, for example, (p, s), which are expected to be in the joint normal distribution of some given coefficient ro, would there be any routine provided to do this just like stats.norm.rvs? If currently there is no such routine, would there be some workaround to achieve this by combining the existing routines? I am not very good at numerical probabilities... Any suggestions would be much appreciated:) Regards, Parvel
Parvel Gu wrote:
Hi all,
Since I am fresh to use SciPy for simulation, I am not sure about how to get a sequence of pairs of random variables in a joint normal distribution.
I have read the info docs in the module scipy.stats. It seems that only the normal distribution rvs for single variable is provided. Thus one sequence of random variables could be get by calling stats.norm.rvs() for times (right?). But if I want pairs of random variables, for example, (p, s), which are expected to be in the joint normal distribution of some given coefficient ro, would there be any routine provided to do this just like stats.norm.rvs?
If currently there is no such routine, would there be some workaround to achieve this by combining the existing routines? I am not very good at numerical probabilities...
In [1]: from numpy import random In [2]: random.multivariate_normal? Type: builtin_function_or_method Base Class: <type 'builtin_function_or_method'> Namespace: Interactive Docstring: Return an array containing multivariate normally distributed random numbers with specified mean and covariance. multivariate_normal(mean, cov) -> random values multivariate_normal(mean, cov, [m, n, ...]) -> random values mean must be a 1 dimensional array. cov must be a square two dimensional array with the same number of rows and columns as mean has elements. The first form returns a single 1-D array containing a multivariate normal. The second form returns an array of shape (m, n, ..., cov.shape[0]). In this case, output[i,j,...,:] is a 1-D array containing a multivariate normal. -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
Hi, Thanks a lot. And I am still puzzled about the input arguments, mean and cov. So take my current problem for example. I am expecting the random variable P and S, which follow a joint normal distribution with (Mu)p=(Mu)s=0.5 (the mean?), and (Sigma)p=(Sigma)s=0.4 (the variance), and a coefficient ro = 0.8. According to the function multivariate_normal(mean, cov), only the matrixes of mean and cov are provided as input. Mapping to my problem, the mean could be [0.5, 0.5]. and the cov matrix is supposed to be [ cov(p,p), cov(p,s) cov(s,p), cov(s,s) ] Is it indicated that we have to get each cov(p, s) with some formula like ro = cov(p, s) / (sqrt(Dp) * sqrt(Ds)) = 0.8 then fill the result into the cov matrix? Regards, Parvel On 10/29/07, Robert Kern <robert.kern@gmail.com> wrote:
In [1]: from numpy import random
In [2]: random.multivariate_normal? Type: builtin_function_or_method Base Class: <type 'builtin_function_or_method'> Namespace: Interactive Docstring: Return an array containing multivariate normally distributed random numbers with specified mean and covariance.
multivariate_normal(mean, cov) -> random values multivariate_normal(mean, cov, [m, n, ...]) -> random values
mean must be a 1 dimensional array. cov must be a square two dimensional array with the same number of rows and columns as mean has elements.
The first form returns a single 1-D array containing a multivariate normal.
The second form returns an array of shape (m, n, ..., cov.shape[0]). In this case, output[i,j,...,:] is a 1-D array containing a multivariate normal.
-- Robert Kern
"I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco _______________________________________________ SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user
Parvel Gu wrote:
Hi,
Thanks a lot. And I am still puzzled about the input arguments, mean and cov.
So take my current problem for example. I am expecting the random variable P and S, which follow a joint normal distribution with (Mu)p=(Mu)s=0.5 (the mean?), and (Sigma)p=(Sigma)s=0.4 (the variance), and a coefficient ro = 0.8.
Careful there. The Greek letter sigma is usally reserved for the standard deviation, the square root of variance.
According to the function multivariate_normal(mean, cov), only the matrixes of mean and cov are provided as input. Mapping to my problem, the mean could be [0.5, 0.5]. and the cov matrix is supposed to be [ cov(p,p), cov(p,s) cov(s,p), cov(s,s) ]
Is it indicated that we have to get each cov(p, s) with some formula like ro = cov(p, s) / (sqrt(Dp) * sqrt(Ds)) = 0.8 then fill the result into the cov matrix?
Yes. -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
Hi, Thanks so much for the reply. I could go furthur in this problem. And here is the following issue: In my understanding, the continues calls to random.multivariate_normal would get a serials of random vairables which are all following the joint normal distribution. Then what about how to reset the random generating in each iteration? Assuming that I have to simulate something for 50 iterations then to get the average value. In each iteration I need a serial of random pairs which follows the joint normal distrubition. The interations are expected to be independent to each other. Thus: In iteration 1: serial = [] in sub loop serial.append(random.multivariate_normal(m, cov)) ... And in iteration 2, I want a fresh serial which should not be affected by the previous one s. Is there any problem to empty the serial then still call random.multivariate_normal(m, cov) ? In my understanding if there is no refresh or something performed to the internal random generator, the random.multivariate_normal would continue to generate variables in the context of the previous ones. Br, Parvel
Parvel Gu wrote:
Hi,
Thanks so much for the reply. I could go furthur in this problem.
And here is the following issue:
In my understanding, the continues calls to random.multivariate_normal would get a serials of random vairables which are all following the joint normal distribution. Then what about how to reset the random generating in each iteration?
Assuming that I have to simulate something for 50 iterations then to get the average value. In each iteration I need a serial of random pairs which follows the joint normal distrubition. The interations are expected to be independent to each other.
Thus:
In iteration 1: serial = [] in sub loop serial.append(random.multivariate_normal(m, cov)) ...
And in iteration 2, I want a fresh serial which should not be affected by the previous one s. Is there any problem to empty the serial then still call random.multivariate_normal(m, cov) ? In my understanding if there is no refresh or something performed to the internal random generator, the random.multivariate_normal would continue to generate variables in the context of the previous ones.
Just keep sampling. The samples will be as independent of each other as you can possibly get them. Reseeding the generator (the only thing I can imagine that you meant by "refresh or something") will only make things worse. -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
participants (2)
-
Parvel Gu -
Robert Kern