5

I want to know how to generate the same random (Normal Distribution) numbers in numpy as I do in MATLAB.

As an example when I do this in MATLAB

RandStream.setGlobalStream(RandStream('mt19937ar','seed',1));
rand
ans =
    0.417022004702574

Now I can reproduce this with numpy:

import numpy as np
np.random.seed(1)
np.random.rand()
0.417022004702574

Which is nice, but when I do this with normal distribution I get different numbers.

RandStream.setGlobalStream(RandStream('mt19937ar','seed',1));
randn
ans =
    -0.649013765191241

And with numpy

import numpy as np
np.random.seed(1)
np.random.randn()
1.6243453636632417

Both functions say in their documentation that they draw from the standard normal distribution, yet give me different results. Any idea how I can adjust my python/numpy to get the same numbers as MATLAB.

Because someone marked this as a duplicate: This is about normal distribution, as I wrote in the beginning and end. As I wrote uniform distribution works fine, this is about normal distribution. None of the answers in the linked thread help with normal distribution.

  • 2
    One possible hack to get same normal distribution samples could be to generate samples from uniform distribution and then convert them into normal distribution: http://stackoverflow.com/questions/75677/converting-a-uniform-distribution-to-a-normal-distribution – Prakhar Jul 19 '16 at 14:47
  • See also: http://stackoverflow.com/questions/3722138/is-it-possible-to-reproduce-randn-of-matlab-with-numpy?rq=1 (It showed up in the "Related" sidebar; might as well make it "Linked".) – Warren Weckesser Jul 19 '16 at 15:57

1 Answers1

2

My guess would be that the matlab and numpy may use different methods to get normal distribution of random numbers (which are obtained from uniform numbers in some way).

You can avoid this problem by writing a box-muller method to generate the random numbers yourself. For python,

import numpy as np

# Box-muller normal distribution, note needs pairs of random numbers as input
def randn_from_rand(rand):

    assert rand.size == 2

    #Use box-muller to get normally distributed random numbers
    randn = np.zeros(2)
    randn[0] = np.sqrt(-2.*np.log(rand[0]))*np.cos(2*np.pi*rand[1])
    randn[1] = np.sqrt(-2.*np.log(rand[0]))*np.sin(2*np.pi*rand[1])

    return randn


np.random.seed(1)
r = np.random.rand(2)
print(r, randn_from_rand(r))

which gives,

(array([ 0.417022  ,  0.72032449]), array([-0.24517852, -1.29966152]))

and for matlab,

% Box-muller normal distribution, note needs pairs of random numbers as input
function randn = randn_from_rand(rand)

    %Use box-muller to get normally distributed random numbers
    randn(1) = sqrt(-2*log(rand(1)))*cos(2*pi*rand(2));
    randn(2) = sqrt(-2*log(rand(1)))*sin(2*pi*rand(2));

which we call with

RandStream.setGlobalStream(RandStream('mt19937ar','seed',1));
r = [rand, rand]
rn = randn_from_rand(r)

with answer,

r =
    0.4170    0.7203
rn =
   -0.2452   -1.2997

Note, you can check the output is normally distributed, for python,

import matplotlib.pyplot as plt

ra = []
np.random.seed(1)
for i in range(1000000):
    rand = np.random.rand(2)
    ra.append(randn_from_rand(rand))


plt.hist(np.array(ra).ravel(),100)
plt.show()

which gives,

enter image description here

Ed Smith
  • 12,716
  • 2
  • 43
  • 55
  • Thanks for the answer. But I can not change the MATLAB code. It looks like MATLAB is using the ziggurat, I will try to reimplement that in python. –  Jul 19 '16 at 16:02
  • 1
    I think unless your implementation is exactly the same, you may not get reproducible random numbers, according to this link https://github.com/numpy/numpy/issues/5810 the code for MATLAB is proprietary, which suggests it probably isn't simple. Could you not simply redefine randn in MATLAB and use the same code? – Ed Smith Jul 19 '16 at 16:08