统计理论与方法第二讲的内容，火山GG布置了作业，因为要迭代，如果用excel的话要操作到死，何况题目里去掉一个样本就要再做一次，如果一步一步人工操作的话那工作量真是无法想象。恰逢要学R，就凭自己的入门级水平努力了几天，凑出了这个粗糙的函数。有了这个函数，有其他数据要用一样方法拟合直接扔进去就好了哇哈哈。虽然粗糙，不过怎么说也是自己第一个正式函数，值得纪念吖...
说说过程，先是按照老师给的方法一步一步的按照给定数据写下来，然后run通过；第二步是将迭代部分用循环语句实现；第三步是编函数，实现任意数据该算法拟合。
说说经验，程序里出错最多的几个地方：1."："的左右要注意，很多时候要括号，因为冒号的优先级高；2.下意识的错误--2k、3N这样的错误，要记得表示成2*k、3*N...3.明确每个对象的类型，向量就不要用[,i]的索引，数据框就不要x[i]；4.赋值是变量的要特别小心，由其是分组不确定组长时。
下面就是成果啦：
说说过程，先是按照老师给的方法一步一步的按照给定数据写下来，然后run通过；第二步是将迭代部分用循环语句实现；第三步是编函数，实现任意数据该算法拟合。
说说经验，程序里出错最多的几个地方：1."："的左右要注意，很多时候要括号，因为冒号的优先级高；2.下意识的错误--2k、3N这样的错误，要记得表示成2*k、3*N...3.明确每个对象的类型，向量就不要用[,i]的索引，数据框就不要x[i]；4.赋值是变量的要特别小心，由其是分组不确定组长时。
下面就是成果啦：
Res<-function(x,precision1,precision2){
if(dim(x)[2]==3){x1<-x[,1];x<-x[,c(2,3)]} ##把输入的数据框一般化并按自变量大小且升序排列
o<-order(x[,1])
x<-x[o,]
k<-dim(x)[1]%/%3 ##确定三组的组长
remainder<-dim(x)[1]%%3
if(remainder==0){L<-x[1:k,];M<-x[(k+1):(2*k),];R<-x[(2*k+1):(3*k),];}
if(remainder==1){L<-x[1:k,];M<-x[(k+1):(2*k+1),];R<-x[(2*k+2):(3*k+1),]}
if(remainder==2){L<-x[1:(k+1),];M<-x[(k+2):(2*k+1),];R<-x[(2*k+2):(3*k+2),]}
amountL<-dim(L)[1]
amountM<-amountL+dim(M)[1]
amountR<-amountM+dim(R)[1]
XL<-summary(L[,1])[["Median"]] ##开始拟合
XM<-summary(M[,1])[["Median"]]
XR<-summary(R[,1])[["Median"]]
YL<-summary(L[,2])[["Median"]]
YM<-summary(M[,2])[["Median"]]
YR<-summary(R[,2])[["Median"]]
b<-numeric();a<-numeric();c<-numeric();d<-numeric()
b[1]<-(YR-YL)/(XR-XL)
a[1]<-1/3*((YL-b[1]*XL)+(YM-b[1]*XM)+(YR-b[1]*XR))
kexi<-(x[,2]-(a[1]+b[1]*x[,1])) ##求残差
savekexi<-kexi
x<-cbind(x,kexi)
L<-x[1:amountL,] ##下面把残差对自变量回归
M<-x[(amountL+1):amountM,]
R<-x[(amountM+1):amountR,]
XL<-summary(L[,1])[["Median"]]
XM<-summary(M[,1])[["Median"]]
XR<-summary(R[,1])[["Median"]]
kexiL<-summary(L$kexi)[["Median"]]
kexiM<-summary(M$kexi)[["Median"]]
kexiR<-summary(R$kexi)[["Median"]]
d[1]<-(kexiR-kexiL)/(XR-XL)
c[1]<-1/3*((kexiL-d[1]*XL)+(kexiM-d[1]*XM)+(kexiR-d[1]*XR))
i=1 ##迭代
while((abs(c[i])>precision1)|abs((d[i])>precision2)){
b[i+1]<-b[i]+d[i]
a[i+1]<-a[i]+c[i]
kexi<-(x[,2]-(a[i+1]+b[i+1]*x[,1]))
x$kexi<-kexi
savekexi<-cbind(savekexi,kexi)
L<-x[1:amountL,]
M<-x[(amountL+1):amountM,]
R<-x[(amountM+1):amountR,]
kexiL<-summary(L$kexi)[["Median"]]
kexiM<-summary(M$kexi)[["Median"]]
kexiR<-summary(R$kexi)[["Median"]]
i=i+1
d[i]<-(kexiR-kexiL)/(XR-XL)
c[i]<-1/3*((kexiL-d[i]*XL)+(kexiM-d[i]*XM)+(kexiR-d[i]*XR))
}
No.kexi<-paste("kexi",1:i) #整理数据、合并、显示
colnames(savekexi)<-No.kexi
All<-cbind(x1,x[,c(1,2)],savekexi)
argument<-cbind(a,b,c,d)
list(All, argument)
}
发表评论