机器算法验证 - 预测异方差数据的方差 - 吾爱随笔录

预测异方差数据的方差

机器算法验证回归 spss 方差残差异方差

2022-01-20 15:24:50

我正在尝试对异方差数据进行回归，其中我试图预测误差方差以及线性模型的平均值。像这样的东西：

\begin{aligned} y (x, t) & = \bar{y} (x, t) + ξ (x, t), \\ ξ (x, t) & \sim N (0, σ (x, t)), \\ \bar{y} (x, t) & = y_{0} + a x + b t, \\ σ (x, t) & = σ_{0} + c x + d t . \end{aligned}

$\begin{align}\\ y\left(x,t\right) &= \bar{y}\left(x,t\right)+\xi\left(x,t\right),\\ \xi\left(x,t\right) &\sim N\left(0,\sigma\left(x,t\right)\right),\\ \bar{y}\left(x,t\right) &= y_{0}+ax+bt,\\ \sigma\left(x,t\right) &= \sigma_{0}+cx+dt. \end{align}$

换句话说，数据由在 x 和 t 的不同值下对 ( ,的重复测量组成。我假设这些测量值由一个“真实”平均值组成，它是和的线性函数，以及加性高斯噪声，其标准偏差（或方差，我还没有决定）也线性地取决于。和进行更复杂的依赖——线性形式没有很强的理论动机——但我不想在这个阶段让事情变得过于复杂。） $y(x,t)$ $x$ $t$ $\bar{y}(x,t)$ $x$ $t$ $\xi(x,t)$ $x,t$ $x$ $t$

我知道这里的搜索词是“heteroscedasticity”，但到目前为止我所能找到的只是关于如何减少/删除它以更好地预测的讨论，但在尝试预测方面却一无所获在自变量方面。我想用置信区间（或贝叶斯等价物）估计和，如果在 SPSS 中有一种简单的方法就更好了！我该怎么办？谢谢。 $\bar{y}$ $\sigma$ $y_0, a, b, \sigma_0, c$ $d$

3个回答

我认为你的第一个问题是不再是正态分布，数据需要如何转换为同方差取决于究竟是什么是。例如，如果，则误差为比例型，回归前应取 y 数据的对数，或者从普通最小二乘法调整回归（OLS的加权最小二乘（将回归更改为最小化比例类型误差）。类似地，如果，则必须取对数的对数并将其回归。 $N\left(0,\sigma\left(x,t\right)\right)$ $\sigma\left(x,t\right)$ $\sigma\left(x,t\right)= ax+bt$ $1/y^2$ $\sigma\left(x,t\right)= e^{a x+b t}$

我认为对错误类型的预测很少被覆盖的原因是首先进行任何旧的回归（呻吟，通常是普通最小二乘法，OLS）。并且从残差图中，即，一个人观察残差形状，一个人绘制数据的频率直方图，并查看它。然后，如果残差是一个向右开口的扇形光束，则尝试比例数据建模，如果直方图看起来像指数衰减，则可以尝试倒数、等，以此类推，以求平方根、平方、取幂，取指数-y。 $model-y$ $1/y$

现在，这只是一个小故事。更长的版本包括更多类型的回归，包括泰尔中值回归、戴明二元回归和最小化不适定问题误差的回归，这些误差与最小化的传播误差没有特别的曲线拟合优度关系。最后一个是一个巨大的，但是，看这个举个例子。因此，人们试图获得的答案会产生很大的不同。通常，如果想要建立变量之间的关系，常规 OLS 不是首选方法，泰尔回归将是一种快速而肮脏的改进。OLS 仅在 y 方向上最小化，因此斜率太浅，截距太大而无法确定变量之间的潜在关系。换句话说，OLS 在给定 x 的情况下给出了 ay 的最小误差估计，它没有给出 x 如何随 y 变化的估计。当 r 值非常高（0.99999+）时，使用的回归几乎没有什么区别，y 中的 OLS 与 x 中的 OLS 大致相同，但是当 r 值很低时，y 中的 OLS 与x中的OLS。

总之，很大程度上取决于最初进行回归分析的动机是什么。这决定了所需的数值方法。做出选择后，残差就具有与回归目的相关的结构，需要在更大的背景下进行分析。

STATS BREUSCH PAGAN 扩展命令既可以测试残差的异方差性，也可以将其估计为部分或全部回归量的函数。

解决此类问题的一般方法是最大化数据的（正则化）可能性。

在您的情况下，对数似然看起来像

L L (y_{0}, a, b, σ_{0}, c, d) = \sum_{i = 1}^{n} \log ϕ (y_{i}, y_{0} + a x_{i} + b t_{i}, σ_{0} + c x_{i} + d t_{i})

$LL(y_0, a, b, \sigma_0, c, d) = \sum_{i=1}^n \log \phi(y_i, y_0 + a x_i + b t_i, \sigma_0 + c x_i + d t_i)$ 在哪里

ϕ (x, μ, σ) = \frac{1}{\sqrt{2 π} σ} e^{- \frac{(x - μ)^{2}}{2 σ^{2}}}

$\phi(x, \mu, \sigma) = \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(x-\mu)^2}{2\sigma^2}}$

您可以将此表达式编码到您最喜欢的统计包中的函数中（我更喜欢 Python、R 或 Stata，因为我从未在 SPSS 中进行过编程）。然后你可以将它提供给数值优化器，它会估计最优值 $\hat{\theta}$ 你的参数 $\theta=(y_0, a, b, \sigma_0, c, d)$ .

如果你需要置信区间，这个优化器还可以估计 Hessian 矩阵 $H$ 的 $\theta$ （二阶导数）在最优值附近。最大似然估计理论说，对于大 $n$ 的协方差矩阵 $\hat{\theta}$ 可以估计为 $H^{-1}$ .

这是 Python 中的示例代码：

import scipy
import numpy as np

# generate toy data for the problem
np.random.seed(1) # fix random seed
n = 1000 # fix problem size
x = np.random.normal(size=n)
t = np.random.normal(size=n)
mean = 1 + x * 2 + t * 3
std = 4 + x * 0.5 + t * 0.6
y = np.random.normal(size=n, loc=mean, scale=std)

# create negative log likelihood
def neg_log_lik(theta):
    est_mean = theta[0] + x * theta[1] + t * theta[2]
    est_std = np.maximum(theta[3] + x * theta[4] + t * theta[5], 1e-10)
    return -sum(scipy.stats.norm.logpdf(y, loc=est_mean, scale=est_std))

# maximize
initial = np.array([0,0,0,1,0,0])
result = scipy.optimize.minimize(neg_log_lik, initial)
# extract point estimation
param = result.x
print(param)
# extract standard error for confidence intervals
std_error = np.sqrt(np.diag(result.hess_inv))
print(std_error)

请注意，您的问题表述可能会产生负面影响 $\sigma$ ，我不得不通过蛮力替换太小来保护自己免受它的影响 $\sigma$ 和 $10^{-10}$ .

代码产生的结果（参数估计及其标准误差）是：

[ 0.8724218   1.75510897  2.87661843  3.88917283  0.63696726  0.5788625 ]
[ 0.15073344  0.07351353  0.09515104  0.08086239  0.08422978  0.0853192 ]

您可以看到估计值接近其真实值，这证实了此模拟的正确性。

其它你可能感兴趣的问题

上一篇UMP 测试何时不存在？下一篇多元生物时间序列：VAR和季节性