利用newton迭代法求解三参数威布尔分布的极大似然估计值，python代码

利用牛顿迭代法求解三参数威布尔分布（Weibull distribution）的极大似然估计值，通常需要以下几个步骤：

首先，我们需要威布尔分布的概率密度函数 (PDF) 和其对数概率密度函数 (log PDF)，因为牛顿法优化的是目标函数的负梯度。

import numpy as np

def weibull_pdf(x, alpha, beta, loc=0):
    return (alpha * beta**alpha) / x**(alpha + 1) * np.exp(-(beta / x)**alpha)

def weibull_logpdf(x, alpha, beta, loc=0):
    return -np.log(alpha * beta**alpha) - (alpha + 1) * np.log(x) + (alpha * np.log(beta) + alpha * np.log(x))

# 假设我们有一个数据集 data
data = ...

# 初始化猜测值，一般选择数据的平均值作为初始点
initial_params = [np.mean(data), np.std(data), data.min()]

def newton_raphson_iterate(params, data, tolerance=1e-6, max_iterations=100):
    # 牛顿迭代的核心函数
    def gradient(params):
        alpha, beta, loc = params
        gradient_alpha = sum(-weibull_logpdf(x, alpha, beta, loc) for x in data)
        gradient_beta = sum(-((alpha + 1) * np.log(x) - alpha * np.log(beta)) for x in data)
        gradient_loc = sum(-np.log(x - loc) for x in data if x > loc)  # 去除小于loc的数据
        return gradient_alpha, gradient_beta, gradient_loc

    def hessian(params):
        alpha, beta, loc = params
        hessian_alpha = -sum(weibull_logpdf(x, alpha, beta, loc)**2 for x in data)
        hessian_beta = -sum(((alpha + 1) * np.log(x) - alpha * np.log(beta))**2 for x in data)
        hessian_loc = -sum(np.heaviside(loc - x, 0) * (-1 / (x - loc)**2) for x in data if x > loc)  # 去除小于loc的数据
        return hessian_alpha, hessian_beta, hessian_loc

    params_history = []
    current_params = params.copy()
    for _ in range(max_iterations):
        gradient_val = gradient(current_params)
        hessian_val = hessian(current_params)
        step = np.linalg.solve(hessian_val, -gradient_val)
        next_params = current_params + step
        params_history.append(next_params)

        if abs(sum((current_params - next_params)**2)) < tolerance:
            break
        else:
            current_params = next_params

    return next_params, params_history

estimated_params, convergence_curve = newton_raphson_iterate(initial_params, data)