How to efficiently generate a straight line with random slope and intercept in Python?

Question

Consider a very basic Monte Carlo simulation of a straight line y = m * x + b, e.g. To visualize the effect of uncertainty in the parameters m and b. m and b are both sampled from a normal distribution. Coming from a MATLAB-background, I would write this as

import numpy as np
import matplotlib.pyplot as plt

x = np.arange(start=0, stop=5, step=0.1)

n_data = len(x)
n_rnd = 1000

m = np.random.normal(loc=1, scale=0.3, size=n_rnd) 
b = np.random.normal(loc=5, scale=0.3, size=n_rnd)

y = np.zeros((n_data, n_rnd))  # pre-allocate y

for realization in xrange(n_rnd):
    y[:,realization] = m[realization] * x + b[realization]

plt.plot(x, y, "k", alpha=0.05);

This does produce the desired output, but I kinda feel like there must be a more "Pythonic" way to do this. Am I wrong? If not, could anyone provide me with some code example for how to do this more efficiently?

To give an example what I am looking for: In MATLAB this could easily be written without the loop using bsxfun(). Is there something similar in Python, or maybe even a package for things like these?

Answer 1

You can use numpy array broadcasting to create your y array in one step as shown below.

import numpy as np
import matplotlib.pyplot as plt

x = np.arange(start=0, stop=5, step=0.1)

n_data = len(x)
n_rnd = 1000

m = np.random.normal(loc=1, scale=0.3, size=n_rnd) 
b = np.random.normal(loc=5, scale=0.3, size=n_rnd)

y = m * x[:, np.newaxis] + b

for val in y.transpose():
    plt.plot(x, val, alpha=0.05)

# Or without the iteration:
# plt.plot(x, y, alpha=0.05)

plt.show()

x[:, np.newaxis] forces x to become a column vector of shape (50, 1) as opposed to (50,) which means that the broadcasting works.

You can then iterate over the numpy array directly (rather than iterating over it's index) but you have to transpose the array (using y.transpose()) otherwise for each iteration you'll get the x-value for each 1000 random numbers.