应用回归分析课作业摸索

应用回归分析课作业摸索

写在前面

其实还没看懂你什么意思，只能先摸索一下，学点别人的东西。

作业描述

写一篇论文，通过比较ols和稳健回归（自己找1-2种稳健方法）的结果，得出稳健回归更好的结论。生成一个表格：
 

• case1：误差项服从正态分布
• case2：误差项服从柯西分布
• 预期结果：
  • 1<2，1<3：就是ols更好才对
  • 二<一，三<一：用柯西则相反
• 具体怎么做（数值模拟）：
  • 先通过y=β0+β1x+E，自己设定β0和β1，生成数据点（xi，yi）【i大概做50/200个左右】
  • 用（xi，yi）去估β0，β1，看看哪种方法估得好（搞个图之类的）
  • 还得重复1000次，取对β0和β1的平均
• 步骤：
  • 第一步先随机生成数据点

# 假设现在β0=0.5，β1=2，生成100个点set.seed(20)
x<-rnorm(100)
e<-rnorm(100)
y<-0.5+2*x+esummary(y)
plot(x,y)

• 第二步：搞个循环，分别去做ols和稳健回归

• 第三步：生成表格，把东西放在一起

感觉找到能用的代码：

k = 100000 # 定义实验次数
beta_x1 = c() # 定义空列
beta_x2 = c()
for (i in 1:k) {beta1 = 0.8 beta0 = 1.5 # 设置真实系数x = rnorm(100,4,2) # 产生随机数error = rnorm(100,0,1) # 产生随机误差y = beta0+beta1*x+errordata = data.frame(x,y) # 构建数据框ols = lm(y~x1+x2,data = data)ols = summary(ols)beta1_r = (ols$coefficients[2]-beta1)^2 # 取估计的系数beta2_r = (ols$coefficients[3]-beta2)^2beta_x1 = c(beta_x1,beta1_r) # 追加成向量beta_x2 = c(beta_x2,beta2_r)}
MSEx1 = sum(beta_x1)/length(beta_x1) # 求MES的公式
MSEx2 = sum(beta_x2)/length(beta_x2)
message('x1的MES为：',MSEx1)
message('x2的MES为：',MSEx2)

0917的进度：好像是已经把函数弄出来了

但是其实是2*3的条件（误差项可以服从两种不同的分布，回归会用三种不同的方法）

在考虑怎么写循环更方便一点

参考的链接

R语言进行数值模拟：模拟泊松回归模型的数据 – 拓端 (tecdat)

感觉这个步骤挺像的，不过泊松回归是什么...?

开始学习别人的代码

验证回归模型的首选方法是模拟来自它们的数据，并查看模拟数据是否捕获原始数据的相关特征。

感兴趣的基本特征是平均值（和作业要求的好像一样）

标准回归模型假设存在将预测变量与结果相关联的真实/固定参数。但是，当我们执行回归时，我们只估计这些参数。因此，回归软件返回表示系数不确定性的标准误差。

# 别人用的代码（注：下面的“我”是原作者本人）# 我模拟了两个预测变量,使用50的小样本。
n<-50 
set.seed(18050518)# xc的系数为0.5 ，xb的系数为1 。我对预测进行取幂，并使用该rpois()函数生成泊松分布结果#指数预测并传递给rpois（）【这里是不是没有步骤】
summary(dat)glm(formula = y ~ xc + xb, family = poisson, data = dat)# 估计的系数与人口模型相距不太远，.21代表截距而不是0，.46而不是.5，而0.81而不是1。# 接下来模拟模型中的数据，我想要10,000个模拟数据集，为了捕捉回归系数的不确定性，我假设系数来自多元正态分布，估计系数作为均值，回归系数的方差 – 协方差矩阵作为多元正态分布的方差 – 协方差矩阵。coefs <- mvrnorm(n = 10000, mu = coefficients(fit.p), Sigma = vcov(fit.p))# 检查模拟系数与原始系数的匹配程度。
coefficients(fit.p)
colMeans(coefs) #模拟系数的平均值# 标准误
sqrt(diag(vcov(fit.p)))
apply(coefs, 2, sd) #模拟系数的标准偏差# 下一步是模拟模型中的数据。我们通过将模拟系数的每一行乘以原始预测变量来实现。然后我们传递预测：
#每种情况一行，每组模拟系数一行#模型矩阵与系数的乘积，取幂，
#然后用于模拟泊松分布的结果for (i in 1:nrow(coefs)) {sim.dat[, i] <- rpois(n, exp(fit.p.mat %*% coefs[i ]))
}
rm(i, fit.p.mat)# 现在一个是完成模拟，将模拟数据集与原始数据集至少比较结果的均值和方差：
c(mean(dat$y), var(dat$y))#原始结果的均值和方差
c(mean(colMeans(sim.dat)), mean(apply(sim.dat, 2, var)))#10,000个模拟结果的平均值和变量的平均值# 模拟结果的平均值略高于原始数据，平均方差更高。平均而言，可以预期方差比平均值更偏离目标。方差也将与一些极高的值正偏差，同时，它的界限为零，因此中位数可能更好地反映了数据的中心：
# 中位数方差更接近原始结果的方差。# 得到模拟均值与方差的分布
### 绘制10,000个模拟数据集中的一些数据集的直方图并将其与原始结果的直方图进行比较也是有用的。还可以测试原始数据和模拟数据集中xb = 1和xb = 0之间结果的平均差异。回到基础R，它具有simulate()执行相同操作的功能：
sim.default <- simulate(fit.p, 10000)
此代码相当于：
sim.default <- replicate(10000, rpois(n, fitted(fit.p)))
fitted(fit.p)是响应尺度的预测，或指数线性预测器，因为这是泊松回归。因此，我们将使用模型中的单组预测值来重复创建模拟结果。c(mean(colMeans(sim.default)), mean(apply(sim.default, 2, var)),[1] 2.020036 3.931580 3.810612
与忽略系数不确定性时相比，均值和方差更接近原始结果的均值和方差。