import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.stats import multivariate_normal
def Gaussian_Distribution(N=2, M=1000, m=0, sigma=1):
    mean = np.zeros(N) + m
    cov = np.eye(N) * sigma

    data = np.random.multivariate_normal(mean, cov, M)
    Gaussian = multivariate_normal(mean=mean, cov=cov)
    
    return data, Gaussian

M = 100
data, Gaussian = Gaussian_Distribution(N=2, M=M, sigma=0.1)
X, Y = np.meshgrid(np.linspace(-1,1,M), np.linspace(-1,1,M))
d = np.dstack([X,Y])
Z = Gaussian.pdf(d).reshape(M,M)


fig = plt.figure(figsize=(6,4))
ax = Axes3D(fig)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='seismic', alpha=0.8)

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.savefig('123.png')
plt.show()