• Archive by category "R语言实践"

Blog Archives

线性和非线性的最小二乘回归

算法为王系列文章,涵盖了计算机算法,数据挖掘(机器学习)算法,统计算法,金融算法等的多种跨学科算法组合。在大数据时代的背景下,算法已经成为了金字塔顶的明星。一个好的算法可以创造一个伟大帝国,就像Google。

算法为王的时代正式到来….

环球彩票关于 作者:

  • 张丹(Conan), 程序员/Quant: Java,R,Nodejs
  • blog: http://xiyueit.com
  • email: bsspirit@gmail.com

转载请注明出处:
http://xiyueit.com/r-lsm-regression

前言

很多经典的机器学习模型都是在最小二乘法的基础上发展起来的,理解最小二乘法的原理,对于环球彩票环球彩票我 们 掌握机器学习的高级算法是非常重要的,尤其是回归模型中。最小二乘法是主要是用来做函数拟合,或者求函数极值等问题。

本文以回归模型为例,介绍在最小二乘法在线性和非线性模型中的应用。

目录

  1. 最小二乘法介绍
  2. 最小二乘求解
  3. 线性最小二乘回归
  4. 非线性最小二乘回归

1. 最小二乘法介绍

最小二乘法(Least Squares Method,简记为LSE),源于天文学和测地学上的应用需要,是勒让德( A. M. Legendre)于1805年在其著作《计算慧星轨道的新环球彩票方法 》中提出的,它的主要思想是通过计算理论值与观测值之差的平方和最小,求解未知的参数。

最小二乘法的公式为:

公式解读:

  • H, 目标值,观测值与理论值之差的平方和
  • H’, 最小化的H值
  • y, 观测值
  • yi, 理论值
  • argmin, 最小化函数

环球彩票环球彩票我 们 先了把术语进行统一定义,观测值通常说的是环球彩票环球彩票我 们 的样本数据,理论值就是假设拟合函数产生的预测数据,目标函数是判断一个模型的好不好的标准,进行最小二乘法建立目标函数,让目标函数最小化。

以上图为线性拟合为例,观测值为y1,y2,…,y7,理论值为y1′,y2′,….,y7’,是线上的点并与观测值的x相同,红色线就是环球彩票环球彩票我 们 要做的拟合直线,环球彩票环球彩票我 们 可以画出无数条这样的直线,判断标准是让这条线上的点的理论值与观测值的误差平方和最小,即最小二乘法。

