RStudio

packageを使用しないでやってみた(局所線形回帰)

局所線形回帰(Introduction to Nonparametric Regression, Takezawa Kunio, p.49)を実装してみました。Gaussianを重みとした平均2乗誤差を最小化する方法です。各(x,y)値ごとに傾きと切片を求める方法です。

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),x=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82))

a=b=array(0,dim=c(1,nrow(data)))
h=1

for(j in 1:nrow(data)){
 val_a1=(sum(data$x*exp(-((data$x-data$x[j])^2)/(2*h^2)))*sum(data$y*exp(-((data$x-data$x[j])^2)/(2*h^2)))-sum(data$x*data$y*exp(-((data$x-data$x[j])^2)/(2*h^2)))*sum(exp(-((data$x-data$x[j])^2)/(2*h^2))) )
 val_a2=sum(data$x*exp(-((data$x-data$x[j])^2)/(2*h^2)))^2-sum(data$x^2*exp(-((data$x-data$x[j])^2)/(2*h^2)))*sum(exp(-((data$x-data$x[j])^2)/(2*h^2)))

 a[1,j]=val_a1/val_a2

 val_b1=sum(data$x*data$y*exp(-((data$x-data$x[j])^2)/(2*h^2)))*sum(data$x*exp(-((data$x-data$x[j])^2)/(2*h^2)))-sum(data$x^2*exp(-((data$x-data$x[j])^2)/(2*h^2)))*sum(data$y*exp(-((data$x-data$x[j])^2)/(2*h^2)))

val_b2=sum(data$x*exp(-((data$x-data$x[j])^2)/(2*h^2)))^2-sum(exp(-((data$x-data$x[j])^2)/(2*h^2)))*sum(data$x^2*exp(-((data$x-data$x[j])^2)/(2*h^2))) 

b[1,j]=val_b1/val_b2

 }
X=a*data$x+b
result=data.frame(num=1:nrow(data),y=data$y,predict=t(X))

plot(result$num,result$y,xlim=c(1,nrow(data)),ylim=c(min(result$y,result$predict),max(result$y,result$predict)),xlab="sample number",ylab="values",type="o",col=2)
par(new=T)
plot(result$num,result$predict,xlim=c(1,nrow(data)),ylim=c(min(result$y,result$predict),max(result$y,result$predict)),xlab="sample number",ylab="values",type="o",col=3)