kernel smoothing


```{r}

#kernel smoothing(nadaraya-watson estimator)

library(dplyr)

kernel=function(x){
  
  return(exp(-(x)^2/(2*3^2)))

}

h_n=0.1

val1=-3;val2=8;x=seq(val1,val2,h_n)

kernel_data=data.frame(x=seq(val1,val2,h_n),K_h_n=0,F_h_n=0,m_h=0,y=cos(seq(val1,val2,h_n)^2))

n=nrow(kernel_data)

for(j in 1:nrow(kernel_data)){
  
kernel_data$K_h_n[j]=kernel((kernel_data$x[j])/h_n)/h_n  

kernel_data$F_h_n[j]=sum(kernel((kernel_data$x[j]-x)/h_n))/(h_n*n)

kernel_data$m_h[j]=(sum(kernel((kernel_data$x[j]-x)/h_n)*kernel_data$y)/(h_n))/(sum(kernel((kernel_data$x[j]-x)/h_n))/(h_n))

}

sigma_vec=c()

for(j in 1:nrow(kernel_data)){
  
value=sum((kernel((kernel_data$x[j]-x)/h_n)/h_n)*(kernel_data$y-kernel_data$m_h[kernel_data$x==kernel_data$x[j]])^2)/(sum(kernel((kernel_data$x[j]-x)/h_n))/(h_n))

value=value/n

sigma_vec=c(sigma_vec,value)
  
}

#c_K=1

CUP=c();CLO=c()

alpha=0.05;c_k=2

for(j in 1:nrow(kernel_data)){

  CUP=c(CUP,kernel_data$m_h[j]+qnorm(1-alpha)*c_k*sigma_vec[j]/sqrt(n*h_n*kernel_data$F_h_n[j]))  
  
  CLO=c(CLO,kernel_data$m_h[j]-qnorm(1-alpha)*c_k*sigma_vec[j]/sqrt(n*h_n*kernel_data$F_h_n[j]))
  
}

kernel_data=kernel_data %>% mutate(CUP=CUP,CLO=CLO)

#predict values and test values

plot(kernel_data$x,kernel_data$y,type="o",col=2,xlim=c(val1,val2),ylim=c(min(c(kernel_data$y,kernel_data$m_h)),max(c(kernel_data$y,kernel_data$m_h))),xlab="x",ylab="y",main="kernel smoothing")

par(new=T)

plot(kernel_data$x,kernel_data$m_h,type="o",col=3,xlim=c(val1,val2),ylim=c(min(c(kernel_data$y,kernel_data$m_h)),max(c(kernel_data$y,kernel_data$m_h))),xlab="x",ylab="y",main="kernel smoothing")

```

```{r}
#predict values and confidence interval

plot(kernel_data$x,kernel_data$CUP,type="o",col=3,xlim=c(val1,val2),ylim=c(min(c(kernel_data$y,kernel_data$m_h)),max(c(kernel_data$y,kernel_data$m_h))),xlab="x",ylab="y",main="kernel smoothing")

par(new=T)

plot(kernel_data$x,kernel_data$m_h,type="o",col=2,xlim=c(val1,val2),ylim=c(min(c(kernel_data$y,kernel_data$m_h)),max(c(kernel_data$y,kernel_data$m_h))),xlab="x",ylab="y",main="kernel smoothing")

par(new=T)

plot(kernel_data$x,kernel_data$CLO,type="o",col=3,xlim=c(val1,val2),ylim=c(min(c(kernel_data$y,kernel_data$m_h)),max(c(kernel_data$y,kernel_data$m_h))),xlab="x",ylab="y",main="kernel smoothing")


```
```{r}

#semi-parametric reg model

library(dplyr)


n=50

Y=sample(80:100,n,replace=T)

z1=n;z2=10

Z=array(0,dim=c(z1,z2))

for(j in 1:ncol(Z)){
  
Z[,j]=sample(1:20,n,replace = T)  
  
}

w1=n;w2=15

W=array(0,dim=c(n,w2))

for(j in 1:w2){
  
W[,j]=sample(30:50,n,replace = T)  
  
}


P_W=W%*%solve(t(W)%*%W)%*%t(W)

spline=0

if(spline==1){
  
  spline_func=function(x){
    
  mat=array(0,dim=c(n,n))
  
  for(i in 1:n){
  for(j in 1:n){
    
    
mat[i,j]=1-exp(-mean(abs(x[i,]-x[j,])))
    
  }  
  }
    
  return(mat/c(apply(mat,1,mean)) ) 
  }
  
  W_h=spline_func(Z)
  
  Z_bar=(diag(1,n)-W_h)%*%Z
  
  Y_bar=(diag(1,n)-W_h)%*%Y
  
  b=solve(t(Z_bar)%*%Z_bar)%*%t(Z_bar)%*%Y_bar
  
  m_h=W_h%*%(Y-Z%*%b)
  
  predict=Z%*%b+m_h
  
}else{
  
  b=solve(t(Z)%*%t(diag(1,n)-P_W)%*%(diag(1,n)-P_W)%*%Z)%*%t(Z)%*%t(diag(1,n)-P_W)%*%(diag(1,n)-P_W)%*%Y

  r=solve(t(W)%*%W)%*%t(W)%*%P_W%*%(Y-Z%*%b)
  
  predict=Z%*%b+W%*%r

}



#m_h=W%*%(Y-Z%*%b)

plot(c(1:n),Y,xlim=c(1,n),ylim=c(min(c(Y,predict)),max(c(Y,predict))),col=2,type="o",ylab="values")

par(new=T)

plot(c(1:n),predict,xlim=c(1,n),ylim=c(min(c(Y,predict)),max(c(Y,predict))),col=3,type="o",ylab="values")

cor(Y,predict)

```

