X = seq(41,60,by=1) #data x
Y = c (196,200,195,207,201,195,202,203,202,198,207,201,196,188,205,199,199,204,204,202) #data y
band= 1
width= 10 #pembobotan yg digunakan 1-10
p=1 #orde polinomial yg digunakan 1
u=1 #lower
v=1 #upperMPL<-function(x,eps=1e-20)
{
x<-as.matrix(x)
xsvd<-svd(x)
diago<-xsvd$d[xsvd$d>eps]
if(length(diago)==1)
{
xplus<-as.matrix(xsvd$v[,1])%*%t(as.matrix(xsvd$u[,1])/diago)
}
else
{
xplus<-
xsvd$v[,1:length(diago)]%*%diag(1/diago)%*%t(xsvd$u[,1:length(diago)])
}
return(xplus)
}
ker<-function(x)
{
kergauss <- (1/sqrt(2 * pi)) * exp((-1/2) * x^2)
return(kergauss)
}gcv_LocalPolynomial = function(X,Y,band,width,p,u,v){
n=length(X)
x0=matrix(seq((min(X)-u),(max(X)+v),by=1),ncol=1)
z=length(x0)
q=c(seq(band,width,by=1))
k=length(q)
M=matrix(nrow=(z*k),ncol=4)
for (j in 1:z){
a=rep(x0[j],n)
X_x0=X-a
h=c(seq(band,width,by=1))
s=length(h)
Gcv<-matrix(nrow=s,ncol=1)
Mse<-matrix(nrow=s,ncol=1)
X_nol=rep(x0[j],s)
for (m in 1:s){
matrix_X=matrix(0,ncol=p+1,nrow=n)
for (i in 1:(p+1))
matrix_X[,i]=X_x0^(i-1)
W=matrix(0,ncol=n,nrow=n)
for (g in 1:n)
W[g,g]=ker(X_x0[g]/h[m])
inverse=MPL(t(matrix_X)%*%W%*%matrix_X)
beta=inverse%*%t(matrix_X)%*%W%*%Y
m_hat=matrix_X%*%beta
H_hat=matrix_X%*%inverse%*%t(matrix_X)%*%W
MSE=(t(Y-m_hat)%*%(Y-m_hat))/n
I=matrix(0,ncol=n,nrow=n)
for (g in 1:n)
I[g,g]=1
GCV=(n^2*MSE)/(sum(diag(I-H_hat)))^2
Gcv[m]=GCV
Mse[m]=MSE
}
R<-matrix(c(X_nol,h,Gcv,Mse),nrow=s)
if (j==1){
M[1:s,]=R}
else {
M[(((j-1)*s)+1):(j*s),]=R}
}
sort.M<-M[order(M[,3]),]
S<-sort.M[1:10,]
cat("Running program untuk order =",p+1,"\n")
cat("=========================================","\n")
cat (" Xo h GCV \n")
cat("=========================================","\n")
print(S)
cat("=========================================","\n")
}gcv_LocalPolynomial(X,Y,band,width,p,u,v)## Running program untuk order = 2
## =========================================
## Xo h GCV
## =========================================
## [,1] [,2] [,3] [,4]
## [1,] 40 10 24.71303 20.01756
## [2,] 59 10 24.71416 20.01847
## [3,] 58 10 24.71468 20.01889
## [4,] 41 10 24.71489 20.01906
## [5,] 41 9 24.71641 20.02029
## [6,] 58 9 24.71758 20.02124
## [7,] 42 9 24.71799 20.02157
## [8,] 57 10 24.71890 20.02231
## [9,] 60 10 24.71896 20.02236
## [10,] 42 10 24.72011 20.02329
## =========================================Berdasarkan hasil output diatas dapat kita perhatikan bahwa x yang optimal adalah x0 = 40 untuk order 2 dengan h = 10 dan nilai dari GCV = 24.71303 dimana ini merupakan nilai GCV yang paling rendah.
LocalPolynomial<-function(p,x0,h){
n=length(X)
X_x0=X-x0
matrix_X=matrix(0,ncol=p+1,nrow=n)
for (i in 1:(p+1))
matrix_X[,i]=X_x0^(i-1)
W=matrix(0,ncol=n,nrow=n)
for (g in 1:n)
W[g,g]=ker(X_x0[g]/h)
inverse=MPL(t(matrix_X)%*%W%*%matrix_X)
beta=inverse%*%t(matrix_X)%*%W%*%Y
m_hat=matrix_X%*%beta
H_hat=matrix_X%*%inverse%*%t(matrix_X)%*%W
MSE=(t(Y-m_hat)%*%(Y-m_hat))/n
I=matrix(0,ncol=n,nrow=n)
for (g in 1:n)
I[g,g]=1
GCV=(n^2*MSE)/(sum(diag(I-H_hat)))^2
cat("=========================================","\n")
cat(" Estimasi Parameter Lokal ","\n")
print(beta)
cat("=========================================","\n")
cat("dengan nilai X0 =",x0,"\n")
cat("Polinomial order =",p+1,"\n")
cat("Bandwidth =",h,"\n")
cat("GCV =",GCV,"\nMSE =",MSE, "\n")
cat("=========================================","\n")
cat(" Estimasi Y")
cat("\n================\n")
print(m_hat)
cat("================\n")
win.graph()
plot(X,Y,type="p",xlim=c(min(X),max(X)),ylim=c(min(Y),max(Y)),xlab=" ",ylab=" ")
par(new=T)
plot(X,m_hat,type="l",col="red",xlim=c(min(X),max(X)),ylim=c(min(Y),max(Y)),xlab="USIA",ylab="KADAR KOLESTEROL")
}LocalPolynomial(1,40,10)## =========================================
## Estimasi Parameter Lokal
## [,1]
## [1,] 199.0792997
## [2,] 0.1117838
## =========================================
## dengan nilai X0 = 40
## Polinomial order = 2
## Bandwidth = 10
## GCV = 24.71303
## MSE = 20.01756
## =========================================
## Estimasi Y
## ================
## [,1]
## [1,] 199.1911
## [2,] 199.3029
## [3,] 199.4147
## [4,] 199.5264
## [5,] 199.6382
## [6,] 199.7500
## [7,] 199.8618
## [8,] 199.9736
## [9,] 200.0854
## [10,] 200.1971
## [11,] 200.3089
## [12,] 200.4207
## [13,] 200.5325
## [14,] 200.6443
## [15,] 200.7561
## [16,] 200.8678
## [17,] 200.9796
## [18,] 201.0914
## [19,] 201.2032
## [20,] 201.3150
## ================Untuk lebih jelas lagi bisa kunjungi refrensi berikut :
YT : https://www.youtube.com/watch?v=3BfL65KXmHs
Ampun…. dah pusing banget coy :)
がんばって ください Renaldi－さん。。。