I've been trying to write a Python code to generate a 2 dimensional matrix of likelihood values for a grid of x̄ and σ values using 10 values as a set of measurements. x̄ goes from 0.5 to 1.5 in steps of 0.01 and σ goes from 0.01 to 0.3 in steps of 0.01. After that I want to create a contour plot of the distribution.

The likelihood function I'm using is Gaussian Likelihood. So far I have:

List1 = open("DataSet1.dat").readlines()
List2 = [float(i) for i in List1]
x = List2[0:10]

gaussian_pmf = lambda x, mu, sigma: (1 / (sigma*np.sqrt(2*np.pi)))*np.exp((-
(x-mu)**2)/(2*(sigma**2)))

mu = np.arange(0.5,1.51,0.01)
sigma = np.arange(0.01,1.02,0.01)

for i in range(0,len(mu)):
    distribution = []
    for i in range (0,len(x)):
        distribution.append(gaussian_pmf(x[i],mu,sigma))

plt.contour(distribution)
plt.show()

X is an array of the 10 data values I need but I'm not sure what I am doing wrong here. plt.contour requires a 2d array and that is what I'm supposed to create but I'm not quite sure how. Any help would be much appreciated!

solution

Your target function gaussian_pmf is 3-dimensional, so you'll have a separate contour plot for each mu, i.e. a hundred plots (is that really what you want?). Here's how you can do each one of them:

x = np.arange(0.5, 1.5, 0.01)
mu = np.arange(0.5, 1.51, 0.01)
sigma = np.arange(0.01, 1.02, 0.01)

# Create a meshgrid of (x, sigma)
X, Y = np.meshgrid(x, sigma)
Z = gaussian_pmf(X, mu[5], Y)

plt.contour(X, Y, Z)
plt.xlim([0.5, 0.7])
plt.ylim([0.0, 0.2])
plt.show()

Note that Z computes all values at once, since both X and Y are [101, 101] (meshed x and sigma).

The result plot:

You might also see these helpful examples.