标签:value 距离 line one mamicode 没有 cal call 压力
第二章
2.12
(1)拟合模型:
> library(openxlsx) #加载包 openxlsx > data = read.xlsx("22_data.xlsx",sheet = 2) #read.xlsx 函数读入数据 > x = data[,1] > y = data[,2] > res = lm(y~x) #构造线性回归模型函数 > res #结果 Call: lm(formula = y ~ x) Coefficients: (Intercept) x #得出线性回归模型 y = -6.332 + 9.208 x -6.332 9.208 > summary(res) #打印方差分析,系数的估计值及其检验。 Call: lm(formula = y ~ x) Residuals: #残差统计量,残差第一四分位数(1Q)和第三分位数(3Q)有大约相同的幅度,意味着有较对称的钟形分布 Min 1Q Median 3Q Max -2.5629 -1.2581 -0.2550 0.8681 4.0581 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -6.33209 1.67005 -3.792 0.00353 ** #截距的点估计值及其检验 x 9.20847 0.03382 272.255 < 2e-16 *** #自变量系数的点估计值及其检验 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.946 on 10 degrees of freedom Multiple R-squared: 0.9999, Adjusted R-squared: 0.9999 #相关系数与调整的相关系数 F-statistic: 7.412e+04 on 1 and 10 DF, p-value: < 2.2e-16 #模型的显著性检验(F检验)
(2)根据上面程序结果,自变量具有显著性,模型具有显著性。
(3) 不能支持管理员的观点,根据构造的线性回归模型,平均环境温度增加1°F,平均月水蒸气消耗量将增加 9208+lb ,达不到10000lb.
(4) 使用58°F的平均环境温度构造一个月中水蒸气消耗量的99%预测区间:
> library(openxlsx) > data = read.xlsx("22_data.xlsx",sheet = 2) > x = data[,1] > y = data[,2] > fun = function(x) #计算预测值函数 + { + y = -6.332 + 9.208*x + } > y_pred = fun(x) #计算所有数的预测值 > s_y0_pred = function(x0,x,y,n) #构造计算预测值标准差的函数 + { + n = 12 + y_pred = fun(x) + sse = sum((y_pred - y)*(y_pred - y)) + se = sqrt(sse/(n-2)) + se * sqrt(1/n + ((x0-mean(x))^2)/sum((x-rep(mean(x),n))^2)) + } > x0 = 58 ; n = 12 > y0_pred = fun(x0) #当环境温度为58°F,对应的因变量预测值 > s = s_y0_pred(x0,x,y,n) > print(c(y0_pred-qt(0.995,n-2)*s,y0_pred+qt(0.995,n-2)*s)) #输出结果 [1] 525.5666 529.8974
2.13
a.做出散点图
> data = read.xlsx("22_data.xlsx",sheet = 1) > x = data[,2] > y = data[,1] > plot(x,y,main = "散点图",xlab = "index",ylab = "days")
> abline(lm(y~x))
b.估计预测方程
> lm(y~x) Call: lm(formula = y ~ x) Coefficients: (Intercept) x -193.0 15.3
预测方程为:y = *-193.0 + 15.3 x
c.进行回归显著性检验
> summary(lm(y~x)) Call: lm(formula = y ~ x) Residuals: Min 1Q Median 3Q Max -41.70 -21.54 2.12 18.56 36.42 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -192.984 163.503 -1.180 0.258 #p值大于0.05 x 15.296 9.421 1.624 0.127 #p值大于0.05 , 回归变量与响应变量没有显著相关性 Residual standard error: 23.79 on 14 degrees of freedom Multiple R-squared: 0.1585, Adjusted R-squared: 0.09835 F-statistic: 2.636 on 1 and 14 DF, p-value: 0.1267
根据上述结果,指数与天数并没有显著相关性。
d.计算并画出95%置信带与95%预测带
> sx = sort(x)
> #计算置信区间 > conf = predict(fm,data.frame(x = sx),interval = "confidence") > #计算预测区间 > pred = predict(fm,data.frame(x=sx),interval = "prediction") > plot(x,y,ylim = c(0,150),xlab = "index",ylab = "days",main = "95%预测带、置信带") > abline(fm) > lines(sx,conf[,2],col = "red") > lines(sx,conf[,3],col = "red") > lines(sx,pred[,2],col = "blue") > lines(sx,pred[,3],col = "blue")
2.14
a.散点图
b.估计预测方程
> fm = lm(y~x) > fm Call: lm(formula = y ~ x) Coefficients: (Intercept) x 0.6714 -0.2964
预测方程为:y = 0.6714 - 0.2964 x
c.数据分析
> summary(fm) Call: lm(formula = y ~ x) Residuals: Min 1Q Median 3Q Max -0.20464 -0.10634 0.02196 0.08527 0.20643 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.6714 0.1595 4.209 0.00563 ** x -0.2964 0.2314 -1.281 0.24754 #p值大于0.05 ,该自变量没有显著相关 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.15 on 6 degrees of freedom Multiple R-squared: 0.2147, Adjusted R-squared: 0.08382 R^2 = 0.2147 F-statistic: 1.64 on 1 and 6 DF, p-value: 0.2475 #整个模型不具有显著性。
d.计算并画出95%置信带和95%预测带
> plot(x,y,main = "散点图",xlab = "比率",ylab = "黏度",ylim = c(-0.1,1)) > sx = sort(x) > conf = predict(fm,data.frame(sx),interval = "confidence") > pred = predict(fm,data.frame(sx),interval = "prediction") > abline(fm) > lines(sx,conf[,2],col = "red") #绘制置信下限 > lines(sx,conf[,3],col = "red") #绘制置信上限 > lines(sx,pred[,2],col = "blue") #绘制预测下限 > lines(sx,pred[,3],col = "blue") #绘制预测上限
2.15
a.估计预测方程
Call: lm(formula = y ~ x) Coefficients: (Intercept) x 1.281511 -0.008758
预测方程为: y = 1.281511 - 0.008758 x
b.全面分析此模型
> fm = lm(y~x) > summary(fm) Call: lm(formula = y ~ x) Residuals: Min 1Q Median 3Q Max -0.043955 -0.035863 -0.009305 0.019900 0.069559 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.2815107 0.0468683 27.34 1.58e-07 *** x -0.0087578 0.0007284 -12.02 2.01e-05 *** #根据 p 值,自变量温度极显著 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.04743 on 6 degrees of freedom Multiple R-squared: 0.9602, Adjusted R-squared: 0.9535 F-statistic: 144.6 on 1 and 6 DF, p-value: 2.007e-05 #根据 p 值,整个回归模型是显著的
c.画95%置信带、预测带
2.16
先画出散点图:
从散点图可以看出容量与压力之间具有明显的线性关系,我们构造一元线性模型:
> fm = lm(y~x) > fm Call: lm(formula = y ~ x) Coefficients: (Intercept) x -290.707 2.346
估计预测模型为: y = -290.707 + 2.346x
再对模型进行检验:
> summary(fm) Call: lm(formula = y ~ x) Residuals: Min 1Q Median 3Q Max -4.3276 -0.9227 0.0773 1.2676 2.9577 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -2.907e+02 1.355e+00 -214.6 <2e-16 *** x 2.346e+00 4.007e-04 5855.4 <2e-16 *** #该自变量极显著 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.741 on 31 degrees of freedom Multiple R-squared: 1, Adjusted R-squared: 1 F-statistic: 3.429e+07 on 1 and 31 DF, p-value: < 2.2e-16 #整个回归模型极显著
2.17
> x = data[,2] > y = data[,1] > n = length(x) > plot(x,y) > > fm = lm(y~x) #一元回归模型 > abline(fm) > coef(fm) (Intercept) x 163.930734 1.579647 > summary(fm) Call: lm(formula = y ~ x) Residuals: Min 1Q Median 3Q Max -1.41483 -0.91550 -0.05148 0.76941 2.72840 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 163.9307 2.6551 61.74 < 2e-16 *** x 1.5796 0.1051 15.04 1.88e-10 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.291 on 15 degrees of freedom Multiple R-squared: 0.9378, Adjusted R-squared: 0.9336 F-statistic: 226 on 1 and 15 DF, p-value: 1.879e-10 > anova(fm) #方差分析 Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(>F) x 1 376.92 376.92 226.04 1.879e-10 *** Residuals 15 25.01 1.67 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > > rm(list = ls())
2.18
library(openxlsx) data = read.xlsx("2.18.xlsx",sheet = 1) x = data[,2] y = data[,3] n = length(x) plot(x,y) fm = lm(y~x) #一元回归 coef(fm) #输出回归系数 summary(fm) anova(fm) #构造此数据的95%置信带与预测带 sx = sort(x) conf = predict(fm,data.frame(x = sx),interval = "confidence") pred = predict(fm,data.frame(x = sx),interval = "prediction") abline(fm) lines(sx,conf[,2],col = ‘red‘) lines(sx,conf[,3],col = ‘red‘) lines(sx,pred[,2],col = ‘blue‘) lines(sx,pred[,3],col = ‘blue‘) plot(x,y,main = "%95置信带与95%预测带",xlab = "花费钱数",ylab = "每周挣回的印象",ylim=c(-100,200)) rm(list = ls())
第三章
3.8
#a.拟合co2产量y与总溶剂量x6和氢消耗量x7关系的多元回归模型
library(openxlsx) data = read.xlsx("3.8.xlsx",sheet = 1) data y = data[,1] #响应变量 x = data[,c(7,8)] #回归变量 fm = lm(y~.,x) #多元线性回归 summary(fm) anova(fm) #检验显著性 summary(fm) #d confint(fm) #e x6 = data[,7] fm1 = lm(y~x6) summary(fm1) anova(fm1) confint(fm1,level = 0.95) rm(list = ls())
3.9
library(openxlsx) data = read.xlsx("3.9.xlsx",sheet = 1) y = data[,1] x = data[,c(2,5)] #a.拟合多元回归模型 fm = lm(y~.,x) coef(fm) #b,c 检验回归显著性() anova(fm) summary(fm) #e #检验是否有潜在的多重共线性 r2 = 0.6367 vif = 1/(1-r2) rm(list = ls())
3.10
> #3.10 > library(openxlsx) > data = read.xlsx("3.10.xlsx",sheet = 1) > #a > y = data[,6] > x = data[,c(1:5)] > fm = lm(y~.,x) > coef(fm) (Intercept) Clarity Aroma Body Flavor 3.9968648 2.3394535 0.4825505 0.2731612 1.1683238 Oakiness -0.6840102 > #b,c > summary(fm) Call: lm(formula = y ~ ., data = x) Residuals: Min 1Q Median 3Q Max -2.85552 -0.57448 -0.07092 0.67275 1.68093 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.9969 2.2318 1.791 0.082775 . Clarity 2.3395 1.7348 1.349 0.186958 Aroma 0.4826 0.2724 1.771 0.086058 . Body 0.2732 0.3326 0.821 0.417503 Flavor 1.1683 0.3045 3.837 0.000552 *** Oakiness -0.6840 0.2712 -2.522 0.016833 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.163 on 32 degrees of freedom Multiple R-squared: 0.7206, Adjusted R-squared: 0.6769 F-statistic: 16.51 on 5 and 32 DF, p-value: 4.703e-08 > #d > xx = data[,c(2,4)] > fm1 = lm(y~.,xx) > summary(fm1) Call: lm(formula = y ~ ., data = xx) Residuals: Min 1Q Median 3Q Max -2.19048 -0.60300 -0.03203 0.66039 2.46287 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 4.3462 1.0091 4.307 0.000127 *** Aroma 0.5180 0.2759 1.877 0.068849 . Flavor 1.1702 0.2905 4.027 0.000288 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.229 on 35 degrees of freedom Multiple R-squared: 0.6586, Adjusted R-squared: 0.639 F-statistic: 33.75 on 2 and 35 DF, p-value: 6.811e-09 > AIC(fm) #优先考虑的模型应是AIC值最小的那一个 [1] 126.7552 > AIC(fm1) [1] 128.3761 > #e > conf = confint(fm) > conf1 = confint(fm1) > conf = as.matrix(conf) > conf1 = as.matrix(conf1) > > conf[5,2]-conf[5,1] [1] 1.240414 > conf1[3,2]-conf[3,1] [1] 1.83241 > > rm(list = ls()) >
3.11
> #3.11 > library(openxlsx) > data = read.xlsx("3.11.xlsx",sheet = 1) > y = data[,6] > x = data[,c(1:5)] > #a > fm = lm(y~.,x) > coef(fm) (Intercept) x1 x2 x3 5.207905e+01 5.555556e-02 2.821429e-01 1.250000e-01 x4 x5 1.776357e-16 -1.606498e+01 > #b,c > summary(fm) Call: lm(formula = y ~ ., data = x) Residuals: Min 1Q Median 3Q Max -12.250 -4.438 0.125 5.250 9.500 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 5.208e+01 1.889e+01 2.757 0.020218 * x1 5.556e-02 2.987e-02 1.860 0.092544 . x2 2.821e-01 5.761e-02 4.897 0.000625 *** x3 1.250e-01 4.033e-01 0.310 0.762949 x4 1.776e-16 2.016e-01 0.000 1.000000 x5 -1.606e+01 1.456e+00 -11.035 6.4e-07 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 8.065 on 10 degrees of freedom Multiple R-squared: 0.9372, Adjusted R-squared: 0.9058 F-statistic: 29.86 on 5 and 10 DF, p-value: 1.055e-05 > #d > xx = data[,c(2,5)] > fm1 = lm(y~.,xx) > summary(fm1) Call: lm(formula = y ~ ., data = xx) Residuals: Min 1Q Median 3Q Max -15.375 -4.188 -0.875 3.438 12.625 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 80.13461 5.69146 14.080 3.01e-09 *** x2 0.28214 0.05883 4.796 0.000349 *** x5 -16.06498 1.48659 -10.807 7.26e-08 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 8.236 on 13 degrees of freedom Multiple R-squared: 0.9149, Adjusted R-squared: 0.9018 F-statistic: 69.89 on 2 and 13 DF, p-value: 1.107e-07 > AIC(fm) [1] 118.6885 > AIC(fm1) [1] 117.5552 > > #e > confint(fm) 2.5 % 97.5 % (Intercept) 9.99688896 94.1612109 x1 -0.01100273 0.1221138 x2 0.15378045 0.4105053 x3 -0.77353688 1.0235369 x4 -0.44926844 0.4492684 x5 -19.30879739 -12.8211665 > confint(fm1) #温度:x2 2.5 % 97.5 % (Intercept) 67.8389462 92.4302647 x2 0.1550559 0.4092298 x5 -19.2765650 -12.8533989 >
第四章
例4.1 根据例3.1数据,输出残差,标准化残差,学生化残差,press残差,外部学生化残差 表格
library(openxlsx) #例4.1--------------------------------------------------------- #处理数据 data = read.xlsx("3.1.xlsx",sheet = 1) data = data[,c(2,3,4)] names(data)=c("time","cases","distance") y = data$time x1 = data$cases x2 = data$distance #线性回归 fm = lm(y~x1+x2) #残差 ei = residuals(fm) View(ei) #标准化残差(1) di = rstandard(fm) View(di) #计算mse的函数 mse = function(ei,p) #ei是残差向量,p是回归参数个数 { n = length(ei) sse = sum(ei**2) mse = sse/(n-p) return(mse) } di_ = ei/sqrt(mse(ei,3))#标准化残差(2) View(di_) #学生化残差(1) ri = rstudent(fm) View(ri) #计算帽子矩阵,并提取对角线元素 H = function(X) #X是回归向量矩阵 { h = X%*%solve(t(X)%*%X)%*%t(X) hii = diag(h) return(hii) } X = cbind(1,x1,x2) hii = H(X) #计算hii View(hii) ri_ = ei/sqrt(mse(ei,3)*(1-hii)) #学生化残差(2) View(ri_) #计算PRESS统计量 e_i = ei/(1-hii) #计算e(i) View(e_i) #外部学生化残差 ti = function(ei,X) #输入残差回归变量矩阵 { p = ncol(X) #回归参数个数 n = length(ei) #数据个数 hii = H(X) #帽子矩阵主对角线元素 s2_i = ((n-p)*mse(ei,p) -(ei**2)/(1-hii)) / (n-p-1) #计算S(i)^2 ans = ei / sqrt(s2_i*(1-hii)) return(ans) } ti = ti(ei,X) View(ti) #计算PRESS统计量 press = function(ei,X) { hii = H(X) res = sum((ei/(1-hii))**2) #View(res) } Press = (ei/(1-hii))**2 View(Press) PRESS = press(ei,X) #输出PRESS统计量 #将所有残差数据写入表格 Num = seq(1,length(ei)) mydata = cbind(Num,ei,di_,ri_,hii,e_i,ti,Press) class(mydata) View(mydata) write.xlsx(mydata,"C:\\Users\\86130\\Desktop\\mydata.xlsx")
#例4.2 ti #外部学生化残差 View(ti) n = length(ti) #数据个数 order = rank(ti) #rank函数返回ti按升序排序之后的序号 Pi = (order-1/2)/n #累积概率 plot(ti,Pi,xlim=c(-3,5)) #画正态概率图 fm_tP = lm(Pi~ti) #线性回归模型 abline(fm_tP) #添加回归线 #例4.3 #画残差与拟合值y_i的残差图 plot(fitted(fm),ti) #fitted(fm)返回模型fm的预测值 abline(h = 0) #添加直线y=0 #例4.4 #画残差与回归变量的残差图 par(mfrow =c(1,2)) plot(x1,ti,xlab = "箱数",ylab = "外部学生化残差") abline(h=0) #h:y轴 v:x轴 plot(x2,ti,xlab = "距离",ylab = "外部学生化残差") abline(h=0) #例4.5 #画偏回归图 #回归变量x1的偏回归图 lm.y_x2 = lm(y~x2) lm.x1_x2 = lm(x1~x2) plot(resid(lm.x1_x2),resid(lm.y_x2),xlab = "箱数",ylab = "时间") #回归变量x2的偏回归图 lm.y_x1 = lm(y~x1) lm.x2_x1 = lm(x2~x1) plot(resid(lm.x2_x1),resid(lm.y_x1),xlab = "距离",ylab = "时间",pch = 10) #例4.6 #计算PRESS的预测R^2 R_pred = function(X,y) { hii = H(X) ei = resid(lm(y~X[,2]+X[,3])) PRESS = sum((ei/(1-hii))**2) sst = sum((y-mean(y))**2) ans = 1-PRESS/sst return(ans) } R_pred(X,y) #例4.7 data = read.xlsx("2.1.xlsx",sheet = 1) names(data) = c("order","y","x") x = data$x y = data$y X = cbind(1,x) fm = lm(y~x) #绘制正态概率图 plot_ZP = function(ti) #输入外部学生化残差 { n = length(ti) order = rank(ti) #按升序排列,t(i)是第order个 Pi = (order-1/2)/n #累积概率 plot(ti,Pi,xlim=c(-3,3),xlab = "学生化残差",ylab = "百分比") #画正态概率图 } ei = resid(fm) ti = ti(ei,X) #计算外部学生化残差ti plot_ZP(ti) #绘制正态概率图 plot(fitted(fm),ti) #绘制残差与所预测y_pred的残差图 abline(h = 0) #绘制除去5,6两点的正态概率图 data = data[-c(5,6),] x = data$x y = data$y X = cbind(1,x) fm1 = lm(y~x) #线性模型 ei = resid(fm1) ti = ti(ei,X) #计算外部学生化残差ti plot(fitted(fm1),ti) #绘制残差与所预测y_pred的残差图 abline(h = 0) #例4.8 data = read.xlsx("4.8.xlsx",sheet = 1) x = data$x y = data$y fm = lm(y~x) #线性回归 a = anova(fm) #方差分析 sst = sum(a[2]) #总平方和 ssg = a[1,2] #回归平方和 sse = a[2,2] #残差平方和 level_x = c(table(x)>1) #查看哪些自变量重复 #进行失拟检验 library(rsm) #加载rsm包用于失拟检验 lm.rsm<-rsm(y~FO(x)) loftest(lm.rsm) #调用失拟检验函数loftest rm(list = ls())
例4.10 通过近邻点估计纯误差
#例4.10 data = read.xlsx("3.1.xlsx",sheet = 1) #导入数据 names(data)=c("order","time","cases","distance") y = data$time #准备数据 x1 = data$cases x2 = data$distance fm = lm(y~x1 + x2) #线性回归 coef(fm) b_cases = coef(fm)[2] #beta1 b_distance = coef(fm)[3] #beta2 y_pred = predict(fm) #计算预测值 ei = resid(fm) #残差 new_data = cbind(data,y_pred,ei) #构建新数据 new_data = new_data[order(new_data$y_pred),] #按照预测值升序排序 a = anova(fm) #方差分析 sse = a[3,2] #残差平方和 mse = a[3,3] #残差均方和 #计算Dii‘ Di_i = function(i,i_,mse,beta1,beta2,new_data) #i第i个点,i_第i_个点,data数据集 { one = beta1*(new_data$cases[i]-new_data$cases[i_])/sqrt(mse) two = beta2*(new_data$distance[i]-new_data$distance[i_])/sqrt(mse) ans = one**2 + two**2 return(ans) } #定义一个数据框用来存储结果 σ_ans = data.frame( i = numeric(0), #观测值i i_ = numeric(0), #观测值i_ Dii = numeric(0), #Di_i delta = numeric(0) #E|ei-ei_| ) #计算相邻k个点的两点的 Di_i,i,i_,Delta残差 for (k in c(1:4)) { for (i in c(1:24)) { if (i+k>25) break D = Di_i(i,i+k,mse,b_cases,b_distance,new_data) #计算相邻k个点的两点的Di_i E = abs(new_data$ei[i]-new_data$ei[i+k]) #计算相邻k个点的两点的Delta残差 another = data.frame( i = new_data$order[i], i_ = new_data$order[i+k], Dii = D, delta = E ) σ_ans = rbind(σ_ans,another) #合并两个数据框 } } names(σ_ans) = c("i","i_","Dii^2","Delta") #重命名最后的数据框 σ_ans = σ_ans[order(σ_ans$Dii^2),] #按照Di_i对数据框进行排序 row.names(σ_ans) = c(1:90) #对每一行进行编号 #计算累计标准差 std = numeric(0) #存储累计标准差 sum_Delta = 0 #存储累计Delta残差 for (i in 1:90) { sum_Delta = sum_Delta + σ_ans$Delta[i] #0.886/m*Σ(Delta) res = 0.886/i*sum_Delta std[i] = res } σ_ans = cbind(std,σ_ans)
4.16
#######################自己编的函数,运行后直接调用####################### #计算mse的函数 mse = function(ei,p) #ei是残差向量,p是回归参数个数 { n = length(ei) sse = sum(ei**2) mse = sse/(n-p) return(mse) } #计算帽子矩阵,并提取对角线元素 H = function(X) #X是回归向量矩阵 { h = X%*%solve(t(X)%*%X)%*%t(X) hii = diag(h) return(hii) } #外部学生化残差 ti = function(ei,X) #输入残差回归变量矩阵 { p = ncol(X) #回归参数个数 n = length(ei) #数据个数 hii = H(X) #帽子矩阵主对角线元素 s2_i = ((n-p)*mse(ei,p) -(ei**2)/(1-hii)) / (n-p-1) #计算S(i)^2 ans = ei / sqrt(s2_i*(1-hii)) return(ans) } #计算PRESS统计量 press = function(ei,X) #X是自变量的设计矩阵 { hii = H(X) res = sum((ei/(1-hii))**2) #View(res) } #计算PRESS的预测R^2 R_pred = function(X,y) { hii = H(X) ei = resid(lm(y~X[,2]+X[,3])) PRESS = sum((ei/(1-hii))**2) sst = sum((y-mean(y))**2) ans = 1-PRESS/sst return(ans) } #绘制正态概率图 plot_ZP = function(ti) #输入外部学生化残差 { n = length(ti) order = rank(ti) #按升序排列,t(i)是第order个 Pi = (order-1/2)/n #累积概率 plot(ti,Pi,xlim=c(-3,3),xlab = "学生化残差",ylab = "百分比") #画正态概率图 } #进行失拟检验 library(rsm) #加载rsm包用于失拟检验 lm.rsm<-rsm(y~FO(x)) loftest(lm.rsm) #调用失拟检验函数loftest #计算Dii‘ Di_i = function(i,i_,mse,beta1,beta2,new_data) #i第i个点,i_第i_个点,data数据集 { one = beta1*(new_data$cases[i]-new_data$cases[i_])/sqrt(mse) two = beta2*(new_data$distance[i]-new_data$distance[i_])/sqrt(mse) ans = one**2 + two**2 return(ans) }
#4.16 #a.残差的正态概率图 data = read.xlsx(‘3.12.xlsx‘,sheet = 1) #导入数据 y = data$y x1 = data$x1 x2 = data$x2 X = cbind(1,x1,x2) #处理数据 fm = lm(y~x1+x2) #线性回归模型 ei = resid(fm) #计算残差 ti = ti(ei,X) #ti()自制求外部学生化残差函数 plot_ZP(ti) #plot_zp()自制绘制正态概率图函数 #为什么要编写函数? #1.这些题目都是重复的代码操作 #2.如果是想多次重复打代码来熟悉,大可不必,因为会忘的。 #正态概率图有一个异常点,order(ti) 返回第一小的是第28号点 #b.残差与响应变量预测值的残差图 plot(fitted(fm),ti) #残差图表明残差包含在一条水平带中,模型不存在明显的缺点。 #c. #模型fm的PRESS统计量 press_fm = press(ei,X) #新模型fm1的PRESS统计量 fm1 = lm(y~x2) ei = resid(fm1) X = cbind(1,x2) press_fm1 = press(ei,X) #press()自制求press统计量函数 #选择press统计量小的模型
标签:value 距离 line one mamicode 没有 cal call 压力
原文地址:https://www.cnblogs.com/jiaxinwei/p/13546326.html