packageを使わないでやってみた(回帰分析におけるブートストラップ)

標本データの中からいくつかの決められた数のデータを複数回抽出することにより、標本データに偏りがないかを調査します。今回は線形回帰におけるパラメータの誤差がデータの抽出の仕方によりどれくらい変化が起こるかを調べるプログラムを書きました。

```{r}

#bootstrap 


library(dplyr)
data=data.frame(num=1:20,y=c(167,167.5,168.4,172,155.3,151.4,163,174,168,160.4,164.7,171,162.6,164.8,163.3,167.6,169.2,168,167.4,172),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56))

times=1000

X=cbind(data$x1,data$x2)

Y=data$y

a=solve(t(X)%*%X)%*%t(X)%*%Y;b=mean(Y)-sum(a*apply(X,2,mean))

params=array(0,dim=c(times,(ncol(X)+1)))

num=10

AIC=c();BIC=c()

nums_matrix=array(0,dim=c(times,nrow(X)))

for(j in 1:times){
  
nums=sample(c(1:nrow(X)),num)

nums_matrix[j,]=nums

X_sub=array(0,dim=c(1,ncol(X)))

Y_sub=c()

for(i in 1:length(nums)){
  
X_sub=rbind(X_sub,X[nums[i],])  

Y_sub=c(Y_sub,Y[nums[i]])
  
}
  
X_sub=tail(X_sub,nrow(X_sub)-1)

a_sub=solve(t(X_sub)%*%X_sub)%*%t(X_sub)%*%Y_sub;

b_sub=mean(Y_sub)-sum(a_sub*apply(X_sub,2,mean));

sigma_hat=(1/num)*sum((Y_sub-(X_sub%*%a_sub+b_sub))^2)

AIC=c(AIC,num*log(2*pi*sigma_hat)+num+2*(ncol(X)+2))

BIC=c(BIC,num*log(2*pi*sigma_hat)+num+ncol(X)*log(num))

params[j,]=c(a_sub,b_sub)
  
}


#b=apply(params,2,mean)-c(a,b)

sigma=array(0,dim=c(1,length(c(a,b))))

for(j in 1:length(sigma)){
  
sigma[1,j]=(1/(times-1))*sum((params[,j]-apply(params,2,mean)[j])^2)  
  
}

w=c()

for(j in 1:length(AIC)){
  
diff_AIC=AIC[j]-min(AIC)

all_AIC=(AIC-min(AIC))
  
w=c(w,exp(-(1/2)*(diff_AIC))/sum(exp(-(1/2)*all_AIC)))  
  
}

parameters_AIC=t(params)%*%w/sum(w)

w=c()

for(j in 1:length(BIC)){
  
diff_BIC=BIC[j]-min(BIC)

all_BIC=(BIC-min(BIC))
  
w=c(w,exp(-(1/2)*(diff_BIC))/sum(exp(-(1/2)*all_BIC)))  
  
}

parameters_BIC=t(params)%*%w/sum(w)

#parameters[length(parameters)]=mean(params[,length(parameters)])


```

```{r}

#params;G(x)

sep=0.0001;k=1;

inv_G_data=data.frame(alpha=seq(0,1,sep),inv_G=0)

G_data=data.frame(x=seq(-0.1,0.1,sep),G=0)

x=(1/times)*(params[,k]-c(a,b)[k])

for(j in 1:nrow(G_data)){
  
G_data$G[j]=length(x[x<=G_data$x[j]])/length(x)  

}

x=data.frame(num=1:length(x),x=x);

x=x[order(x$x),];x=x %>% dplyr::select(-num);x=data.frame(num=1:nrow(x),x)

for(j in 1:nrow(inv_G_data)){
  
if(floor(times*inv_G_data$alpha[j])>1){
  
inv_G_data$inv_G[j]=as.numeric(x$x[x$num==floor(times*inv_G_data$alpha[j])] )

}

}

inv_G_data=inv_G_data %>% dplyr::filter(inv_G!=0)

G_data=G_data %>% dplyr::filter(G!=0) %>% dplyr::filter(G!=1)


plot(G_data$x,G_data$G,col=2,type="o",main=paste0("num:",num),xlab="x",ylab="G")

plot(inv_G_data$inv_G,inv_G_data$alpha,col=3,type="o",main="inv_G",xlab="quantiles",ylab="alpha")

```

```{r}

#モデル平均化法


predict1=X%*%parameters_AIC[1:2]+parameters_AIC[3]

predict2=X%*%parameters_BIC[1:2]+parameters_BIC[3]

result=data.frame(n=1:length(predict),Y=Y,predict1=predict1,predict2=predict2,predict3=(X%*%a+b))

plot(result$n,result$Y,col=2,type="o",xlab="n",ylab="values",ylim=c(min(Y,result$predict1),max(Y,result$predict1)),main=paste0("model_ave_AIC_cor:",cor(Y,result$predict1)))

par(new=T)

plot(result$n,result$predict1,col=3,type="o",xlab="n",ylab="values",ylim=c(min(Y,result$predict1),max(Y,result$predict1)),main=paste0("model_ave_AIC_cor:",cor(Y,result$predict1)))

plot(result$n,result$Y,col=2,type="o",xlab="n",ylab="values",ylim=c(min(Y,result$predict2),max(Y,result$predict2)),main=paste0("model_ave_BIC_cor:",cor(Y,result$predict2)))

par(new=T)

plot(result$n,result$predict2,col=3,type="o",xlab="n",ylab="values",ylim=c(min(Y,result$predict2),max(Y,result$predict2)),main=paste0("model_ave_BIC_cor:",cor(Y,result$predict2)))

plot(result$n,result$Y,col=2,type="o",xlab="n",ylab="values",ylim=c(min(Y,result$predict3),max(Y,result$predict3)),main=paste0("multiple_reg_cor:",cor(Y,result$predict3)))

par(new=T)

plot(result$n,result$predict3,col=3,type="o",xlab="n",ylab="values",ylim=c(min(Y,result$predict3),max(Y,result$predict3)),main=paste0("multiple_reg_cor:",cor(Y,result$predict3)))

#モデル平均化法AIC
(1/nrow(result))*sum((Y-predict1)^2)

#モデル平均化法BIC
(1/nrow(result))*sum((Y-predict2)^2)

#線形回帰
(1/nrow(result))*sum((Y-(X%*%a+b))^2)

```