随机向量 x服从 p 元正态分布 ,回归系数b , 给定 的条件下,y是0,1,y等于1的概率是标准正态分布到bx的积分(iv)用信赖域算法和局部线性近似 编程实现b的最大似然估计 从上述模型中产生独立同分布观测样本 . python代码(不使用minize函数和optimistic包
时间: 2024-01-22 15:17:28 浏览: 41
以下是使用信赖域算法和局部线性近似实现b的最大似然估计的Python代码:
```python
import numpy as np
# 生成随机数据
np.random.seed(42)
n = 1000
p = 3
x = np.random.normal(size=(n, p))
b_true = np.array([0.5, -1.0, 2.0])
y_prob = np.exp(np.dot(x, b_true)) / (1 + np.exp(np.dot(x, b_true)))
y = np.random.binomial(1, y_prob)
# 定义似然函数
def log_likelihood(b, x, y):
y_prob = np.exp(np.dot(x, b)) / (1 + np.exp(np.dot(x, b)))
return np.sum(y*np.log(y_prob) + (1-y)*np.log(1-y_prob))
# 定义局部线性近似函数
def linear_approximation(b, x, y, h=0.1):
grad = np.zeros_like(b)
hess = np.zeros((len(b), len(b)))
ll = log_likelihood(b, x, y)
for i in range(len(b)):
b_plus = b.copy()
b_minus = b.copy()
b_plus[i] += h
b_minus[i] -= h
grad[i] = (log_likelihood(b_plus, x, y) - log_likelihood(b_minus, x, y)) / (2*h)
for j in range(i, len(b)):
b_pp = b.copy()
b_pm = b.copy()
b_mp = b.copy()
b_mm = b.copy()
b_pp[i] += h
b_pp[j] += h
b_pm[i] += h
b_pm[j] -= h
b_mp[i] -= h
b_mp[j] += h
b_mm[i] -= h
b_mm[j] -= h
hess[i, j] = (log_likelihood(b_pp, x, y) - log_likelihood(b_pm, x, y) -
log_likelihood(b_mp, x, y) + log_likelihood(b_mm, x, y)) / (4*h**2)
if i != j:
hess[j, i] = hess[i, j]
return ll, grad, hess
# 定义信赖域算法函数
def trust_region(x, y, initial_b, delta=1.0, eta=0.1, max_iter=100):
b = initial_b.copy()
ll, grad, hess = linear_approximation(b, x, y)
for i in range(max_iter):
# 在信赖域内求解二次规划问题
p = -np.linalg.solve(hess, grad)
pred_ll = log_likelihood(b+p, x, y)
actual_ll, _, _ = linear_approximation(b+p, x, y)
rho = (actual_ll - ll) / (pred_ll - ll)
if rho < 0.25:
delta *= 0.25
elif rho > 0.75 and np.abs(np.linalg.norm(p)-delta) < 1e-6:
delta = min(2*delta, 100)
if rho > eta:
b += p
ll = actual_ll
grad, hess = linear_approximation(b, x, y)
if np.linalg.norm(grad) < 1e-6:
break
return b
# 进行最大似然估计
b_init = np.zeros(p)
b_mle = trust_region(x, y, b_init)
print("最大似然估计的回归系数为:", b_mle)
```
在这段代码中,我们首先生成了随机数据,然后定义了似然函数和局部线性近似函数。接下来,我们使用信赖域算法函数`trust_region`来计算最大似然估计的回归系数。在每个迭代步骤中,我们首先使用当前的回归系数计算局部线性近似函数的梯度和海森矩阵。然后,我们在当前的信赖域内求解二次规划问题,以获得一个搜索方向。如果搜索方向指向了更好的解决方案,我们就将回归系数更新为搜索方向的值。最后,我们检查梯度是否已经足够小,如果是,则停止迭代。在这个例子中,我们使用了默认的参数`delta=1.0`和`eta=0.1`,并将最大迭代次数设置为100。最后输出的是最大似然估计的回归系数`b_mle`。
阅读全文