weibull分布的em代码

以下是使用EM算法拟合Weibull分布的代码：

import numpy as np
from scipy.stats import weibull_min

# 生成模拟数据
np.random.seed(42)
n_samples = 1000
true_shape = 1.5
true_scale = 2.0
y = weibull_min.rvs(c=true_shape, scale=true_scale, size=n_samples)

# EM算法拟合Weibull分布
shape = 1.0
scale = np.mean(y)
for i in range(100):
    # E步骤
    log_likelihoods = weibull_min.logpdf(y, c=shape, scale=scale)
    weights = np.exp(log_likelihoods)
    weights /= weights.sum()
    
    # M步骤
    scale = np.sum(weights * y) / np.sum(weights)
    shape = np.sum(weights * np.log(y / scale)) / np.sum(weights)
    
    print(f"Iteration {i+1}: shape={shape:.4f}, scale={scale:.4f}")
    
# 输出拟合参数
print(f"True shape={true_shape:.4f}, estimated shape={shape:.4f}")
print(f"True scale={true_scale:.4f}, estimated scale={scale:.4f}")

在上面的代码中，我们首先生成了1000个来自Weibull分布的模拟数据。然后使用EM算法拟合Weibull分布，其中E步骤计算每个样本属于每个分布的概率权重，M步骤使用加权最小二乘法来更新Weibull分布的形状和尺度参数。最后，我们输出了真实参数与估计参数的比较结果。