如何画出似然函数图像,蒙特卡洛仿真的图像,并加入克拉美罗下限,给出画图的代码
时间: 2024-10-19 15:11:22 浏览: 11
在Python中,我们可以使用matplotlib库来绘制似然函数、蒙特卡洛仿真图像以及克拉梅罗下限。以下是一个简单示例,假设我们有一个二维正态分布的数据点,并想基于该数据拟合一个参数(比如均值和方差)。
首先,我们需要导入必要的库并定义一些数据和函数:
```python
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal as mvn
# 定义数据点
data = np.random.multivariate_normal([0, 0], [[1, 0], [0, 1]], size=100)
# 定义似然函数
def log_likelihood(mean, cov, data):
return -0.5 * (len(data) * np.log(2*np.pi) + np.linalg.slogdet(cov)[1] + (data - mean).T @ np.linalg.inv(cov) @ (data - mean))
# 定义蒙特卡洛模拟
def monte_carlo(mean, cov, num_samples):
samples = np.random.multivariate_normal(mean, cov, num_samples)
return samples
# 克拉梅罗下限
def cramer_rao_bound(dim, n):
return dim / n
# 参数初始值
mean_initial = [0, 0]
cov_initial = np.eye(2)
# 绘制似然函数
mean_values = np.linspace(-3, 3, 100)
cov_values = np.linspace(0, 4, 100)
X, Y = np.meshgrid(mean_values, cov_values)
LL = np.zeros(X.shape)
for i in range(X.shape[0]):
for j in range(X.shape[1]):
LL[i,j] = log_likelihood([X[i,j], Y[i,j]], np.array([[cov_values[j], 0], [0, cov_values[j]]]), data)
plt.contourf(X, Y, LL, levels=np.arange(np.min(LL), np.max(LL)+1, 0.5))
plt.colorbar(label='Log-Likelihood')
plt.title("Likelihood Surface")
plt.xlabel("Mean x")
plt.ylabel("Mean y")
plt.show()
# 绘制蒙特卡洛仿真
num_samples = 1000
samples = monte_carlo(mean_initial, cov_initial, num_samples)
plt.scatter(samples[:, 0], samples[:, 1], alpha=0.5)
plt.title("Monte Carlo Simulation")
plt.xlabel("Sample Mean x")
plt.ylabel("Sample Mean y")
# 绘制克拉梅罗下限
crb = cramer_rao_bound(2, len(data))
plt.plot([-3, 3], [0]*2, 'k--', label=f"CRB (x): {crb}")
plt.plot([0]*2, [-3, 3], 'k--', label=f"CRB (y): {crb}")
plt.legend()
plt.show()
```
阅读全文