自编函数

其他文章用的有自己编的，一起放到这里

#######################自己编的函数，运行后直接调用#####################
#计算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)
}

#学生化残差
ri = function(ei,X)
{
  p = ncol(X) #回归系数个数
  mse = mse(ei,p) #残差均方
  hii = H(X)  #返回帽子矩阵对角线元素
  ans = ei/sqrt(mse*(1-hii)) #学生化残差
  return(ans)
}

#外部学生化残差
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,xlab = "学生化残差",ylab = "百分比")  #画正态概率图
  #添加回归线
  fm = lm(Pi~ti)
  abline(fm)
}

#进行失拟检验
# 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)
}