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

以下是使用信赖域算法和局部线性近似方法实现 b 的最大似然估计的 Python 代码，同时生成独立同分布观测样本：

import numpy as np
from scipy.stats import norm

np.random.seed(123)

# 生成数据
p = 5
n = 100
X = np.random.normal(size=(n,p))
beta = np.random.normal(size=p)
y = np.random.binomial(1, norm.cdf(X.dot(beta)))

# 定义负对数似然函数及其导数和海森矩阵
def neg_loglik(b, X, y):
    mu = norm.cdf(X.dot(b))
    return -np.sum(y*np.log(mu) + (1-y)*np.log(1-mu))

def grad_neg_loglik(b, X, y):
    mu = norm.cdf(X.dot(b))
    return X.T.dot(mu - y)

def hess_neg_loglik(b, X, y):
    mu = norm.cdf(X.dot(b))
    W = np.diag(mu*(1-mu))
    return X.T.dot(W).dot(X)

# 信赖域算法
delta = 0.1
eta = 0.1
b0 = np.zeros(p)
max_iter = 100
tol = 1e-6

for i in range(max_iter):
    negloglik = neg_loglik(b0, X, y)
    grad = grad_neg_loglik(b0, X, y)
    hess = hess_neg_loglik(b0, X, y)
    
    step = np.linalg.solve(hess, -grad)
    actual_reduction = neg_loglik(b0 + step, X, y) - negloglik
    predicted_reduction = -grad.dot(step) - 0.5*step.dot(hess).dot(step)
    rho = actual_reduction / predicted_reduction
    
    if rho < 0.25:
        delta = delta / 4
    elif rho > 0.75:
        delta = min(2*delta, 1e6)
    
    if rho > eta:
        b0 = b0 + step
        
    if np.linalg.norm(step) < tol:
        break

print("MLE of beta:", b0)

# 生成观测样本
X_new = np.random.normal(size=(n,p))
y_new = np.random.binomial(1, norm.cdf(X_new.dot(b0)))