```{r}

n=50

Y=sample(80:100,n,replace=T)

z1=n;z2=10

Z=array(0,dim=c(z1,z2))

for(j in 1:ncol(Z)){
  
Z[,j]=sample(1:200,n,replace = T)  
  
}

kernel=function(x){
  
  val=c()
  
  for(i in 1:n){
    
  val=c(val,1-1/(1+(mean(abs(x-Z[i,])))))  
    
  }
  
  return(val)
  
}

W_fun=function(x){
  
return(apply(x,1,sum)/sum(x))  
  
}


m_h=c()

for(i in 1:n){
  
 m_h=c(m_h,mean(kernel(Z[i,])*Y)) 
  
}

W_X=W_fun(Z)

p_h=mean(W_X*(Y-m_h)^2)

K0=1

G=p_h*exp(2/n)

```
```{r}

#kernel reg p.55


library(dplyr)

n=100;k=7

Z=array(0,dim=c(n,k))

Y=array(0,dim=c(n,k))

X=c()


for(j in 1:n){
  
val=rpois(k,sample(10:14,1))

Z[j,]=val/sum(val)

val=rpois(k,sample(19:31,1))

Y[j,]=val

X=c(X,rpois(1,sample(1:10,1)))
  
}

#Z=apply(Z,2,mean)

h=0.1

kernel=function(x,y){
  
return(exp(-(x-y)^2/(2*h^2))/(2*pi))  
  
}

h_s=function(x,s){
  
  return(sum(kernel(x,X)*Z[,s])/sum(kernel(x,X)))
  
  }

E_y_x=array(0,dim=c(n,k))

#V_E_y=array(0,dim=c(n,k))

for(j in 1:k){
  
val1=c();val2=c();val3=c();val4=c();val5=c()

for(i in 1:n){
  
val1=c(val1,sum(Z[i,j]*kernel(X,X[i])*Y[i,j]))  
  
val2=c(val2,sum(Z[i,j]*kernel(X,X[i])))

#val3=c(val3,sum(Z[i,j]*kernel(X,X[i])*Y[i,j]^2))

#val4=c(val4,sum(Z[i,j]*kernel(X,X[i]))*h_s(X[i],j))

#val5=c(val5,(h_s(X[i],j)-1)/h_s(X[i],j))

}

E_y_x[,j]=val1/val2;

#V_E_y[,j]=val3/val4+val5*(val1/val2)^2 

}




E_y=c()

for(j in 1:k){
  
E_y=c(E_y,mean(Z[,j]*Y[,j]+(1-Z[,j])*E_y_x[,j])) 
  
}

kernel_mat=array(0,dim=c(n,n))

for(j in 1:n){
  
kernel_mat[j,]=kernel(X,X[j])  
  
}



h_mat=array(0,dim=c(n,k))

for(j in 1:k){
for(l in 1:n){
  
h_mat[l,j]=h_s(X[l],j)  
  
}  
}

cov_E_y_mat=array(0,dim=c(k,k))

for(j in 1:k){
for(i in 1:k){  
  
cov_E_y_mat[i,j]=mean(E_y_x[,i]*E_y_x[,j])-E_y[i]*E_y[j] 
   
}
}

```

```{r}


#mixed effects model p.180

library(dplyr);library(mvtnorm)


n=500

Y=sample(80:100,n,replace=T)

z1=n;z2=10

Z=array(0,dim=c(z1,z2))

for(j in 1:ncol(Z)){

Z[,j]=sample(1:20,n,replace = T)  

}

w1=n;w2=15

W=array(0,dim=c(n,w2))

for(j in 1:w2){

W[,j]=sample(30:50,n,replace = T)  

}


beta=solve(t(Z)%*%Z)%*%t(Z)%*%Y

sigma=array(1,dim=c(n,n))+diag(1,n)

pthi=array(2,dim=c(w2,w2))+diag(1,w2)

ite=10

for(j in 1:ite){
  
y_adj=Y-Z%*%beta

u=solve(t(W)%*%solve(sigma)%*%W+solve(pthi))%*%t(W)%*%solve(sigma)%*%y_adj
  
y_adj=Y-W%*%u

beta=solve(t(Z)%*%solve(sigma)%*%Z)%*%t(Z)%*%solve(sigma)%*%y_adj  

predict=Z%*%beta+W%*%u

sigma=cov(rmvnorm(1000000,c(predict),sigma))
 
pthi=cov(rmvnorm(1000000,rep(0,w2),pthi))

print(sum(abs(predict-Y)))
 
}

```