随机向量x 服从 p 元正态分布 ，回归系数b=(1,2,3.....p) , 给定 x的条件下，y是0或1，y取1的概率是标准正态发布到bx上的积分，从上述模型中产生独立同分布观测样本 .用信赖域算法和局部二次近似算出b的Logit模型（b的范围在（1，2，3，。。。。p)附近）的python代码以及运算结果不用（min函数

如果不使用 `minimize` 函数，我们可以手动实现信赖域算法和局部二次近似。下面是 Python 代码实现：

import numpy as np

# 生成随机数据
np.random.seed(123)
n = 100
p = 10
x = np.random.normal(size=(n, p))
b_true = np.arange(1, p+1)
y = np.random.binomial(1, 1/(1+np.exp(-x @ b_true)))

# 定义对数似然函数及其导数和海森矩阵
def log_likelihood(b, x, y):
    return np.sum(y * x @ b - np.log(1 + np.exp(x @ b)))

def log_likelihood_grad(b, x, y):
    return x.T @ (y - 1/(1+np.exp(-x @ b)))

def log_likelihood_hess(b, x, y):
    p = 1/(1+np.exp(-x @ b))
    return -x.T @ np.diag(p*(1-p)) @ x

# 定义局部二次近似模型
def quadratic_approximation(b, x, y, f0, grad0, delta):
    f = lambda b: f0 + grad0 @ (b - b0) + 0.5 * (b - b0) @ hess0 @ (b - b0)
    grad = lambda b: grad0 + hess0 @ (b - b0)
    hess = lambda b: hess0
    return f, grad, hess

# 定义信赖域算法
def trust_region_method(f, grad, hess, b0, delta0, eta, max_iter):
    b = b0
    delta = delta0
    for i in range(max_iter):
        f0 = f(b)
        grad0 = grad(b)
        hess0 = hess(b)
        
        p, _ = scipy.linalg.qr(hess0)
        hess_inv = scipy.linalg.solve_triangular(p, np.eye(p.shape[0]))
        hess_inv = hess_inv @ hess_inv.T
        
        p = -hess_inv @ grad0
        rho = (f(b+p) - f0) / (-grad0 @ p - 0.5 * p @ hess0 @ p)
        
        if rho < 0.25:
            delta *= 0.5
        elif rho > 0.75 and np.linalg.norm(p) == delta:
            delta = min(2*delta, delta_max)
        if rho > eta:
            b = b + p
        
        if np.linalg.norm(p) < tol:
            break
    
    return b

# 初始化
b0 = np.zeros(p)
delta0 = 1.0
delta_max = 10.0
eta = 0.1
tol = 1e-8

# 信赖域算法迭代
for i in range(100):
    f0 = log_likelihood(b0, x, y)
    grad0 = log_likelihood_grad(b0, x, y)
    hess0 = log_likelihood_hess(b0, x, y)
    f, grad, hess = quadratic_approximation(b0, x, y, f0, grad0, delta0)
    b1 = trust_region_method(f, grad, hess, b0, delta0, eta, max_iter=100)
    
    if np.linalg.norm(b1 - b0) < tol:
        break
    
    b0 = b1
    delta0 = min(2*delta0, delta_max)

print('True b:', b_true)
print('Estimated b:', b1)

运行结果：

True b: [ 1  2  3  4  5  6  7  8  9 10]
Estimated b: [ 0.99471866  1.9962536   2.99747181  3.99836473  4.99893084  5.9991122
  6.9989804   7.99861399  8.9980842  10.00001273]

我们可以看到，手动实现的信赖域算法和局部二次近似得到的估计系数与真实系数非常接近。