I'm producing a multivariate probability density function below. This works fine but I'm hoping to normalise the Z value so it elicits a value between 0 and 1.
To achieve this I want to divide the distribution value at the mean so it's always 1 at the mean and lower elsewhere. I understand the sum of all values will be greater than 1.
I'm diving Z but the sum of Z but when printing the values, they still are outside my intended normalised range.
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
# Our 2-dimensional distribution will be over variables X and Y
N = 60
X = np.linspace(-3, 3, N)
Y = np.linspace(-3, 4, N)
X, Y = np.meshgrid(X, Y)
# Mean vector and covariance matrix
mu = np.array([0., 1.])
Sigma = np.array([[ 1. , -0.5], [-0.5, 1.5]])
# Pack X and Y into a single 3-dimensional array
pos = np.empty(X.shape + (2,))
pos[:, :, 0] = X
pos[:, :, 1] = Y
def multivariate_gaussian(pos, mu, Sigma):
"""Return the multivariate Gaussian distribution on array pos.
pos is an array constructed by packing the meshed arrays of variables
x_1, x_2, x_3, ..., x_k into its _last_ dimension.
"""
n = mu.shape[0]
Sigma_det = np.linalg.det(Sigma)
Sigma_inv = np.linalg.inv(Sigma)
N = np.sqrt((2*np.pi)**n * Sigma_det)
# This einsum call calculates (x-mu)T.Sigma-1.(x-mu) in a vectorized
# way across all the input variables.
fac = np.einsum('...k,kl,...l->...', pos-mu, Sigma_inv, pos-mu)
return np.exp(-fac / 2) / N
# The distribution on the variables X, Y packed into pos.
Z = multivariate_gaussian(pos, mu, Sigma)
#normalise Z so range is 0-1
Z = Z/sum(Z)
print(Z)
# Create a surface plot and projected filled contour plot under it.
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_surface(X, Y, Z, rstride=3, cstride=3, linewidth=1, antialiased=True,
cmap=cm.magma)
cset = ax.contourf(X, Y, Z, zdir='z', offset=-0.15, cmap=cm.magma)
# Adjust the limits, ticks and view angle
ax.set_zlim(-0.15,1)
ax.view_init(27, -21)
plt.show()
If you want to normalise Z, you need to divide it not by the sum, but by the maximum over all its values. Thus you ensure that the new maximum is 1:
# normalise Z so range is 0-1
Z = Z / np.max(Z)
# show summary statistics for the result
import pandas as pd
print(pd.Series(Z.flatten()).describe())
count 3.600000e+03
mean 1.605148e-01
std 2.351826e-01
min 6.184228e-08
25% 7.278911e-03
50% 4.492385e-02
75% 2.135538e-01
max 1.000000e+00
dtype: float64
Since you have a Gaussian distribution and you only changed the scaling, the maximum will still be at the mean x and y values. Note, however, that Z is now no longer a probability density function.