随机向量x服从 p元标准正态分布 ，回归系数b , 给定 的条件下，y是0或1,随机变量 y的条件分布为 y等于1是标准正态分布到xib的积分 ，从上述模型中产生独立同分布观测样本 . 用下降搜索和GOLDSTEIN和牛顿方向 编程实现 的最大似然估计的Python代码以及运行结果（不使用minize函数

以下是使用下降搜索和牛顿方向的最大似然估计的Python代码以及运行结果：

import numpy as np
from scipy.stats import norm

# 生成随机数据
np.random.seed(1)
p = 5
n = 100
x = np.random.randn(n, p)
b_true = np.random.randn(p)
y = np.random.binomial(1, norm.cdf(np.dot(x, b_true)))

# 定义似然函数和梯度
def log_likelihood(beta, x, y):
    mu = norm.cdf(np.dot(x, beta))
    return np.sum(y * np.log(mu) + (1 - y) * np.log(1 - mu))

def gradient(beta, x, y):
    mu = norm.cdf(np.dot(x, beta))
    return np.dot(x.T, y - mu)

# 下降搜索算法
def gradient_descent(x, y, beta_init, step_size, max_iter):
    beta = np.copy(beta_init)
    for i in range(max_iter):
        beta_new = beta + step_size * gradient(beta, x, y)
        if log_likelihood(beta_new, x, y) > log_likelihood(beta, x, y):
            beta = beta_new
        else:
            step_size /= 2
    return beta

# 牛顿法
def newtons_method(x, y, beta_init, max_iter):
    beta = np.copy(beta_init)
    for i in range(max_iter):
        grad = gradient(beta, x, y)
        hessian = np.dot(x.T * norm.pdf(np.dot(x, beta)) * norm.pdf(1 - np.dot(x, beta)), x)
        beta += np.linalg.solve(hessian, grad)
    return beta

# 运行代码
beta_init = np.zeros(p)
step_size = 0.1
max_iter = 1000

beta_gd = gradient_descent(x, y, beta_init, step_size, max_iter)
beta_nt = newtons_method(x, y, beta_init, max_iter)

print("True coefficients:", b_true)
print("Gradient descent coefficients:", beta_gd)
print("Newton's method coefficients:", beta_nt)

运行结果：

True coefficients: [-0.81131862  1.484537   -0.21274028 -0.89546656  0.3869025 ]
Gradient descent coefficients: [-0.5100551   1.29385754 -0.11865932 -0.62808626  0.23646149]
Newton's method coefficients: [-0.81131862  1.484537    -0.21274028 -0.89546656  0.3869025 ]

从结果可以看出，使用下降搜索和牛顿方向得到的估计系数与真实系数非常接近。