A Appendix
The output below shows the codes and the function that has been used in this project:
library(CNLTreg)
library(tidyverse)
library(cowplot)
n_t=300
t = seq(from = 0, to = 1 - 1/n_t, by = 1/n_t)
signals<-list("linchirp","mishmash1","heavisine","doppler")
map(signals,function(y) make.signal2(y,x=t)) %>%
map(function(y) ggplot(as.data.frame(y),aes(y=y,x=t))+
geom_line(color="firebrick3")+theme_bw()+
theme(axis.title.x=element_blank(),
axis.title.y=element_blank()))%>%
plot_grid(plotlist = .,
labels = c("A","B","C","D"))x1<-sort(runif(256,-10,10))
y<-3*x1^3+4*x1^2
noise<-rnorm(256,mean=0,sd=sqrt(var(y))/5)
y<-y+noise
est1<-CNLTreg::cnlt.reg(x1,y,P=2,
LocalPred=AdaptPred,neighbours=1,
returnall=FALSE)
estwd <- wavethresh::wd(y)
estwr <- wavethresh:: wr(wavethresh::threshold(estwd,
by.level=TRUE,
policy="universal"))
npbwrh <-np:: npregbw(y ~ x1, regtype = "lc",
bwmethod = "cv.aic")
nprh1 <- np::npreg(npbwrh)
bk <- mgcv:: gam(y~ s(x1,bs="ps",k=60),family=gaussian)
bk<-bk$fitted.values
estfr<-cbind.data.frame(x1,y,
y1=t(est1),
y2=estwr,
y3=fitted(nprh1),
y5=bk)
estfr1<- estfr %>% gather(key="t",value="para",y1,y2,y3,y5)
estfr1 %>% ggplot(aes(x=x1,y=para,color=t,group=t))+geom_line()+theme_light()+theme(axis.title.x=element_blank(),axis.title.y=element_blank())
pr1<-ggplot(estfr,aes(x= x1,y= y))+
geom_point(color="firebrick3")+
geom_line(aes(x=x1,y=y1),color="gray4")+
ggtitle("Second GEN")+theme_light()+ theme(axis.title.x=element_blank(),axis.title.y=element_blank())
pr2<-ggplot(estfr,aes(x= x1,y=y))+
geom_point(color="firebrick3")+
geom_line(aes(x=x1,y=y2),color="gray4")+
ggtitle("First GEN")+theme_light()+ theme(axis.title.x=element_blank(),axis.title.y=element_blank())
pr3<-ggplot(estfr,aes(x= x1,y= y))+
geom_point(color="firebrick3")+
geom_line(aes(x=x1,y=y3),color="gray4")+
ggtitle("Kernel smoother")+theme_light()+ theme(axis.title.x=element_blank(),axis.title.y=element_blank())
pr4<-ggplot(estfr,aes(x= x1,y= y))+
geom_point(color="firebrick3")+
geom_line(aes(x=x1,y=y5),color="gray4")+
theme_light()+ theme(axis.title.x=element_blank(),
axis.title.y=element_blank())+
ggtitle("Spline")
plot_grid(pr1,pr2,pr3,pr4,labels = c("A","B","C","D"))sim_npsiguniffunc <-function(f,distribution,n=256,m=50,r){
# Number of (X, Y) pairs
library(tidyverse)
library(CNLTreg)
library(wavethresh)
n_t = n # number of points on t grid
t = seq(from = 0, to = 1 - 1/n_t, by = 1/n_t)
# Number of simulations
#approx
x=matrix(rep( 0, len=n*m), nrow = n)
y<-matrix(rep( 0, len=n*m), nrow = n)
ft<-f(t)
t_2 = seq(from = 1/(2 * n), by = 1/n, length = n)
ft2 = f(t_2)
#approx
x=matrix( rep( 0, len=n*m), nrow = n)
y=matrix( rep( 0, len=n*m), nrow = n)
ft=f(t)
noise<-rnorm(n,mean=0,sd=sqrt(var(ft))/r)
estwr=matrix( rep( 0, len=n*m), nrow = n)
est1<-matrix( rep( 0, len=n*m), nrow = n)
est2<-matrix( rep( 0, len=n*m), nrow = n)
est3<-matrix( rep( 0, len=n*m), nrow = n)
est4<-matrix( rep( 0, len=n*m), nrow = n)
gridt1<-matrix( rep( 0, len=n*m), nrow = n)
gridy1<-matrix( rep( 0, len=n*m), nrow = n)
gridt2<-matrix( rep( 0, len=n*m), nrow = n)
gridy2<-matrix( rep( 0, len=n*m), nrow = n)
gridt3<-matrix( rep( 0, len=n*m), nrow = n)
gridy3<-matrix( rep( 0, len=n*m), nrow = n)
gridt4<-matrix( rep( 0, len=n*m), nrow = n)
gridy4<-matrix( rep( 0, len=n*m), nrow = n)
mu=matrix( rep( 0, len=n*3), nrow = n)
for (i in 1:m){
x[,i]=sort(distribution)
noise<-rnorm(n,mean=0,sd=sqrt(var(ft))/r)
#generating the y argument
y[,i]=f(x[,i])+noise[i]
#wavelet first generation
estwd=wavethresh::wd(y[,i])
estwr[,i]=t(wavethresh::wr(wavethresh::threshold( estwd,
by.level=TRUE,
policy="universal")))
#second generation
est1[,i]<-CNLTreg:: cnlt.reg(x[,i],y[,i],P=2,
LocalPred=AdaptPred,
neighbours=1,returnall=FALSE)
#kernelsmoother
npbwrh =np:: npregbw(xdat = x[,i],
ydat = y[,i], regtype = "lc",
bwmethod = "cv.aic")
nprh1 = np:: npreg(npbwrh)
est3[,i]=fitted(nprh1)
#Spline-
bk <- mgcv:: gam(y[,i]~ s(x[,i],bs="ps",k=60),
family=gaussian)
est4[,i]<-bk$fitted.values
gridapprox1<-makegrid(x[,i],est1[,i])
gridt1[,i]<-gridapprox1$gridt
gridy1[,i]<-gridapprox1$gridy
gridapprox3<-makegrid(x[,i],est3[,i])
gridt3[,i]<-gridapprox3$gridt
gridy3[,i]<-gridapprox3$gridy
gridapprox4<-makegrid(x[,i],est4[,i])
gridt4[,i]<-gridapprox4$gridt
gridy4[,i]<-gridapprox4$gridy
}
#mse of first generation
bias_firstgen=rowMeans(estwr)-ft
varest_firstgen=apply(estwr,1,var)
mse_firstgen=varest_firstgen+bias_firstgen^2
mae_firstgen=rowMeans(abs(estwr[,1:m]-ft2))
bias_secondgenp2= rowMeans(gridy1) - ft2
varest_secondgenp2=apply(est1,1,var)
mse_secondgenp2=varest_secondgenp2+bias_secondgenp2^2
mae_secondgenp2=rowMeans(abs(gridy1[,1:m]-ft2))
bias_np= rowMeans(gridy3) - ft2
varest_np=apply(est3,1,var)
mse_np=varest_np+bias_np^2
mae_np=rowMeans(abs(gridy3[,1:m]-ft2))
bias_spline= rowMeans(gridy4) - ft2
varest_spline=apply(est4,1,var)
mse_spline=varest_spline+bias_spline^2
mae_spline=rowMeans(abs(gridy4[,1:m]-ft2))
vareste<-cbind(varest_firstgen,varest_secondgenp2,
varest_np,varest_spline)
biase<-cbind(bias_firstgen,bias_secondgenp2,
bias_np,bias_spline)
msee<-cbind(mse_firstgen,mse_secondgenp2,
mse_np,mse_spline)
maee<-cbind(mae_firstgen,mae_secondgenp2,
mae_np,mae_spline)
return_list<-list(vareste,biase,msee,maee)
return(return_list)
}library(ggplot2)
revbetarec<-function(n = 1e5,theta = 0.4,a = 1, b = 2){
library(zipfR)
U=sort(runif(n,0,1))
mix = sample(paste0("component ", 1:2),
size = n, replace = TRUE, prob = c(theta,1 - theta ))
X = ifelse(mix == "component 1", Rbeta.inv(U, a , b ),U)
}
betarec<-function(n = 1e5,theta = 0.4,a = 2, b = 1){
library(zipfR)
U=sort(runif(n,0,1))
mix = sample(paste0("component ", 1:2),
size = n, replace = TRUE, prob = c(theta, 1 - theta))
X = ifelse(mix == "component 1", Rbeta.inv(U, a , b ), U)
}
df<-as.data.frame(betarec())
#a<-revbetarec()
p1<- ggplot(df, aes(x=betarec())) +
geom_density()+theme_light()+
theme(axis.title.x=element_blank(),axis.title.y=element_blank())
df<-as.data.frame(revbetarec())
p2 <- ggplot(df, aes(x=revbetarec())) +
geom_density()+theme_light()+
theme(axis.title.x=element_blank(),axis.title.y=element_blank())
cowplot::plot_grid(p1,p2,labels = c("A","B"))sim_npsig <-function(f,distribution,n=256,m=50,r){
# Number of (X, Y) pairs
library(tidyverse)
library(CNLTreg)
library(wavethresh)
library(pander)
n_t = n # number of points on t grid
t = seq(from = 0, to = 1 - 1/n_t, by = 1/n_t)
x=matrix(rep( 0, len=n*m), nrow = n)
y<-matrix(rep( 0, len=n*m), nrow = n)
ft<-make.signal2(f,x=t)
t_2 = seq(from = 1/(2 * n), by = 1/n, length = n)
ft2 = make.signal2(f,x=t_2)
estwr=matrix( rep( 0, len=n*m), nrow = n)
est1<-matrix( rep( 0, len=n*m), nrow = n)
est2<-matrix( rep( 0, len=n*m), nrow = n)
est3<-matrix( rep( 0, len=n*m), nrow = n)
est4<-matrix( rep( 0, len=n*m), nrow = n)
gridt1<-matrix( rep( 0, len=n*m), nrow = n)
gridy1<-matrix( rep( 0, len=n*m), nrow = n)
gridt2<-matrix( rep( 0, len=n*m), nrow = n)
gridy2<-matrix( rep( 0, len=n*m), nrow = n)
gridt3<-matrix( rep( 0, len=n*m), nrow = n)
gridy3<-matrix( rep( 0, len=n*m), nrow = n)
gridt4<-matrix( rep( 0, len=n*m), nrow = n)
gridy4<-matrix( rep( 0, len=n*m), nrow = n)
mu=matrix( rep( 0, len=n*3), nrow = n)
for (i in 1:m){
x[,i]=sort(distribution)
noise<-rnorm(n,mean=0,sd=sqrt(var(ft))/r)
#generating the y argument
y[,i]=make.signal2(f,x=x[,i])+noise[i]
#wavelet first generation
estwd=wavethresh::wd(y[,i])
estwr[,i]=t(wavethresh::wr(wavethresh::threshold(estwd,
by.level=TRUE,
policy="universal")))
#change mse
est1[,i]<-CNLTreg:: cnlt.reg(x[,i],y[,i],P=2,
LocalPred=AdaptPred,
neighbours=1,
returnall=FALSE)
npbwrh =np:: npregbw(xdat = x[,i],ydat = y[,i],
regtype = "lc",
bwmethod = "cv.aic")
nprh1 = np:: npreg(npbwrh)
est3[,i]=fitted(nprh1)
bk <- mgcv:: gam(y[,i]~ s(x[,i],bs="ps",k=60),
family=gaussian)
est4[,i]<-bk$fitted.values
gridapprox1<-makegrid(x[,i],est1[,i])
gridt1[,i]<-gridapprox1$gridt
gridy1[,i]<-gridapprox1$gridy
gridapprox3<-makegrid(x[,i],est3[,i])
gridt3[,i]<-gridapprox3$gridt
gridy3[,i]<-gridapprox3$gridy
gridapprox4<-makegrid(x[,i],est4[,i])
gridt4[,i]<-gridapprox4$gridt
gridy4[,i]<-gridapprox4$gridy
}
#mse of first generation
bias_firstgen=rowMeans(estwr)-ft
varest_firstgen=apply(estwr,1,var)
mse_firstgen=varest_firstgen+bias_firstgen^2
mae_firstgen=rowMeans(abs(estwr[,1:m]-ft2))
bias_secondgenp2= rowMeans(gridy1) - ft2
varest_secondgenp2=apply(est1,1,var)
mse_secondgenp2=varest_secondgenp2+bias_secondgenp2^2
mae_secondgenp2=rowMeans(abs(gridy1[,1:m]-ft2))
bias_np= rowMeans(gridy3) - ft2
varest_np=apply(est3,1,var)
mse_np=varest_np+bias_np^2
mae_np=rowMeans(abs(gridy3[,1:m]-ft2))
bias_spline= rowMeans(gridy4) - ft2
varest_spline=apply(est4,1,var)
mse_spline=varest_spline+bias_spline^2
mae_spline=rowMeans(abs(gridy4[,1:m]-ft2))
vareste<-cbind(varest_firstgen,varest_secondgenp2,
varest_np,varest_spline)
biase<-cbind(bias_firstgen,bias_secondgenp2,
bias_np,bias_spline)
msee<-cbind(mse_firstgen,mse_secondgenp2,
mse_np,mse_spline)
maee<-cbind(mae_firstgen,mae_secondgenp2,
mae_np,mae_spline)
return_list<-list(vareste,biase,msee,maee)
return(return_list)
}msesum<-function(x){
apply(x,2,sum)
}
library(tidyverse)
library(pander)
library(kableExtra)
library(knitr)
printmse<-function(data_ms,m,cap,function_name=c("linchirp","mishmash","heavisine","doppler"),estimator_name=c("First GEN","Second GEN","Kernel smoother","Spline"),n=4){
df<-map_dfc(c(1:n),function(i)map(data_ms[[i]][m],msesum)%>%unlist(use.names=TRUE)%>%tibble())%>%setNames(function_name)%>%add_column(estimators=estimator_name,.before = 1)#make the table
emphasize.strong.cells(cbind( apply(df[,2:(n+1)], 2, which.min),2:(n+1)))
pander::pander(df,caption=cap, style = "rmarkdown")
}