机器算法验证 - 数据具有不确定性的线性模型，使用 R - 吾爱随笔录

数据具有不确定性的线性模型，使用 R

机器算法验证 r 最小二乘错误传播

2022-03-02 10:56:34

假设我有一些不确定的数据。例如：

例如，不确定性的性质可以是重复测量或实验，或测量仪器的不确定性。

我想用 R 拟合一条曲线，通常我会用lm. 但是，当它给我拟合系数的不确定性以及预测区间的不确定性时，这并没有考虑到数据的不确定性。查看文档，该lm页面具有以下内容：

...权重可用于表示不同的观测值具有不同的方差...

所以这让我觉得这可能与它有关。我知道手动执行此操作的理论，但我想知道是否可以使用该lm功能执行此操作。如果没有，是否有任何其他功能（或包）能够做到这一点？

编辑

看到一些评论，这里是一些澄清。举个例子：

x <- 1:10
y <- c(131.4,227.1,245,331.2,386.9,464.9,476.3,512.2,510.8,532.9)
mod <- lm(y ~ x + I(x^2))
summary(mod)

给我：

Residuals:
    Min      1Q  Median      3Q     Max 
-32.536  -8.022   0.087   7.666  26.358 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  39.8050    22.3210   1.783  0.11773    
x            92.0311     9.3222   9.872 2.33e-05 ***
I(x^2)       -4.2625     0.8259  -5.161  0.00131 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 18.98 on 7 degrees of freedom
Multiple R-squared:  0.986, Adjusted R-squared:  0.982 
F-statistic: 246.7 on 2 and 7 DF,  p-value: 3.237e-07

所以基本上，我的系数是 a=39.8±22.3，b=92.0±9.3，c=-4.3±0.8。现在让我们说对于每个数据点，错误是 20。我将weights = rep(20,10)在lm调用中使用，我得到了这个：

Residual standard error: 84.87 on 7 degrees of freedom

但系数上的标准误差不会改变。

手动，我知道如何使用矩阵代数计算协方差矩阵并将权重/误差放在那里，并使用它推导置信区间。那么有没有办法在 lm 函数本身或任何其他函数中做到这一点？

1个回答

这种类型的模型实际上在某些科学分支（例如物理学）和工程学中比“正常”线性回归更常见。因此，在诸如之类的物理工具中ROOT，进行这种类型的拟合是微不足道的，而线性回归并不是本机实现的！物理学家倾向于将其称为“拟合”或卡方最小化拟合。

正态线性回归模型假设每个测量值然后它最大化似然或等效地它的对数因此名称最小二乘——最大化似然是与最小化平方和相同，并且是一个不重要的常数，只要它是常数。对于具有不同已知不确定性的测量，您需要最大化 $\sigma$

L \propto \prod_{i} e^{- \frac{1}{2} {(\frac{y_{i} - (a x_{i} + b)}{σ})}^{2}}

$L \propto \prod_i e^{-\frac{1}{2} \left( \frac{y_i-(ax_i+b)}{\sigma} \right)^2}$

\log (L) = c o n s t a n t - \frac{1}{2 σ^{2}} \sum_{i} (y_{i} - (a x_{i} + b))^{2}

$\log(L) = \mathrm{constant} - \frac{1}{2\sigma^2} \sum_i (y_i-(ax_i+b))^2$

σ

$\sigma$

L \propto \prod e^{- \frac{1}{2} {(\frac{y - (a x + b)}{σ_{i}})}^{2}}

$L \propto \prod e^{-\frac{1}{2} \left( \frac{y-(ax+b)}{\sigma_i} \right)^2}$ 或者等价于它的对数所以，您实际上希望通过逆方差对测量值进行加权，而不是方差。这是有道理的——更准确的测量具有更小的不确定性，应该给予更多的权重。请注意，如果此权重是恒定的，它仍然会从总和中扣除。因此，它不会影响估计值，但会影响标准误差，取自的二阶导数。

\log (L) = c o n s t a n t - \frac{1}{2} \sum {(\frac{y_{i} - (a x_{i} + b)}{σ_{i}})}^{2}

