增加观测值的数量有 R 抛出随机系数 - 数值稳定性问题？-解网

问：

我有这个代码

rm(list=ls())
N = 20000
xvar <- runif(N, -10, 10) 
e <- rnorm(N, mean=0, sd=1)
yvar <- 1 + 2*xvar + e
plot(xvar,yvar)
lmMod <- lm(yvar~xvar)
print(summary(lmMod))

我预计系数类似于 [1,2]。

取而代之的是，R 不断向我抛出没有统计学意义且不适合模型的随机数，$R^2$ 真的很低。我只是看不出我做错了什么。下面是一个示例输出：N =20000

Call:
lm(formula = yvar ~ xvar)

Residuals:
   Min     1Q Median     3Q    Max 
-47.23  -9.10   1.24  11.23  23.74 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  0.03163    0.08291   0.381  0.70286   
xvar         0.04290    0.01427   3.006  0.00265 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 11.73 on 19998 degrees of freedom
Multiple R-squared:  0.0009635, Adjusted R-squared:  0.0009135 
F-statistic: 19.29 on 1 and 19998 DF,  p-value: 1.131e-05

但是，如果我输入 N=200 或 N=2000，它就可以工作。系数与实数相似，在实数的两个标准差内，我得到的 $R^2$ 值高达 99%，并且系数在统计学上都是显著的，为 $p<<0.01$。

这是怎么回事？为什么增加观测值的数量会使回归恶化？R 是否暗中遇到数值稳定性问题？

我在 Kubuntu 19.04 上运行 R 3.6.0。使用 --vanilla 选项在命令行上运行 R 也会出现同样的问题。

编辑：这是输出sessioninfo()

> sessionInfo()
R version 3.6.0 (2019-04-26)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 19.04

Matrix products: default
BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/libmkl_rt.so

Random number generation:
 RNG:     Mersenne-Twister 
 Normal:  Inversion 
 Sample:  Rounding 

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=it_IT.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=it_IT.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=it_IT.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=it_IT.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
[1] compiler_3.6.0 tools_3.6.0

r 线性回归浮点精度数值稳定性

library(flexiblas)

mkl <- "$HOME/intel/oneapi/mkl/latest/lib/intel64/libmkl_rt.so"
idx <- flexiblas_load_backend(mkl)
flexiblas_switch(idx)
sessionInfo()
print(flexiblas_current_backend())

rm(list=ls())
N = 100000
xvar <- runif(N, -10, 10) 
e <- rnorm(N, mean=0, sd=1)
yvar <- 1 + 2*xvar + e
plot(xvar,yvar)
lmMod <- lm(yvar~xvar)
print(summary(lmMod))

返回正确的系数。

增加观测值的数量有 R 抛出随机系数 - 数值稳定性问题？

Increasing the number of observations have R throw random coefficients - Numerical stability problem?

评论

评论

评论