Skip to content
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal

# Define the covariance matrix and mean vector
D = np.array([[4, 2],  # Covariance matrix
              [2, 3]])
mean = np.array([0, 0])  # Mean vector

# Generate random samples
n_samples = 1000
samples = np.random.multivariate_normal(mean, D, n_samples)

# Extract random effects: b_i1 (intercepts) and b_i2 (slopes)
b_i1 = samples[:, 0]
b_i2 = samples[:, 1]

# Plot the samples
plt.figure(figsize=(8, 6))
plt.scatter(b_i1, b_i2, alpha=0.5, edgecolor='k', label='Random Effects')
plt.axhline(0, color='gray', linestyle='--', linewidth=0.8)
plt.axvline(0, color='gray', linestyle='--', linewidth=0.8)
plt.xlabel('$b_{i1}$ (Random Intercept)')
plt.ylabel('$b_{i2}$ (Random Slope)')
plt.title('2D Normal Distribution of Random Effects')
plt.legend()
plt.grid(alpha=0.3)
plt.show()