$\log(L) = \mathrm{constant} - \frac{1}{2} \sum \left( \frac{y_i-(ax_i+b)}{\sigma_i} \right)^2$

1 / σ_{i}^{2}

$1/\sigma_i^2$

\log (L)

$\log(L)$

但是，在这里我们遇到了物理学/科学与统计学之间的另一个区别。通常在统计中，您期望两个变量之间可能存在相关性，但很少会是准确的。另一方面，在物理学和其他科学中，如果不是因为讨厌的测量误差（例如，而不是），您通常会期望相关性或关系是准确的。您的问题似乎更多地属于物理/工程案例。因此，对与您的测量结果和权重相关的不确定性的解释与您想要的并不完全相同。它会考虑权重，但它仍然认为有一个整体 $F=ma$ $F=ma+\epsilon$ lm $\sigma^2$ 考虑回归误差，这不是你想要的——你希望你的测量误差是唯一的误差。（的lm解释的最终结果是，只有权重的相对值很重要，这就是为什么您作为测试添加的恒定权重没有效果的原因）。这里的问答有更多细节：

lm 权重和标准误

那里的答案中有几个可能的解决方案。特别是，那里的匿名答案建议使用

vcov(mod)/summary(mod)$sigma^2

基本上，lm根据估计的缩放协方差矩阵，并且您想要撤消此操作。然后，您可以从校正后的协方差矩阵中获得所需的信息。试试这个，但如果可以使用手动线性代数，请尝试仔细检查。请记住，权重应该是逆方差。 $\sigma$

编辑

如果你经常做这种事情，你可能会考虑使用ROOT（这似乎在本地做这件事lm，glm但不做）。下面是一个简短的示例，说明如何在ROOT. 首先，ROOT可以通过 C++ 或 Python 使用，而且下载和安装量很大。您可以使用 Jupiter notebook 在浏览器中尝试，点击此处的链接，选择右侧的“Binder”和左侧的“Python”。

import ROOT
from array import array
import math
x = range(1,11)
xerrs = [0]*10
y = [131.4,227.1,245,331.2,386.9,464.9,476.3,512.2,510.8,532.9]
yerrs = [math.sqrt(i) for i in y]
graph = ROOT.TGraphErrors(len(x),array('d',x),array('d',y),array('d',xerrs),array('d',yerrs))
graph.Fit("pol2","S")
c = ROOT.TCanvas("test","test",800,600)
graph.Draw("AP")
c.Draw()

我已将平方根作为值的不确定性。拟合的输出是 $y$

Welcome to JupyROOT 6.07/03

****************************************
Minimizer is Linear
Chi2                      =       8.2817
NDf                       =            7
p0                        =      46.6629   +/-   16.0838     
p1                        =       88.194   +/-   8.09565     
p2                        =     -3.91398   +/-   0.78028

并产生了一个不错的情节：

ROOT fitter 还可以处理值中的不确定性，这可能需要对. 如果有人知道在 R 中执行此操作的本地方法，我将有兴趣学习它。 $x$ lm

第二次编辑

@Wolfgang 上一个问题的另一个答案给出了一个更好的解决方案：包中的rma工具metafor（我最初将该答案中的文本解释为它没有计算截距，但事实并非如此）。将测量值 y 中的方差简化为 y：

> rma(y~x+I(x^2),y,method="FE")

Fixed-Effects with Moderators Model (k = 10)

Test for Residual Heterogeneity: 
QE(df = 7) = 8.2817, p-val = 0.3084

Test of Moderators (coefficient(s) 2,3): 
QM(df = 2) = 659.4641, p-val < .0001

Model Results:

         estimate       se     zval    pval    ci.lb     ci.ub     
intrcpt   46.6629  16.0838   2.9012  0.0037  15.1393   78.1866   **
x         88.1940   8.0956  10.8940  <.0001  72.3268  104.0612  ***
I(x^2)    -3.9140   0.7803  -5.0161  <.0001  -5.4433   -2.3847  ***

---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

对于我发现的这种类型的回归，这绝对是最好的纯 R 工具。

其它你可能感兴趣的问题

上一篇为什么随机游走不是平稳过程？下一篇什么是交叉验证折叠，或者这个短语没有意义？