SSE=sum((y1-y1')^2+(y2-y2')^2+(y3-y3')^2+(y4-y4')^2+(y5-y5')^2+(y6-y6')^2+(y7-y7')^2)

利用最小二乘法,可以对数据进行线性和非线性拟合,使误差平方和(SSE)或残差平方和最小。如果观测到的误差近似正态分布时,这种环球彩票方法 是非常有效的。环球彩票环球彩票我 们 在线性回归模型中,要进行残差的正态分布检验,就是基于这个目的的,可以参考文章R语言解读一元线性回归模型

2. 最小二乘求解

在线性方程的求解或是数据曲线拟合中,利用最小二乘法求得的解则被称为最小二乘解。最小二乘法(又称最小平环球彩票方法 )是一种数学环球彩票优化 环球彩票技术 ,它通过最小化误差的平方和寻找数据的最佳函数匹配。其他一些环球彩票优化 问题也可通过最小化能量或最大化熵用最小二乘法来表达。

定义数据集x和y,两个变量。

> x<-c(6.19,2.51,7.29,7.01,5.7,2.66,3.98,2.5,9.1,4.2)
> y<-c(5.25,2.83,6.41,6.71,5.1,4.23,5.05,1.98,10.5,6.3)

可视化,画出x,y变量的散点图。

> plot(x,y)

对数据进行最小二乘求解,找到最小二乘解。


> fit<-lsfit(x,y);fit
$coefficients
Intercept         X
0.8310557 0.9004584

$residuals
[1] -1.1548933 -0.2612063 -0.9853975 -0.4332692 -0.8636686  1.0037250  0.6351198 -1.1022017
[9]  1.4747728  1.6870190

$intercept
[1] TRUE

$qr
$qt
[1] -17.1901414   6.2044421  -0.7047339  -0.1530724  -0.5856560   1.2766692   0.9102648
[8]  -0.8295242   1.7584540   1.9625308

$qr
Intercept           X
[1,] -3.1622777 -16.1718880
[2,]  0.3162278   6.8903149
[3,]  0.3162278  -0.2782874
[4,]  0.3162278  -0.2376506
[5,]  0.3162278  -0.0475287
[6,]  0.3162278   0.3936703
[7,]  0.3162278   0.2020970
[8,]  0.3162278   0.4168913
[9,]  0.3162278  -0.5409749
[10,]  0.3162278   0.1701682

$qraux
[1] 1.316228 1.415440

$rank
[1] 2

$pivot
[1] 1 2

$tol
[1] 1e-07

attr(,"class")
[1] "qr"

结果解读:

  • coefficients, 系数的最小二乘估计,包括截距Intercept和自变量X
  • residuals, 残差
  • intercept, 有截距
  • qt, 矩阵QR分解。QR分解法是将矩阵分解成一个正规正交矩阵与上三角形矩阵。环球彩票关于 矩阵的详细操作,请参考文章R语言中的矩阵计算

环球彩票环球彩票我 们 看到lsfit()函数使用矩阵的QR分解的环球彩票方法 ,计算得出了最小二乘估计的截距Intercept和自变量系数。

目标:y=a*x+b,用最小二乘法估计出参数a和b。

根据最小二乘法的矩阵计算公式,环球彩票我 直接跳过了推到过程,得出结论。

公式解释:

  • X,为自变量矩阵
  • Y,为因变量矩阵
  • B,为参数,需要求解
  • XT,X的转置矩阵
  • *, 矩阵乘法
  • -1, 计算逆矩阵

手动计算过程:


# 因变量矩阵
> y1<-matrix(y,ncol=1)  

# 自变量矩阵,增加截距列
> x1<-matrix(c(x,rep(1,length(x))),ncol=2)

# 求解最小二乘解
> solve(t(x1)%*%x1) %*% t(x1) %*% y1
          [,1]
[1,] 0.9004584
[2,] 0.8310557

对应到一元线性回归方程中,y = 0.9004584 * x + 0.8310557,与lsfit()的结果一致。

3. 线性最小二乘回归

接下来,环球彩票环球彩票我 们 先进行线性回归分析,以上面的数据集作为样本,只包含x和y。以普通最小二乘法进行线性回归,计算出预测值并和观测值进行比较,计算预测值和观测值的误差平方和,使得误差函数最小。线性回归模型的详细介绍,可以参考文章R语言解读一元线性回归模型, R语言解读多元线性回归模型

建立线性回归模型


# 建立回归模型
> line<-lm(y~x)

# 查看回归模型
> summary(line)
Call:
lm(formula = y ~ x)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.1549 -0.9550 -0.3472  0.9116  1.6870 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.8311     0.9440   0.880 0.404338    
x             0.9005     0.1698   5.302 0.000726 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.17 on 8 degrees of freedom
Multiple R-squared:  0.7785,	Adjusted R-squared:  0.7508 
F-statistic: 28.12 on 1 and 8 DF,  p-value: 0.0007263

# 残差平方和
> sum(line$residuals^2)
[1] 10.95334

通过查看模型的结果数据,环球彩票环球彩票我 们 可以发现通过T检验自变量x是非常显著,但截距Intercept不显著。通过F检验判断出整个模型的自变量方差是非常显著。通过R^2的相关系数检验可以判断自变量和因变量是不是高度相关的,残差平方和(residual sum-of-squares)为10.95334,效果不好。

接下来,使用模型进行预测,生成预测值。


# 设置x的取值区间
> new <- data.frame(x = seq(min(x),max(x),len = 100))

# 预测结果
> pred<-predict(line,newdata=new)

进行可视化,把观测值与预测值进行展示,以散点表现观测值,以线表示预测值。


> plot(x,y)
> lines(new$x,pred,col="red",lwd=2)

从图中环球彩票环球彩票我 们 可以很容易看出,拟合的效果并不好,与R^2检验不通过是一致的。那么,如果线性模型表现不好,环球彩票环球彩票我 们 可以试试非线性的模型,还是用最小二乘法进行非线性的回归。

4. 非线性最小二乘回归

由于散点的分布并不是线性的,环球彩票环球彩票我 们 再尝试用非线性的函数,进行最小二乘回归的拟合,包括一元二次函数,一元三次函数,指数函数,倒数函数。在R语言中,环球彩票环球彩票我 们 可以使用nls()函数,进行非线性最小二乘回归。

4.1 一元二次函数非线性拟合
定义公式 y ~ a*x+b*x^2+c ,其中x, x^2是一次函数和二次函数,a,b,c分别是参数,环球彩票环球彩票我 们 需要给定参数的初始值。


# 建立回归模型
> nline1 <- nls(y ~ a*x+b*x^2+c, start = list(a = 1,b = 1,c=2))

# 查看模型结果
> nline1
Nonlinear regression model
  model: y ~ a * x + b * x^2 + c
   data: parent.frame()
       a        b        c 
0.008987 0.082188 2.850389 
 residual sum-of-squares: 9.758

Number of iterations to convergence: 1 
Achieved convergence tolerance: 3.024e-06

# 分析模型
> summary(nline1)

Formula: y ~ a * x + b * x^2 + c

Parameters:
  Estimate Std. Error t value Pr(>|t|)
a 0.008987   0.977890   0.009    0.993
b 0.082188   0.088760   0.926    0.385
c 2.850389   2.379763   1.198    0.270

Residual standard error: 1.181 on 7 degrees of freedom

Number of iterations to convergence: 1 
Achieved convergence tolerance: 3.024e-06

通过查看模型的结果数据,环球彩票环球彩票我 们 可以发现通过T检验相关系数a,b,c都不显著,残差平方和(residual sum-of-squares)为9.758,数值过大,效果不好。

把观测值与预测值进行可视化展示,以散点表现观测值,以线表示预测值。


# 模型预测
> npred1<-predict(nline1,newdata=new)

# 可视化
> plot(x,y)
> lines(new$x,npred1,col="blue",lwd=2)

从图中环球彩票环球彩票我 们 可以看出,拟合的是一条曲线,效果不太好。接下来,环球彩票环球彩票我 们 再试试一元三次函数。

4.2 一元三次函数非线性拟合

定义公式 y ~ a*x+b*x^2+c*x^3+d ,其中x, x^2, x^3是一次函数,二次函数和三次函数,a,b,c,d分别是参数,环球彩票环球彩票我 们 需要给定参数的初始值。

使用一元三次函数,进行非线性拟合。


> nline2 <- nls(y ~ a*x+b*x^2+c*x^3+d, start = list(a = 1,b = 1,c=2,d=1))

# 查看模型结果
> nline2
Nonlinear regression model
  model: y ~ a * x + b * x^2 + c * x^3 + d
   data: parent.frame()
      a       b       c       d 
 10.011  -1.822   0.110 -12.516 
 residual sum-of-squares: 3.453

Number of iterations to convergence: 2 
Achieved convergence tolerance: 3.983e-08

# 分析模型
> summary(nline2)

Formula: y ~ a * x + b * x^2 + c * x^3 + d

Parameters:
   Estimate Std. Error t value Pr(>|t|)  
a  10.01051    3.08630   3.244   0.0176 *
b  -1.82152    0.57797  -3.152   0.0198 *
c   0.10998    0.03323   3.310   0.0162 *
d -12.51586    4.88778  -2.561   0.0429 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7586 on 6 degrees of freedom

Number of iterations to convergence: 2 
Achieved convergence tolerance: 3.983e-08

通过查看模型的结果数据,环球彩票环球彩票我 们 可以发现通过T检验相关系数a,b,c,d在0.01的置信区间显著,残差平方和(residual sum-of-squares)为3.453,数值较小,效果应该算是还可以的,能够表达拟合的特征。

把观测值与预测值进行可视化展示,以散点表现观测值,以线表示预测值。


> npred2<-predict(nline2,newdata=new)
> plot(x,y)
> lines(new$x,npred2,col="green",lwd=2)

从图中环球彩票环球彩票我 们 可以看出,拟合的曲线基本穿过观测值的区域,直观上说明效果还不错,与数据检验给出的结论是一致的。

4.3 指数函数非线性拟合

定义公式 y ~ a*exp(b*x)+c ,其中exp(x)是指数函数,a,b,c分别是参数,环球彩票环球彩票我 们 需要给定参数的初始值。


> nline3 <- nls(y ~ a*exp(b*x)+c, start = list(a = 1,b = 1,c=2))

> nline3
Nonlinear regression model
  model: y ~ a * exp(b * x) + c
   data: parent.frame()
     a      b      c 
0.6720 0.2718 2.2087 
 residual sum-of-squares: 8.966

Number of iterations to convergence: 17 
Achieved convergence tolerance: 7.136e-06

> summary(nline3)

Formula: y ~ a * exp(b * x) + c

Parameters:
  Estimate Std. Error t value Pr(>|t|)
a   0.6720     1.2127   0.554    0.597
b   0.2718     0.1764   1.541    0.167
c   2.2087     2.2447   0.984    0.358

Residual standard error: 1.132 on 7 degrees of freedom

Number of iterations to convergence: 17 
Achieved convergence tolerance: 7.136e-06

通过查看模型的结果数据,环球彩票环球彩票我 们 可以发现通过T检验相关系数a,b,c都不显著,残差平方和(residual sum-of-squares)为8.966,数值过大,效果不好。


> npred3<-predict(nline3,newdata=new)
> plot(x,y)
> lines(new$x,npred3,col="yellow",lwd=2)

从图中环球彩票环球彩票我 们 可以看出,指数函数的拟合曲线,和二次函数表现差不多,效果不太好。

4.4 倒数函数非线性拟合

定义公式 y ~ y ~ b/x+c ,其中1/x是倒数函数,b,c分别是参数,环球彩票环球彩票我 们 需要给定参数的初始值。

使用倒数函数,进行拟合。

> nline4 <- nls(y ~ b/x+c, start = list(b = 1,c=2))

# 查看模型结果
> nline4
Nonlinear regression model
  model: y ~ b/x + c
   data: parent.frame()
    b     c 
-17.0   9.5 
 residual sum-of-squares: 15.74

Number of iterations to convergence: 1 
Achieved convergence tolerance: 1.259e-07

> summary(nline4)

Formula: y ~ b/x + c

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
b  -17.003      4.108  -4.139  0.00326 ** 
c    9.500      1.077   8.817 2.15e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.403 on 8 degrees of freedom

Number of iterations to convergence: 1 
Achieved convergence tolerance: 1.259e-07

通过查看模型的结果数据,环球彩票环球彩票我 们 可以发现通过T检验相关系数b,c非常显著,残差平方和(residual sum-of-squares)为15.74,数值较大,效果不好。

把观测值与预测值进行可视化展示,以散点表现观测值,以线表示预测值。


> npred4<-predict(nline4,newdata=new)
> plot(x,y)
> lines(new$x,npred4,col="gray",lwd=2)


从图中环球彩票环球彩票我 们 可以看出,倒数函数的拟合曲线,确实也效果不太好,而且最后一个点有个想离群点,如果没有最后一个点这个效果应该是更好的。

4.5 模型环球彩票优化
根据上面的几种模型验证的结果,环球彩票我 手动进行了环球彩票优化 ,发现高次函数进行拟合效果会更好。定义公式 y ~ a*x^3+b*x^5+c*x^7+e*d^9+e*x^11 ,其中高次函数定义为11次,并去掉了截距。


> nline5 <- nls(y ~ a*x^3+b*x^5+c*x^7+d*x^9+e*x^11,
+               start = list(a=1,b=2,c=1,d=1,e=1))

# 查看模型结果
> nline5
Nonlinear regression model
  model: y ~ a * x^3 + b * x^5 + c * x^7 + d * x^9 + e * x^11
   data: parent.frame()
         a          b          c          d          e 
 2.959e-01 -2.064e-02  5.837e-04 -7.324e-06  3.367e-08 
 residual sum-of-squares: 2.583

Number of iterations to convergence: 3 
Achieved convergence tolerance: 2.343e-07

# 模型分析
> summary(nline5)

Formula: y ~ a * x^3 + b * x^5 + c * x^7 + d * x^9 + e * x^11

Parameters:
    Estimate Std. Error t value Pr(>|t|)   
a  2.959e-01  5.132e-02   5.767   0.0022 **
b -2.064e-02  6.176e-03  -3.341   0.0205 * 
c  5.837e-04  2.468e-04   2.365   0.0644 . 
d -7.324e-06  3.940e-06  -1.859   0.1222   
e  3.367e-08  2.139e-08   1.574   0.1763   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7187 on 5 degrees of freedom

Number of iterations to convergence: 3 
Achieved convergence tolerance: 2.343e-07

通过查看模型的结果数据,环球彩票环球彩票我 们 可以发现通过T检验相关系数a,b,c是显著,残差的误差平方和(residual sum-of-squares)为2.583,数值是最小的,而且比三次函数增长了2上系数,效果应该算不错的,能够表达拟合的特征。

把观测值与预测值进行可视化展示,以散点表现观测值,以线表示预测值。


> npred5<-predict(nline5,newdata=new)
> plot(x,y)
> lines(new$x,npred5,col="purple3",lwd=2)

从图中环球彩票环球彩票我 们 可以看出,拟合的曲线穿过观测值的区域,并且非常贴近部分数据点,直观上说明效果很好,与数据检验给出的结论是一致的。

4.6 总结

环球彩票环球彩票我 们 把上面的1种线性函数和5种非线性的拟合函数的模型评价指标,整理到一个表格中。其中,三次函数和高次函数的拟合效果是不错的,而其他的几种拟合效果,都是不太好的。判断好坏的标准,就是残差的误差平方和,同时参数的T检验显著,才是比较好的拟合模型。

拟合函数 残差平方和 参数T检验显著性 解释性 颜色
一次函数 10.95334 不显著 红色
二次函数 9.758 不显著 蓝色
三次函数 3.453 显著 不好 绿色
指数函数 8.966 不显著 不好 黄色
倒数函数 15.74 显著 灰色
高次函数 2.583 显著 无法解释 紫色

最后,合并所有拟合曲线,画到同一张图上,进行对比。绿色的三次函数拟合曲线,和紫色的高次函数拟合曲线效果最好。

但是,高次函数也带来了的负面问题,一是过拟合,二是业务的可解释性非常差。如果x和y分别带入到业务场景,比如x是货物的重量和y是货物的单价,那么x的11次方,完全是无法理解的,等同于一个黑盒的模型。所以,真正落地项目在使用的时候,要根据业务场景的要求,在准确度,可解释性,复杂度等标准之下,进行权衡。

本文利用最小二乘法,用R语言进行了线性回归和非线性回归的实验。验证了最小二乘法在处理回归问题上是非常有效的,除了回归问题,最小二乘法还能解决分类的问题。最小二乘法作为机器学习的底层算法,是非常值得环球彩票环球彩票我 们 学习和掌握的,明白原理就能掌握一类的高级模型。

转载请注明出处:
http://xiyueit.com/r-lsm-regression

打赏作者

用R语言填充缺失值mice

R的极客理想系列文章,涵盖了R的思想,使用,环球彩票工具 ,创新等的一系列要点,以环球彩票我 个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,环球彩票互联网 ….都在使用R语言。

要成为有理想的极客,环球彩票环球彩票我 们 不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让环球彩票环球彩票我 们 一起动起来吧,开始R的极客理想。

环球彩票关于 作者:

  • 张丹(Conan), 程序员/Quant: Java,R,Nodejs
  • blog: http://xiyueit.com
  • email: bsspirit@gmail.com

转载请注明出处:
http://xiyueit.com/r-na-mice

前言

在进行数据分析时,从真实场景中所获得的数据经常是有缺失的,缺失数据会给环球彩票环球彩票我 们 数据分析造成很大的麻烦。很多的时候,环球彩票环球彩票我 们 需要先对缺失数据进行填充。数据缺失的原因是各种各样的,所以填充环球彩票方法 也有很多种,要找到符合业务特点的填充环球彩票方法 ,不能因为填充而导致原数据集特征的巨变。

本文通过mice包提供的环球彩票方法 ,合理地估算缺失值的分布,再进行填充,确保填充的数据是合理的。

目录

  1. 缺失的数据
  2. 准备数据集
  3. mice包介绍
  4. 填充环球彩票方法 和程序实现

1. 缺失的数据

缺失值的处理,是数据预处理过程中的重要环节,造成数据缺失的原因有很多,环球彩票我 发现的最主要的原因是,业务允许为空值导致数据集有缺失值,还包括比如,数据丢失、存储故障和拒绝填写信息等原因。

如果数据量很大而缺失值量很少,那么环球彩票环球彩票我 们 直接去除这些数据就行了,也不会有太多的影响。缺失值处理主要讨论的是,数据量不太大,而缺失值的数据的影响又比较关键时,这个时候就需要用到一些环球彩票方法 对缺失值数据进行填充。

1.1 缺失机制分类

Little和Ruth(1987)把数据缺失的机制分为三类:

  1. 完全随机缺失(missing completely at random, MCAR)
    所缺失的数据是完全随机的,缺失发生的概率既与已观察到的数据无关,也与未观察到的数据无关。这是一种比较理想的情况。
  2. 随机缺失(missing at random, MAR)
    数据的缺失不是完全随机的。缺失数据发生的概率与所观察到的变量是有关的,而与未观察到的数据的特征是无关的。这是一个比较严重的问题,在这种情况下,环球彩票环球彩票我 们 需要进一步检查数据收集过程,并尝试了解数据为什么丢失。 例如,如果在一项问卷调查中,大多数人没有回答某个问题,他们为什么这么做,是问题不清楚吗?
  3. 不可忽略的缺失(non-ignorable missing ,NIM),亦称为非随机缺失(not missing at random, NMAR),也有研究者将其称为MNAR(missing not at random)。
    缺失数据不仅依赖于其它变量,又依赖于变量本身,这种缺失即为不可忽略的缺失。

假设数据缺失是完全随机缺失(MCAR),如果数据缺失太多,也可能造成很大的问题。 一般来说,数据的缺失量小于数据总量的5%是可以接受的。如果某个特征或样本的缺失数据量超过20%,那么您可考虑是否需要留下该特征或样本。

1.2 填充环球彩票方法

对于缺失值的填充环球彩票方法 分为单变量填充和多变量填充。单变量填补,包括简单随机填充,均值/中位数/分位数填充,回归填充,热平台和冷平台。多变量填补,使用回归插补法。

下面分别介绍一下,这种填充的环球彩票方法 :

简单随机填补:对于每一个缺失值,从已有的该变量数据中随机抽样作为填补值,填补进缺失位置。仅仅考虑到了缺失变量本身,而并没有考虑到相关变量的信息。因此,信息量的利用少。

均值/中位数/分位数填补:用存在缺失值的变量的已有值的均值/中位数/分位数,作为填补值。这种环球彩票方法 显然会导致方差偏小。

回归填补:将缺失变量作为因变量,相关变量(其他变量)作为自变量,进行回归拟合,用预测值作为填补值。用于作为自变量的变量最好是具有完全数据(无缺失)。

热平台和冷平台:热平台法又称匹配插补法,思路是在完全数据样本中,找到一个和具有缺失值的样本相似的完全数据样本,用完全数据样本值作为填充值,其过程有点类似于K阶近邻的思想。冷平台法又称条件均值插补法,思路是先将总体分层(聚类),采用样本所在层(类)的完全数据的均值来替代缺失值。

回归插补法:对缺失变量和完全数据变量拟合多元回归模型来预测缺失值,是多重填补法的一种应用(Multiple Imputation Missing Data)。本文主要就介绍这种环球彩票方法 ,利用mice包进行实现。

2. 准备数据集

为了进行空值的数据处理,环球彩票环球彩票我 们 使用一个R语言中自带的数据集airquality。airquality数据集,记录了1973年5月至9月在纽约进行的每日空气质量测量。

开发所使用的系统环境

  • 环球彩票Win 10 64bit
  • R: 3.6.1 x86_64-w64-mingw32/x64 b4bit

数据集包括153条记录,6列:

  • Ozone: 臭氧浓度(ppb)
  • Solar.R: 太阳辐射量(lang)
  • 环球彩票Win d: 风速(MPH)
  • Temp: 温度(°F)
  • Month: 发生月
  • Day:发生日

# 查看数据集
> head(airquality)
  Ozone Solar.R 环球彩票Win
d Temp Month Day
1    41     190  7.4   67     5   1
2    36     118  8.0   72     5   2
3    12     149 12.6   74     5   3
4    18     313 11.5   62     5   4
5    NA      NA 14.3   56     5   5
6    28      NA 14.9   66     5   6

查看数据统计概率, Ozone列有37条空值,Solar.R列有7条空值,环球彩票Win d列有7条空值。


> summary(airquality)
     Ozone           Solar.R           环球彩票Win
d             Temp           Month            Day      
 Min.   :  1.00   Min.   :  7.0   Min.   : 1.700   Min.   :56.00   Min.   :5.000   Min.   : 1.0  
 1st Qu.: 18.00   1st Qu.:115.8   1st Qu.: 7.400   1st Qu.:72.00   1st Qu.:6.000   1st Qu.: 8.0  
 Median : 31.50   Median :205.0   Median : 9.700   Median :79.00   Median :7.000   Median :16.0  
 Mean   : 42.13   Mean   :185.9   Mean   : 9.958   Mean   :77.88   Mean   :6.993   Mean   :15.8  
 3rd Qu.: 63.25   3rd Qu.:258.8   3rd Qu.:11.500   3rd Qu.:85.00   3rd Qu.:8.000   3rd Qu.:23.0  
 Max.   :168.00   Max.   :334.0   Max.   :20.700   Max.   :97.00   Max.   :9.000   Max.   :31.0  
 NA's   :37       NA's   :7             

计算缺失值比例


# 定义缺失值比例计算函数miss()
> miss <- function(x){sum(is.na(x))/length(x)}

列,缺失值比例


> apply(airquality,2,miss)
     Ozone    Solar.R       环球彩票Win
d       Temp      Month        Day 
0.24183007 0.04575163 0.00000000 0.00000000 0.00000000 0.00000000 

Ozone列的缺少值比例最高,达到24%,环球彩票环球彩票我 们 可以考虑移除这个变量。

行,缺失值比例。


# 计算行的缺失值比例,取前30条记录
> head(apply(airquality,1,miss),30)
  [1] 0.0000000 0.0000000 0.0000000 0.0000000 0.3333333 0.1666667 0.0000000 0.0000000 0.0000000 0.1666667
 [11] 0.1666667 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
 [21] 0.0000000 0.0000000 0.0000000 0.0000000 0.1666667 0.1666667 0.3333333 0.0000000 0.0000000 0.0000000

其中,第5个值缺失值表示第5行缺失比例为33%,表示前6列中有2列是空的。

3. mice包介绍

mice包实现了,多变量缺失数据填补的环球彩票方法 ,基于多元回归模型来预测,每个不完整变量由单独的模型估算,进行多重混合评估。mice包算法可以实现对于连续型,二进制,无序分类和有序分类数据的进行混合。此外,mice包可以处理连续的两级数据,并通过被动插补来保持插补之间的一致性,通过第各种统计诊断图,保证插补数据的质量。

按照mice包把缺失值的处理形成了一套流程,通过链式方程形成多重填补法。

对于缺失值数据的处理,用3个步骤来进行定义。

  1. 填充:mice()函数,从一个包含缺失数据的数据框开始,然后返回一个包含多个完整数据集的对象,每个完整数据集都是通过对原始数据框中的缺失数据进行插而生成的。
  2. 分析:with()函数,可依次对每个完整数据集应用统计模型,分析填充的结果。
  3. 环球彩票优化 :pool()函数,将这些单独的分析结果整合为一组结果,最终模型的标准误和p值,都将准确地反映出由于缺失值和多重插补而产生的不确定性。

4. 填充环球彩票方法 和程序实现

接下来,环球彩票环球彩票我 们 使用mice包,对airquality数据集的缺失值进行填充。操作步骤按,观察缺失值,多重填补法,分析结果环球彩票优化 模型,获得最终结果。

4.1 观察缺失值

对于上面的数据集,环球彩票环球彩票我 们 可以用mice包来观察数据。用md.pattern()函数来判断行的缺失值。


> md.pattern(airquality)
    环球彩票Win
d Temp Month Day Solar.R Ozone   
111    1    1     1   1       1     1  0
35     1    1     1   1       1     0  1
5      1    1     1   1       0     1  1
2      1    1     1   1       0     0  2
       0    0     0   0       7    37 44

数据解读:

  • 左边第一列,为样本数
  • 右边第一列,为累积缺失值个数,0为没有缺失。
  • 最下面一行,每个列有多少个缺失值

举例解释一下:第一行可以理解为,有111条记录,所有列都有值,没有缺失值。第二行可以理解为,有35条记录,Ozone列有缺失值。

生成缺失值矩阵图片。

接下来,以行和列进行配对处理,按4种模式分别展示出缺失值。

  • rr:行有值,列有值
  • rm:行有值,列丢失
  • mr:行丢失,列有值
  • mm:行丢失,列丢失

> md.pairs(airquality)
$rr
        Ozone Solar.R 环球彩票Win
d Temp Month Day
Ozone     116     111  116  116   116 116
Solar.R   111     146  146  146   146 146
环球彩票Win
d      116     146  153  153   153 153
Temp      116     146  153  153   153 153
Month     116     146  153  153   153 153
Day       116     146  153  153   153 153

$rm
        Ozone Solar.R 环球彩票Win
d Temp Month Day
Ozone       0       5    0    0     0   0
Solar.R    35       0    0    0     0   0
环球彩票Win
d       37       7    0    0     0   0
Temp       37       7    0    0     0   0
Month      37       7    0    0     0   0
Day        37       7    0    0     0   0

$mr
        Ozone Solar.R 环球彩票Win
d Temp Month Day
Ozone       0      35   37   37    37  37
Solar.R     5       0    7    7     7   7
环球彩票Win
d        0       0    0    0     0   0
Temp        0       0    0    0     0   0
Month       0       0    0    0     0   0
Day         0       0    0    0     0   0

$mm
        Ozone Solar.R 环球彩票Win
d Temp Month Day
Ozone      37       2    0    0     0   0
Solar.R     2       7    0    0     0   0
环球彩票Win
d        0       0    0    0     0   0
Temp        0       0    0    0     0   0
Month       0       0    0    0     0   0
Day         0       0    0    0     0   0

为了让可视化效果更强一点,环球彩票我 可以引入VIM包,进行可视化。


# 加载VIM包
> library(VIM)
> mice_plot <- aggr(airquality, col=c('navyblue','yellow'),
+                     numbers=TRUE, sortVars=TRUE,
+                     labels=names(airquality), cex.axis=.7,
+                     gap=3, ylab=c("Missing data","Pattern"))

 Variables sorted by number of missings: 
 Variable      Count
    Ozone 0.24183007
  Solar.R 0.04575163
     环球彩票Win
d 0.00000000
     Temp 0.00000000
    Month 0.00000000
      Day 0.00000000

左侧图显示了每个字段中缺失样本数量占比。右侧图每一行代表了一种缺失模式,黄色代表缺失,蓝色代表未缺失,右侧图例表示此模式数量,可与md.pattern()结果对应观察。

画出数据集 airquality 的各列的数据分布。

> matrixplot(airquality)


图中,红色表示缺失值,从黑色到灰色再到白色的过渡色表示值的不同,过渡色越多表示数据取值越分散,过渡色越少表示数据取值越集中。Month列,只有5个颜色,表示只有5个取值。

环球彩票环球彩票我 们 还可以用散点图,来观察两个变量之间的缺失值的关系。画出Ozone和Solar.R列的空缺值的分布散点图。


> marginplot(airquality[c(1,2)], col = c("blue", "red", "orange"))

x轴为Ozone列,y轴为Solar.R列,蓝色点为数据散点图,红色点分别表示Ozone的缺失值和Solar.R列的缺失值位置。

4.2 多重填补法

接下来,使用mice()函数进行数据填充,


# 进行填充
> imputed_Data <- mice(airquality, m=5, maxit = 50, method = 'pmm', seed = 500)
 iter imp variable
  1   1  Ozone  Solar.R
  1   2  Ozone  Solar.R
  1   3  Ozone  Solar.R
  1   4  Ozone  Solar.R
  1   5  Ozone  Solar.R
  2   1  Ozone  Solar.R
  2   2  Ozone  Solar.R
  2   3  Ozone  Solar.R
  2   4  Ozone  Solar.R
  2   5  Ozone  Solar.R

// 省略输出

参数解读:

  • airquality, 原始数据集
  • m, 多重填补法的填补矩阵数,默认为5次
  • maxit,最大迭代次数,默认5次
  • method, 填补用的环球彩票方法 ,pmm为预测均值匹配(predictive mean matching)。用methods(mice) 可以看到有哪些可用的环球彩票方法 。

查看统计结果,输出填充算法,预测矩阵predictorMatrix。填充算法这里使用的是pmm, 除了pmm还有很多种。预测矩阵中,显示为1的位置是会进行填充值的位置,如果不想进行填充,可以修改矩阵中的值为0.


> summary(imputed_Data)
Class: mids
Number of multiple imputations:  5 
Imputation methods:
  Ozone Solar.R    环球彩票Win
d    Temp   Month     Day 
  "pmm"   "pmm"      ""      ""      ""      "" 
PredictorMatrix:
        Ozone Solar.R 环球彩票Win
d Temp Month Day
Ozone       0       1    1    1     1   1
Solar.R     1       0    1    1     1   1
环球彩票Win
d        1       1    0    1     1   1
Temp        1       1    1    0     1   1
Month       1       1    1    1     0   1
Day         1       1    1    1     1   0

查看填补的结果


> imputed_Data$imp
$Ozone
     1   2   3  4   5
5   18  18  18  6   1
10  20  30  16 16  16
25  18  19   8  6  13
26   4   1  28 19  37
27  14   7  11 18  41
32  59  40  44 23  31
33  36   7  13 65  41
// 省略输出

$Solar.R
     1   2   3   4   5
5  215 299  65  98 112
6  222 157 320 332 279
11 183  47 137 135   7
27  36 223  25 238 238
96 273 323 291 294 175
97 287 275 285  31 237
98 323 187 294 320 264

// 省略输出

以$Ozone为例,左边第一列为缺失值的行号,其他每一列都为每次pmm算法生成的填充结果,每一行表示生成的5次结果分别是什么值。

分面板观察,以独立的各个指标进行分组,5组数据的填充情况。


> stripplot(imputed_Data, col=c("grey",mdc(2)),pch=c(1,20))

每一个面板为一个列的原始数据和5次填充的数据,灰色散点为原始数据取值的分布,红色的小点为填充的数据取值的分布。

分面板观察,以Ozone和Solar.R做进行分组组合,观察5组数据的填充情况。


> xyplot(imputed_Data , Ozone ~  Solar.R | .imp, pch=20,cex=1.2)

每一个面板为一个列的原始数据和5次填充的数据,x轴为Solar.R,y轴为Ozone,蓝色散点为原始数据取值的分布,红色的小点为填充的数据取值的分布。

3.3 分析结果,环球彩票优化 模型
使用生成的填补的结果,应用统计模型,对每个完整的数据集进行评价。这里环球彩票环球彩票我 们 以多元线性回归模型作为评价填充值好坏的模型。

使用with()函数,对5组插补数据集进行多元线性回归分析模型,进行T检验,判断数据集中每个变量的有效性。


> fit=with(imputed_Data,lm(Ozone ~ 环球彩票Win
d + Solar.R + Temp))

# 查看统计结果
> summary(fit)
# A tibble: 20 x 5
   term        estimate std.error statistic  p.value
   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
 1 (Intercept) -39.9      19.1        -2.09 3.87e- 2
 2 环球彩票Win
d         -3.60      0.549      -6.55 8.90e-10
 3 Solar.R       0.0572    0.0200      2.87 4.72e- 3
 4 Temp          1.38      0.213       6.46 1.39e- 9
 5 (Intercept) -69.2      19.7        -3.51 5.99e- 4
 6 环球彩票Win
d         -3.24      0.567      -5.71 5.81e- 8
 7 Solar.R       0.0582    0.0204      2.86 4.84e- 3
 8 Temp          1.71      0.218       7.85 7.50e-13
 9 (Intercept) -59.5      19.2        -3.10 2.34e- 3
10 环球彩票Win
d         -3.13      0.550      -5.68 6.85e- 8
11 Solar.R       0.0571    0.0198      2.88 4.54e- 3
12 Temp          1.57      0.215       7.29 1.66e-11
13 (Intercept) -67.4      18.5        -3.64 3.73e- 4
14 环球彩票Win
d         -2.94      0.531      -5.53 1.41e- 7
15 Solar.R       0.0668    0.0189      3.54 5.40e- 4
16 Temp          1.61      0.204       7.87 6.83e-13
17 (Intercept) -38.4      21.6        -1.78 7.69e- 2
18 环球彩票Win
d         -3.99      0.620      -6.44 1.52e- 9
19 Solar.R       0.0523    0.0223      2.34 2.05e- 2
20 Temp          1.44      0.239       6.04 1.15e- 8

多元线性模型y=a*x1+b*x2+c*x3+e,x1=环球彩票Win d,x2=Solar.R,x3=Temp, e=Intercept,上面数据集以每4行为一组,每组分别为多元线性回归的结果。环球彩票关于 多元线性回归的详细介绍,请参考文章R语言解读多元线性回归模型

数据解释

  • term,变量
  • estimate ,参数估计值
  • std.error,残差的标准误差
  • statistic,t统计量
  • p.value,p-value值

4.4 模型评估

使用pool()函数,将5个单独回归模型进行汇总,整理模型的标准误和p值等多个评价参数,进行F检验,判断模型整体拟合的有效性。pooled对象是一个包含这m个统计分析平均结果的列表对象,准确地反映出由于缺失值和多重插补而产生的不确定性。


> pooled=pool(fit)
> pooled
Class: mipo    m = 5 
                estimate         ubar            b            t dfcom        df        riv
(Intercept) -54.87825512 3.864370e+02 2.189907e+02 6.492258e+02   149  19.08909 0.68003032
环球彩票Win
d         -3.37959755 3.186764e-01 1.758552e-01 5.297026e-01   149  19.61464 0.66219610
Solar.R       0.05833371 4.121520e-04 2.745861e-05 4.451024e-04   149 114.74908 0.07994702
Temp          1.54135605 4.758977e-02 1.736383e-02 6.842637e-02   149  30.33948 0.43783767
                lambda        fmi
(Intercept) 0.40477265 0.45866597
环球彩票Win
d        0.39838627 0.45159195
Solar.R     0.07402865 0.08975652
Temp        0.30451120 0.34623283

指标解读:

  • mipo,一种数据结构,pool用于产生mipo对象
  • m, m=5 多组插补数据
  • estimate,估计值,pool汇总的完整数据估计
  • ubar, 估计的内插方差,std.error的平方的均值,越小越好
  • b,估计之间的插补方差, estimate的方差,越小越好
  • t,估计的总方差,ubar + (1 + 1 / m) * b, 越小越好
  • dfcom,完整数据的自由度
  • df,t统计量的自由度,由Barnard-Rubin提供一种自由度的计算环球彩票方法 ,越小越好
  • riv,方差的相对增加,(1 + 1 / m) * b / ubar,越小越好
  • lambda, 缺失值的比例,(1 + 1 / m) * b / t,越小越好
  • fmi, 缺少信息的部分,(riv + 2 / (df + 3)) / (riv + 1)),越小越好

对插补的结果,再进行R^2检验。


> pool.r.squared(fit)
          est     lo 95     hi 95 fmi
R^2 0.5630547 0.4404848 0.6677567 NaN

指标解读:

  • est, 合并的R^2,估计值
  • lo 95,合并的R^2,95%下限
  • hi 95,合并的R^2,95%上限
  • fmi,缺少信息的部分,越小越好

pool的过程,主要就是完成F检验和R^2检验。

3.4 最终结果

最后,根据检验统计量,环球彩票环球彩票我 们 选择以第2次pmm的结果做为填充,获得填充后的数据集。


> completeData <- complete(imputed_Data,2)
> completeData
    Ozone Solar.R 环球彩票Win
d Temp Month Day
1      41     190  7.4   67     5   1
2      36     118  8.0   72     5   2
3      12     149 12.6   74     5   3
4      18     313 11.5   62     5   4
5      18     299 14.3   56     5   5
6      28     157 14.9   66     5   6
7      23     299  8.6   65     5   7
8      19      99 13.8   59     5   8
9       8      19 20.1   61     5   9
10     30     194  8.6   69     5  10

数据和填充过程,如下图所示。

mice包提供了29种填充算法,上面的例子环球彩票我 只是用到了pmm一种环球彩票方法 。大家可以根据不同的业务场景要求和数据的分布情况,用不同环球彩票方法 进行填充,让数据差异变化最小。


> methods(mice)
 [1] mice.impute.2l.bin       mice.impute.2l.lmer      mice.impute.2l.norm      mice.impute.2l.pan      
 [5] mice.impute.2lonly.mean  mice.impute.2lonly.norm  mice.impute.2lonly.pmm   mice.impute.cart        
 [9] mice.impute.jomoImpute   mice.impute.lda          mice.impute.logreg       mice.impute.logreg.boot 
[13] mice.impute.mean         mice.impute.midastouch   mice.impute.norm         mice.impute.norm.boot   
[17] mice.impute.norm.nob     mice.impute.norm.predict mice.impute.panImpute    mice.impute.passive     
[21] mice.impute.pmm          mice.impute.polr         mice.impute.polyreg      mice.impute.quadratic   
[25] mice.impute.rf           mice.impute.ri           mice.impute.sample       mice.mids               
[29] mice.theme       

缺失数据处理,通常都是一种复杂的不确定问题,非常棘手。通过使用mice包,环球彩票环球彩票我 们 可以高效地完成复杂的缺失数据填充,把数据缺失问题转化成统计问题,用统计量来评估插值对于原数据集的影响。非常强大的环球彩票工具 ,让环球彩票环球彩票我 们 少拍脑袋,科学处理数据问题。

转载请注明出处:
http://xiyueit.com/r-na-mice

打赏作者

用R语言解读牛顿冷却定律

R的极客理想系列文章,涵盖了R的思想,使用,环球彩票工具 ,创新等的一系列要点,以环球彩票我 个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,环球彩票互联网 ….都在使用R语言。

要成为有理想的极客,环球彩票环球彩票我 们 不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让环球彩票环球彩票我 们 一起动起来吧,开始R的极客理想。

环球彩票关于 作者:

  • 张丹(Conan), 程序员/Quant: Java,R,Nodejs
  • blog: http://xiyueit.com
  • email: bsspirit@gmail.com

转载请注明出处:
http://xiyueit.com/r-newton-cooling

前言

牛顿冷却定律是经典的热力学模型,描述了温度与时间之间的指数衰减函数。环球彩票环球彩票我 们 通过对基础科学的研究,可以极大地开阔思路,并形成跨学科的解决方案,引用新的思考维度,解决现有场景的难题。站在物理学的制高点,推进数据科学的发展。

目录

  1. 牛顿冷却定律
  2. 数据集
  3. 程序实现

1. 牛顿冷却定律

牛顿冷却定律,是一个热力学的基本定律,指物体所损失的热的速率与物体和其周围环境间的温度差是成比例的。当物体表面与周围存在温度差时,单位时间从单位面积散失的热量与温度差成正比,比例系数称为热传递系数。牛顿冷却定律,是由英国物理学家艾萨克·牛顿爵士所提出的,是一个经验发现的规律。

当物体处于与周围环境不同的温度时,它将逐渐冷却或加热直至温度相等。环球彩票环球彩票我 们 每个人都在环球彩票生活中遇到过,比如煮水沏茶的过程,把凉水烧热到100度的高温,然后沏茶,再等一会儿,到达可以喝的温度。水从100度的高温,自然到达可以喝的温度,就是冷却的过程,即物体温度与周围环境温度热交换的过程。

公式:

公式解释:

  • T: 温度
  • t: 时间
  • T(t):t时刻的温度
  • Ta:代表室温,T(t)-Ta就是当前温度与室温之间的温差。当前温度高于室温时,T(t)-Ta > 0。
  • T0:初始温度
  • k:传热系数,常数,表示室温与降温速率之间的比例关系
  • dT/dt:温度和时间的微分形式,表示物体的温度随时间下降的速度

把上面微分方程进行数学变型,产生可以下面等式。

经过数学变型后,就容易进行计算了。由于k是一个常数,所以环球彩票环球彩票我 们 需要先计算出k值,再用k值计算出每一时刻的模型值,最后用模型值和测试值进行对比。

进行k值的公式变型:

2. 数据集

用数据说话,环球彩票环球彩票我 们 通过一个实验来验证牛顿冷去定律,把牛顿冷却定律计算的数值与真实数据的测量值进行拟合。实验数据来源:Newton’s Law of Cooling An Experimental Investigation

实验开始,准备3个烧杯加热至沸腾,让3个烧杯中的水自然冷却,记录60分钟每分钟3个烧杯中水的温度。3个烧杯水中水量不同,第一个烧杯放入100毫升水,第二个烧杯放入300毫升水,第三个烧杯放入800毫升水。每个烧杯都有自己的温度计,温度计在测量之间保持在烧杯中,确保不会出现温度滞后。每分钟测量每个烧杯中水的温度,环境温度为23°C。

开发环境所使用的系统环境

  • 环球彩票Win 10 64bit
  • R: 3.6.1 x86_64-w64-mingw32/x64 b4bit

记录3个烧杯温度的数据集,共4列,保存文件为Temperature.csv。

  • time_min:分钟,冷却时间,从0分钟到60分钟
  • temp_100ml:温度(摄氏度),100ml水随时间变化的温度
  • temp_300ml:温度(摄氏度),300ml水随时间变化的温度
  • temp_800ml:温度(摄氏度),800ml水随时间变化的温度

> dat
   time_min temp_100ml temp_300ml temp_800ml
1         0        100        100        100
2         1         95         95         96
3         2         82         91         95
4         3         79         87         92
5         4         74         84         90
6         5         70         81         88
7         6         67         78         85
8         7         65         76         83
9         8         61         73         80
10        9         59         71         78
11       10         57         70         76
12       11         56         68         75
13       12         54         66         74
14       13         52         64         73
15       14         51         63         71
16       15         50         61         70
17       16         49         60         68
18       17         48         58         66
19       18         47         58         66
20       19         45         56         65
21       20         45         55         63
22       21         44         55         62
23       22         43         54         61
24       23         42         53         60
25       24         42         52         60
26       25         41         51         59
27       26         41         50         58
28       27         40         49         56
29       28         39         48         56
30       29         38         48         55
31       30         38         47         54
32       31         38         46         53
33       32         38         46         52
34       33         37         45         52
35       34         36         45         51
36       35         36         45         50
37       40         34         42         47
38       45         33         40         45
39       50         31         38         43
40       55         30         37         41
41       60         29         36         40

观察数据集,用散点图画出的测试值的3种情况。


# 加载数据集
> dat<-read.csv("Temperature.csv")

# 以散点图观察数据
> plot(dat$time_min,dat$temp_100ml)
> points(dat$time_min,dat$temp_300ml,col="red")
> points(dat$time_min,dat$temp_800ml,col="blue")


图解释:

  • x轴为时间
  • y轴为维度
  • 黑色散点,100毫升水的温度变化
  • 红色散点,300毫升水的温度变化
  • 蓝色散点,800毫升水的温度变化

初步观察,3个烧杯水开始温度重合,都是100度。100毫升水随时间变化快,800毫升水随时间变化慢。

3. 程序实现

接下来,分别计算出3种情况下的水温冷却的理论值,通过牛顿冷却定律,并验证理论值与观测值之差的差别,是否能通过理论值来描述测量值。

3.1 用牛顿冷却定律计算100毫升水的温度

计算步骤:

  1. 计算1分钟内的k值
  2. 计算每个统计时间的k值,并求平均k值
  3. 根据平均k值,计算每个统计时间温度
  4. 对比理论值和测量值,进行可视化
  5. 计算理论值和测量值的均方根误差
  6. 进行理论值和测量值的统计检验
  7. 得出结论

第一步,计算1分钟内的k值。


# 室温 23度
> Ta <- 23

# 根据公式计算k值, k = ln((T(1)-Ta)/(T(0)-Ta))/t1
> k = log((dat$temp_100ml[2]-Ta)/(dat$temp_100ml[1]-Ta)) / dat$time_min[2]
> k
[1] -0.0671393

当t=1时(1分钟),传热系数k值为-0.0671393。

第二步,计算每个统计时间的k值,并求平均k值


# 每个统计时候的k值
> kv = log((dat$temp_100ml[-1]-Ta)/(dat$temp_100ml[1]-Ta)) / dat$time_min[-1]
> kv
 [1] -0.06713930 -0.13313399 -0.10615124 -0.10299495 -0.09873156 -0.09326930 -0.08659083
 [8] -0.08827741 -0.08447628 -0.08174449 -0.07702708 -0.07581818 -0.07511612 -0.07225721
[15] -0.06986457 -0.06785681 -0.06617233 -0.06476398 -0.06593489 -0.06263815 -0.06187062
[22] -0.06127605 -0.06084202 -0.05830694 -0.05813735 -0.05590129 -0.05594785 -0.05611488
[29] -0.05640535 -0.05452517 -0.05276630 -0.05111735 -0.05165903 -0.05231930 -0.05082446
[36] -0.04864775 -0.04536045 -0.04528728 -0.04359810 -0.04253410

# k平均值
> k<-mean(kv);k
[1] -0.06758501

第三步,根据平均k值,计算每个统计时间温度


# 计算理论温度
> Tt<- Ta + (dat$temp_100ml[1] - Ta)*exp(k*dat$time_min[-1])
> Tt
 [1] 94.96792 90.26469 85.86882 81.76024 77.92015 74.33103 70.97645 67.84111 64.91067 62.17173
[11] 59.61179 57.21915 54.98287 52.89273 50.93919 49.11331 47.40676 45.81174 44.32095 42.92759
[21] 41.62529 40.40809 39.27044 38.20714 37.21333 36.28446 35.41630 34.60487 33.84648 33.13764
[31] 32.47513 31.85591 31.27716 30.73624 30.23066 28.15726 26.67841 25.62362 24.87129 24.33470

第四步,对比理论值和测量值,进行可视化


# 把100ml的理论值和测量值合并数据框
> df100<-data.frame(dat[,1:2],Tt=c(100,Tt))

# 进行可视化,测试值为黑色点,理论值为红色线
> plot(df100$time_min,df100$temp_100ml,ylim=c(0,100))
> lines(df100$time_min,df100$Tt,col="red")

第五步,计算理论值和测量值的均方根误差


> sqrt(sum((df100$Tt - df100$temp_100ml)^2)/nrow(df100))
[1] 4.800323

第六步,进行理论值和测量值的统计检验。

对观测值和理论分别进行正常分步的检验,判断是否属于正态分布。


# 观测值检验
> shapiro.test(df100[,2])

	Shapiro-Wilk normality test

data:  df100[, 2]
W = 0.88026, p-value = 0.0004557

# 理论值检验
> shapiro.test(df100[,3])

	Shapiro-Wilk normality test

data:  df100[, 3]
W = 0.90219, p-value = 0.001925

通过shapiro-wilk检验,以0.05为显著性水平,p-value都小于0.05,所以拒绝原假设,观测值和预测值都不符合正态分布,也就不能使用T检验和F检验,对数据集是分析了。

卡方检验,用于检验2个独立样本的观测值与期望值的差距。


# 卡方检验
> x <- matrix(c(df100[,2],df100[,3]), ncol = 2)
> chisq.test(x)

	Pearson's Chi-squared test

data:  x
X-squared = 10.266, df = 40, p-value = 1

以0.05为显著性水平,p-value=1都大于0.05,所以不能拒绝原假设,观测值和理论值没有明显的差异性,用理论值来近似观测值的环球彩票方法 是合理的。

第七步,得出结论,使用牛顿冷却定律进行理论值计算与实际观测值没有明显的差异性,用理论值来近似观测值的环球彩票方法 是合理的,根绝牛顿冷却定律得出的理论值可以作为有效的热传递的计算环球彩票方法 。

再分别计算300毫升水和800毫升水的理论值和观测值之间的关系。

3.2 用牛顿冷却定律计算300毫升水的温度
300毫升水,计算理论值。


# 300毫升水
> kv = log((dat$temp_300ml[-1]-Ta)/(dat$temp_300ml[1]-Ta)) / dat$time_min[-1]

# k均值
> k<-mean(kv);k
[1] -0.0447057

# 计算每个统计时间的温度
> Tt<- Ta + (dat$temp_300ml[1] - Ta)*exp(k*dat$time_min[-1]);Tt
 [1] 96.63347 93.41413 90.33555 87.39156 84.57629 81.88411 79.30963 76.84771 74.49342 72.24207
[11] 70.08915 68.03036 66.06159 64.17888 62.37850 60.65682 59.01043 57.43601 55.93043 54.49067
[21] 53.11386 51.79725 50.53820 49.33420 48.18284 47.08182 46.02894 45.02208 44.05925 43.13852
[31] 42.25804 41.41606 40.61089 39.84092 39.10461 35.87873 33.29902 31.23605 29.58630 28.26701

# 合并理论值和观测值
> df300<-data.frame(dat[,c('time_min',"temp_300ml")],Tt=c(100,Tt))

# 均方根误差
> sqrt(sum((df300$Tt - df300$temp_300ml)^2)/nrow(df300))
[1] 3.711398

# 可视化
> plot(df300$time_min,df300$temp_300ml,ylim=c(0,100))
> lines(df300$time_min,df300$Tt,col="red")

3.3 用牛顿冷却定律计算800毫升水的温度
800毫升水,计算理论值


> kv = log((dat$temp_800ml[-1]-Ta)/(dat$temp_800ml[1]-Ta)) / dat$time_min[-1]

# k均值
> k<-mean(kv);k
[1] -0.03266545

# 计算每个统计时间的温度
> Tt<- Ta + (dat$temp_800ml[1] - Ta)*exp(k*dat$time_min[-1]);Tt
 [1] 97.52540 95.13032 92.81222 90.56862 88.39712 86.29540 84.26123 82.29244 80.38692 78.54263
[11] 76.75762 75.02998 73.35785 71.73947 70.17309 68.65706 67.18974 65.76959 64.39507 63.06473
[21] 61.77714 60.53093 59.32478 58.15738 57.02750 55.93394 54.87552 53.85111 52.85963 51.90001
[31] 50.97123 50.07230 49.20226 48.36018 47.54516 43.84653 40.70523 38.03729 35.77137 33.84689

# 合并理论值和观测值
> df800<-data.frame(dat[,c('time_min',"temp_800ml")],Tt=c(100,Tt))

# 均方根误差
> sqrt(sum((df800$Tt - df800$temp_800ml)^2)/nrow(df800))
[1] 2.237821

# 可视化
plot(df800$time_min,df800$temp_800ml,ylim=c(0,100))
lines(df800$time_min,df800$Tt,col="red")


3.4 总结
最后,把3种情况的评价指标整理到一起,形成一个表格。

水量 k平均值 均方根误差
100ml -0.06758501 4.800323
300ml -0.0447057 3.711398
800ml -0.03266545 2.237821

800ml的烧瓶的水,均方根误差最小,理论值与观测值拟合最好;同时,传热系数k值最小,说明水量越多传热越慢。100ml的烧杯的水,均方根误差最大,理论值与观测值拟合最差;同时,传热系数k值最大,说明水量越少传热越快。

牛顿冷却定律是一种热力学模型,通过温度与时间之间的函数关系,构建出了一个指数衰减的过程,这种衰减过程不仅能表达热传导的关系,还被应用到了如新闻热度排名领域。对基础科学的公理和定理的研究,可以极大地开阔研究数据的思路,形成跨学科的解决方案,下篇文章环球彩票我 将介绍牛顿冷却定律在基于用户投票的排名算法场景中的应用。

转载请注明出处:
http://xiyueit.com/r-newton-cooling

打赏作者

用R语言解读统计检验-卡方检验

R的极客理想系列文章,涵盖了R的思想,使用,环球彩票工具 ,创新等的一系列要点,以环球彩票我 个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,环球彩票互联网 ….都在使用R语言。

要成为有理想的极客,环球彩票环球彩票我 们 不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让环球彩票环球彩票我 们 一起动起来吧,开始R的极客理想。

环球彩票关于 作者:

  • 张丹(Conan), 程序员/Quant: Java,R,Nodejs
  • blog: http://xiyueit.com
  • email: bsspirit@gmail.com

转载请注明出处:
http://xiyueit.com/r-test-x2

前言

做统计分析R语言是最好的,R语言对于统计检验有非常好的支持。环球彩票我 会分7篇文章,分别介绍用R语言进行统计量和统计检验的计算过程,包括T检验F检验卡方检验P值KS检验AICBIC等常用的统计检验环球彩票方法 和统计量。

本文是第三篇卡方检验,卡方检验关注点在实际值与预测值的差距,对于预测类的模型给出了客观地评价环球彩票方法 。

目录

  1. 卡方检验介绍
  2. 数据集
  3. 四格表法推导过程
  4. 程序实现

1. 卡方检验介绍

卡方检验,又称χ2检验,是一种非参数检验,主要是比较两个以及两个以上样本率以及两个分类变量之间是否具有显著的相关性,其根本思想是统计样本的实际观测值与理论推断值之间的偏离程度。卡方检验,是由英国统计学家Karl Pearson在1900年首次提出的,在《哲学杂志》上发表。

卡方检验,实际观测值与理论推断值之间的偏离程度就决定卡方值的大小,如果卡方值越大,两个变量偏差程度越大;反之,卡方值越小两个变量偏差越小;如果两个变量的完全相等,卡方值就为0,表明观测值与理论值完全符合。

通常卡方检验的主要应用方向为:卡方拟合优度检验,卡方独立性检验。卡方检验可以分为成组比较(不配对资料)和个别比较(配对,或同一对象两种处理的比较)两类。卡方检验试用条件包括,随机样本数据,卡方检验的理论频数不能太小。

卡方检验的计算公式:


公式解释:

  • χ2,为统计量,用于衡量实际值与理论值的差异程度。
  • A,为实际值
  • T,为理论值
  • 自由度,df = n-1

卡方检验有3种推导过程:四格表法的卡方检验,行列表法的卡方检验,列联表法的卡方检验,本文将介绍四格表法的卡方检验。

2. 数据集

卡方检验,对于数据有比较严格的要求,所以环球彩票环球彩票我 们 需要先找到一个合适的数据集,作为测试数据集。

环球彩票环球彩票我 们 以“抛硬币”来作为卡方检验的使用场景。对于1次抛硬币来说,只有2种结果,正面是字,反面是头像,两种结果的出现概率相同都是50%,且每次抛硬币都是独立事件。那么,根据这个抛硬币的结果做出假设,原假设(H0)抛硬币出现正面和反面的概率都是50%,为验证这个假设成立,就是卡方检验。

开发环境所使用的系统环境

  • 环球彩票Win 10 64bit
  • R: 3.6.1 x86_64-w64-mingw32/x64 b4bit

接下来,用随机过程模拟抛硬币的过程,设置抛硬币100次,1为正面,0为反面。


> N=100                 # 抛硬币100次
> set.seed(1)           # 随机种子
> coin<-sample(x=c(0,1), prob=c(0.5,0.5),size = N, replace = TRUE)  # 模拟抛硬币
> coin                  # 打印结果
  [1] 1 1 0 0 1 0 0 0 0 1 1 1 0 1 0 1 0 0 1 0 0 1 0 1 1 1 1 1 0 1 1 0 1 1 0 0 0 1 0 1 0 0 0 0
 [45] 0 0 1 1 0 0 1 0 1 1 1 1 1 0 0 1 0 1 1 1 0 1 1 0 1 0 1 0 1 1 1 0 0 1 0 0 1 0 1 1 0 1 0 1
 [89] 1 1 1 1 0 0 0 0 1 1 0 0

统计抛硬币的结果


> table(coin)  # 统计0和1的出现次数
coin
 0  1 
48 52 

用可视化展示累积的正面出现概率


> cnt <- cumsum(coin)  # 积累求和,正面1出现的次数
> n<-1:N               # 次数向量
> rate<- cnt/1:N       # 累积的正面出现概率
> plot(n,rate,type="o",log="x",   # 画图,x轴做log变型
+      xlim=c(1,N), ylim=c(0.0,1.0), cex.axis=1.5,
+      main="Running Proportion of Heads",cex.axis=1.5)

接下来,环球彩票环球彩票我 们 就以coin作为数据集,进行卡方检验的数据集。

3. 四格表法推导过程

四格表法的卡方检验,用于进行两个样本率或两个样本构成比的比较。

环球彩票环球彩票我 们 可以把抛硬币的过程,分成观察组和期望组,观察组就是上面随机过程的结果,期望组是环球彩票环球彩票我 们 根据概率的公式给出来的,把数据填入下面表格。

正面 反面 合计
观察次数 52 48 100
期望次数 50 50 100
合计 102 98 200

通过简单的计算,环球彩票环球彩票我 们 得出观察组的正面概率为52%和期望组的正面概率为50%,两者的差别可能是误差导致,也有可能是受到抛硬币的次数影响的。

为了确定真实原因,环球彩票环球彩票我 们 先假设抛硬币是符合大数定理的,是不受抛硬币的次数影响的,即抛硬币的次数和概率无关,环球彩票环球彩票我 们 可以得出正面结果的实际是(52+50) / (52+48+50+50)= 0.51 = 51%。

可推到出:

正面 反面 合计
观察次数 100*0.51 =51 100*(1-0.51) = 49 100
期望次数 100*0.51 =51 100*(1-0.51) = 49 100
合计 102 98 200

如果抛硬币次数与正面结果的概率是独立无关的,那么四格表里的理论值和实际值差别应该会很小。

代入公式进行计算。


X^2 = (52- 51)^2 / 51 + (48 - 49)^2 / 49 + (50 - 51)^2 / 51 + (50 - 49)^2 / 49 = 0.08003201

当以四格表法进行卡方检验时,还有一些限制条件,当两个独立样本比较时,需要分成以下3种情况:

  1. 所有的理论数T≥5并且总样本量n≥40,用Pearson卡方进行检验.
  2. 如果理论数T<5但T≥1,并且n≥40,用连续性校正的卡方进行检验.
  3. 如果有理论数T<1或n<40,则用Fisher’s检验.

根据限制条件,由于样本数据量为200大于40,所以需要对原公式进行修正。

修正公式为:


X^2 = (abs(52- 51)-0.5)^2 / 51 + (abs(48 - 49)-0.5)^2 / 49 + (abs(50 - 51)-0.5)^2 / 51 + (abs(50 - 49)-0.5)^2 / 49 
= 0.020008

最后,算出卡方统计量为0.020008。下面环球彩票环球彩票我 们 再用R语言的程序来算一下,看看是否是相同的结果。

4. 程序实现

用程序实现卡方检验,可以直接使用chisq.test()函数。


# 合并观测值和期望值,形成矩阵
> s<-c(table(coin),c(50,50))
> x <- matrix(s, ncol = 2, byrow=TRUE)
> x
     [,1] [,2]
[1,]   48   52
[2,]   50   50

# 卡方检验
> chisq.test(x)

	Pearson's Chi-squared test with Yates' continuity correction

data:  x
X-squared = 0.020008, df = 1, p-value = 0.8875

指标解释:

  • H0:抛硬币出现正面和反面的概率都是50%,与次数无关
  • X-squared统计量:0.020008
  • df,自由度,1
  • p-value值:0.8875

结果解读,以0.05为显著性水平,自由度df=1,统计量X-squared = 0.020008 小于临界值0.455(查表),说明X-squared值不显著,不能拒绝原假设。以0.05为显著性水平,p-value=0.8875大于0.05,所以不能拒绝原假设,说明抛硬币的次数与结果为正面的概率是没有关系的。这个结果与环球彩票环球彩票我 们 构造的数据是一致的,也是符合大数定理的。

查看观察值和期望值。


# 查看观察值
> chi$observed
     [,1] [,2]
[1,]   48   52
[2,]   50   50

# 查看期望值
> chi$expected
     [,1] [,2]
[1,]   49   51
[2,]   49   51

# 残差
> chi$residuals
           [,1]      [,2]
[1,] -0.1428571  0.140028
[2,]  0.1428571 -0.140028

手动计算X^2值和P值,环球彩票关于 P值的详细解释,请查看文章R语言实现统计检验-P值


# 手动计算X^2值
> sr <- rowSums(x)
> sc <- colSums(x)
> n <- sum(x)
> E <- outer(sr, sc, "*")/n;E
     [,1] [,2]
[1,]   49   51
[2,]   49   51
> x2 <- sum((abs(x - E) - 0.5)^2/E);x2
[1] 0.020008

# 手动计算P值
> pchisq(x2, 1, lower.tail = FALSE)
[1] 0.8875147

本文用R语言详细介绍了卡方检验的计算过程,x^2越小,越能验证假设是正确的,越能说明观测值与期望值是一致的,x^2越大,越证明假设是错误的,预测越不准。为预测类的模型提供了评价标准。

转载请注明出处:
http://xiyueit.com/r-test-x2

打赏作者

用R语言解读统计检验-F检验

R的极客理想系列文章,涵盖了R的思想,使用,环球彩票工具 ,创新等的一系列要点,以环球彩票我 个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,环球彩票互联网 ….都在使用R语言。

要成为有理想的极客,环球彩票环球彩票我 们 不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让环球彩票环球彩票我 们 一起动起来吧,开始R的极客理想。

环球彩票关于 作者:

  • 张丹(Conan), 程序员/Quant: Java,R,Nodejs
  • blog: http://xiyueit.com
  • email: bsspirit@gmail.com

转载请注明出处:
http://xiyueit.com/r-test-f

前言

做统计分析R语言是最好的,R语言对于统计检验有非常好的支持。环球彩票我 会分7篇文章,分别介绍用R语言进行统计量和统计检验的计算过程,包括T检验F检验卡方检验P值KS检验AICBIC等常用的统计检验环球彩票方法 和统计量。

本文是第二篇F检验,T检验关注点在均值,而F检验关注点在方差。

目录

  1. F检验介绍
  2. 数据集
  3. F检验实现

1. F检验介绍

F检验法(F-test),初期叫方差比率检验(Variance Ratio),又叫联合假设检验(Joint Hypotheses Test),是英国统计学家Fisher提出的,主要通过比较两组数据的方差,以确定他们的密度是否有显著性差异。至于两组数据之间是否存在系统误差,则在进行F检验并确定它们的密度没有显著性差异之后,再进行T检验

F检验,是一种在零假设(H0)之下,统计值服从F-分布的检验。

样本标准偏差的平方公式:

F统计量计算公式:


公式解释

  • F:统计量,根据自由度查表,当F值小于查表值时没有显著差异,当F值大于等于查表值时有显著差异
  • S1:样本1的标准差
  • S2:样本2的标准差
  • 分子自由度: df=分子的数量-1
  • 分母自由度: df=分母的数量-1

T检验和F检验对比
T检验用来检测数据的准确度(系统误差),F检验用来检测数据的精密度(偶然误差)。在定量分析过程中,常遇到两种情况:一种是样本测量的平均值与真值不一致;另一种是两组测量的平均值不一致。

上述不一致是由于定量分析中的系统误差和偶然误差引起的,因此必须对两组分析结果的准确度或精密度是否存在显著性差异做出判断,两组数据的显著性检验顺序是先F检验后T检验。

T检验是检查两组均值的差异,而F检验是检查多组均值之间的差异。

对于多元线性回归模型,t检验是对于单个变量进行显著性,检验该变量独自对被解释变量的影响。f检验是检验回归模型的显著意义,即所有解释变量联合起来对被解释变量的影响,环球彩票关于 线性回归请参考文章,R语言解读一元线性回归模型R语言解读多元线性回归模型

2. 数据集

F检验,对于数据有比较严格的要求,所以环球彩票环球彩票我 们 需要先找到一个合适的数据集,作为测试数据集。环球彩票我 发现了R语言自带的一个数据集ToothGrowth,是很好的测试数据集,本文接下来的内容,将以这个数据集进行测试,来介绍F检验。

开发环境所使用的系统环境

  • 环球彩票Win 10 64bit
  • R: 3.4.2 x86_64-w64-mingw32/x64 b4bit

数据集ToothGrowth,记录了60只豚鼠的牙齿生长速度,使用2种不同的环球彩票方法 (OJ和VC),每天按3种不同的注射剂量进行注射,对牙齿的生长速度的对比数据,共3列,60条记录。

  • len列,为牙齿长度
  • supp列,为注射环球彩票方法
  • dose列,为注射剂量

查看数据集,打印前10行


> head(ToothGrowth,10)
    len supp dose
1   4.2   VC  0.5
2  11.5   VC  0.5
3   7.3   VC  0.5
4   5.8   VC  0.5
5   6.4   VC  0.5
6  10.0   VC  0.5
7  11.2   VC  0.5
8  11.2   VC  0.5
9   5.2   VC  0.5
10  7.0   VC  0.5

F检验对于数据的正态性非常敏感,环球彩票环球彩票我 们 需要先对选定数据集进行进行正态分布检验。使用Shapiro-Will作为正态分布检验的环球彩票方法 ,原假设H0:样本符合正态分布。


# 按不同的处理环球彩票方法
,进行分组
> len_VC<-ToothGrowth$len[which(ToothGrowth$supp=='VC')]
> len_OJ<-ToothGrowth$len[which(ToothGrowth$supp=='OJ')]

# 正态分布检验
> shapiro.test(len_VC)

	Shapiro-Wilk normality test

data:  len_VC
W = 0.96567, p-value = 0.4284

# 正态分布检验
> shapiro.test(len_OJ)

	Shapiro-Wilk normality test

data:  len_OJ
W = 0.91784, p-value = 0.02359

两个样本的W统计量都接近1,且p-value都大于0.05,不能拒绝原假设,两组样本数据为正态分布。

查看数据的相关性。


> coplot(len ~ dose | supp, data = ToothGrowth, panel = panel.smooth,
       xlab = "ToothGrowth data: length vs dose, given type of supplement")

3. F检验实现

3.1 随机数进行F检验
环球彩票环球彩票我 们 先用一种随机数,来做一下F检验。以正态分布生成2组数据,数量,均值,方差都不同,进行F检验。


# 生成随机数
> set.seed(1)
> x <- rnorm(50, mean = 0, sd = 2)
> y <- rnorm(30, mean = 1, sd = 1)

# 进行F检验
> var.test(x, y)

	F test to compare two variances

data:  x and y
F = 2.6522, num df = 49, denom df = 29, p-value
= 0.006232
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 1.332510 4.989832
sample estimates:
ratio of variances 
          2.652168 

指标解释:

  • H0:原假设2组样本的方差,无显著差异
  • F统计量:2.6522
  • num df,分子自由度,50-1=49
  • denom df,分每自由度,30-1=29
  • p-value值:0.006232
  • 95 percent confidence interval:95%的置信区间
  • ratio of variances:方差比率2.652168

结果解读,以0.05为显著性水平,F = 2.6522大于临界值1.81(查表),F值显著,拒绝原假设。以0.05为显著性水平,p-value=0.006232小于0.05,拒绝原假设,两样本方差有显著性差异。这个结果与环球彩票环球彩票我 们 构造的数据是一致的,样本的方差就是不同的。

3.2 ToothGrowth进行F检验
使用ToothGrowth数据集进行F检验,原假设HO,用VC和OJ两种环球彩票方法 按3种剂量进行注射,对于60只豚鼠的牙齿生长速度的方差,没有显著性差异。


> var.test(len_VC,len_OJ)

	F test to compare two variances

data:  len_VC and len_OJ
F = 1.5659, num df = 29, denom df = 29, p-value
= 0.2331
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.745331 3.290028
sample estimates:
ratio of variances 
          1.565937 

结果解读,以0.05为显著性水平,F=1.5659小于临界值1.90(查表),F值不显著,不能拒绝原假设。以0.05为显著性水平,p-value=0.2331大于0.05,不能拒绝原假设,所以两种环球彩票方法 的3种剂量实验的方差,没有显著性的差异。

环球彩票环球彩票我 们 可以用F值进行显著性差异判断,也可以用p值进行显著性差异判断,他们的作用是一样的。F值判断时,需要用计算所得的F值,与显著性水平查表对比。p值相当于是把F值,进行一种标准化的变型,只和已经定义好的显著性水平比就行了,比如0.05, 0.01, 0.001等几个固定值。

手动计算F值和P值,环球彩票关于 P值的详细解释,请查看文章R语言实现统计检验-P值


# 手动计算T值
> Xn<-length(len_VC)
> Yn<-length(len_OJ)
> Xm<-mean(len_VC)
> Ym<-mean(len_OJ)

# 计算两组样本的偏方差
> fx<-sum((len_VC-Xm)^2)/(Xn-1)
> fy<-sum((len_OJ-Ym)^2)/(Yn-1)

# 计算F值
> fx/fy
[1] 1.565937

# 手动计算P值,双边检验
> p_value<-pf(f_stat,Yn-1,Xn-1)
> p_value<-2*min(p_value, 1 - p_value);p_value
[1] 0.2331433

用F检验测试样本数据的偶然误差,对数据集进行方差齐性检验,从而判断数据是否有显著性差异,为方差分析提供了基本的判别环球彩票方法 ,对于研究数据的波动性是非常有用的。

转载请注明出处:
http://xiyueit.com/r-test-f

打赏作者