Bivariate Normal Distribution

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.colors import LinearSegmentedColormap

x = np.linspace(-3, 3, 200)
y = np.linspace(-3, 3, 200)
X, Y = np.meshgrid(x, y)

# Bivariate normal
mu_x, mu_y = 0, 0
sigma_x, sigma_y = 1, 1
rho = 0.5
Z = (1/(2*np.pi*sigma_x*sigma_y*np.sqrt(1-rho**2))) * np.exp(
    -1/(2*(1-rho**2)) * ((X-mu_x)**2/sigma_x**2 - 2*rho*(X-mu_x)*(Y-mu_y)/(sigma_x*sigma_y) + (Y-mu_y)**2/sigma_y**2)
)

fig, ax = plt.subplots(figsize=(10, 8), facecolor='#ffffff')
ax.set_facecolor('#ffffff')

# Custom colormap: coral tones
colors = ['#fef2f2', '#fca5a5', '#F5276C', '#9C2007']
cmap = LinearSegmentedColormap.from_list('coral', colors, N=256)

cs = ax.contourf(X, Y, Z, levels=25, cmap=cmap)
ax.contour(X, Y, Z, levels=12, colors='#7f1d1d', linewidths=0.5, alpha=0.6)
cbar = plt.colorbar(cs, ax=ax, pad=0.02)
cbar.set_label('Probability Density', color='#374151', fontsize=11)
cbar.ax.yaxis.set_tick_params(color='#374151')
plt.setp(plt.getp(cbar.ax.axes, 'yticklabels'), color='#374151')

ax.set_xlabel('X', fontsize=11, color='#374151', fontweight='500')
ax.set_ylabel('Y', fontsize=11, color='#374151', fontweight='500')
ax.set_title('Bivariate Normal Distribution', fontsize=14, color='#1f2937', fontweight='bold', pad=15)

ax.tick_params(colors='#6b7280', labelsize=9)
for spine in ax.spines.values():
    spine.set_color('#d1d5db')

plt.tight_layout()
plt.show()