随机向量 x服从 p元正态分布 ，回归系数b , 考虑如下的线性回归模型 y=bx+e , 其中随机误差项e 与x 相互独立，且e服从卡方（5），.从上述模型中产生独立同分布观测样本 . 在绝对值损失函数下建立中位数回归模型 (Median) (i) 用原对偶内点算法建立中位数回归的线性优化模型python代码以及运行结果（不用scipy包和pulp包）

中位数回归模型的优化问题可以表示为：

$$\min_{b}\sum_{i=1}^{n}|y_i - b^Tx_i|$$

其中，$x_i \in \mathbb{R}^p$ 是第 $i$ 个样本的特征向量，$y_i \in \mathbb{R}$ 是第 $i$ 个样本的标签，$b \in \mathbb{R}^p$ 是回归系数。

将绝对值函数改写成等价的线性约束形式：

$$ |y_i - b^Tx_i| \leq t_i $$

其中，$t_i \geq 0$ 是辅助变量。于是，中位数回归模型可以表示为以下线性规划问题：

$$ \min_{b,t} \sum_{i=1}^{n}t_i $$
$$ s.t. \ -t_i \leq y_i - b^Tx_i \leq t_i, \ i=1,\ldots,n $$

为了方便求解，可以将上述问题改写为对偶形式：

$$ \max_{\alpha} \sum_{i=1}^{n}y_i\alpha_i - b^T\sum_{i=1}^{n}x_i\alpha_i $$
$$ s.t. \ \sum_{i=1}^{n}\alpha_i = \frac{n}{2}, \ 0 \leq \alpha_i \leq \frac{1}{2}, \ i=1,\ldots,n $$

其中，$\alpha \in \mathbb{R}^n$ 是拉格朗日乘子。

现在，我们可以使用原对偶内点算法来求解上述线性规划问题。具体而言，我们需要先定义目标函数、约束条件和初始化参数：

import numpy as np

def objective_function(alpha, y):
    return -np.dot(y, alpha)

def constraint_function(alpha):
    return np.sum(alpha) - n / 2

def dual_gap(alpha, b, x, y, mu):
    t = y - np.dot(x, b)
    primal_value = np.sum(np.abs(t))
    dual_value = -np.dot(y, alpha) + np.dot(alpha, t) - 0.5 * np.dot(alpha, np.dot(x, x.T).dot(alpha))
    gap = primal_value - dual_value
    return np.abs(gap) / (1 + np.abs(primal_value) + np.abs(dual_value))

n = 1000  # 样本数
p = 10  # 特征维数
mu = 1e-3  # 内点法的参数
x = np.random.randn(n, p)  # 特征矩阵
y = np.random.chisquare(5, n)  # 标签向量
alpha = np.zeros(n)  # 拉格朗日乘子

然后，我们可以使用牛顿迭代法来更新 $b$ 和 $\alpha$：

def newton_step(alpha, b, x, y, mu):
    t = y - np.dot(x, b)
    w = 1 / (2 * t + mu)
    z = alpha + w * t
    q = w.sum()
    s = w.dot(y) - q * mu
    r = np.dot(x.T, (w * x).T)
    b = np.linalg.solve(r, x.T.dot(w * z))
    alpha = np.maximum(0, np.minimum(q * 0.5, z + s / mu))
    return alpha, b

def solve_dual(alpha, b, x, y, mu, max_iter=100, tol=1e-6):
    for i in range(max_iter):
        gap = dual_gap(alpha, b, x, y, mu)
        if gap < tol:
            break
        alpha, b = newton_step(alpha, b, x, y, mu)
    return alpha, b

alpha, b = solve_dual(alpha, np.zeros(p), x, y, mu)

最后，我们可以打印出结果：

print("Regression coefficients:", b)

完整的代码如下：