随机向量 x服从 p 元正态分布 ，回归系数b , 给定 的条件下，y是0,1，y等于1的概率是标准正态分布到bx的积分(iv)用信赖域算法和局部二次近似 编程实现b的最大似然估计 从上述模型中产生独立同分布观测样本 . python代码（不使用minize函数和optimistic包并且产生结果）其中b是（1，2，3.。。。p）附近

以下是使用 Python 实现的最大似然估计代码。代码中使用了 NumPy 库来进行数值计算。

import numpy as np

# 产生样本数据
p = 5 # 维度
n = 1000 # 样本数
x = np.random.normal(size=(n, p))
b_true = np.random.normal(size=p)
y_prob = np.exp(x.dot(b_true)) / (1 + np.exp(x.dot(b_true)))
y = np.random.binomial(1, y_prob)

# 定义似然函数和梯度
def log_likelihood(b):
    xb = x.dot(b)
    return np.sum(y * xb - np.log(1 + np.exp(xb)))

def gradient(b):
    return x.T.dot(y - 1 / (1 + np.exp(-x.dot(b))))

# 定义局部二次近似似然函数
def quadratic_approximation(b, delta):
    f = log_likelihood(b)
    g = gradient(b)
    H = x.T.dot(np.diag(y_prob * (1 - y_prob)).dot(x))
    p = np.linalg.solve(H, -g)
    q = 0.5 * p.dot(H).dot(p)
    if q <= delta:
        return f + p.dot(g) + q
    else:
        return np.nan

# 定义信赖域算法
def trust_region_method(b0, delta=0.1, eta=0.1, max_iter=100):
    b = b0.copy()
    f = log_likelihood(b)
    g = gradient(b)
    for i in range(max_iter):
        H = x.T.dot(np.diag(y_prob * (1 - y_prob)).dot(x))
        p = np.linalg.solve(H, -g)
        q = 0.5 * p.dot(H).dot(p)
        if q <= delta:
            b_new = b + p
            f_new = log_likelihood(b_new)
            rho = (f_new - f) / (-g.dot(p) - 0.5 * p.dot(H).dot(p))
            if rho > eta:
                b = b_new
                f = f_new
                g = gradient(b)
        delta = min(2 * delta, 1)
    return b

# 进行最大似然估计
b0 = np.ones(p)
b_mle = trust_region_method(b0)

# 输出结果
print("真实回归系数：", b_true)
print("估计回归系数：", b_mle)

代码中首先产生了一个样本数据，其中 x 是一个 p 维的正态分布随机向量，y 是一个二分类变量。然后定义了似然函数和梯度函数，以及局部二次近似似然函数。最后使用信赖域算法进行最大似然估计，并输出结果。