Thank you in advance.
Now I want to make comparison of the different bandwidth h in a normal
distribution graph.
This is the table of bandwidth h: thumb rule (normal)--0.00205; thumb
rule(Epanech.)--0.00452; Plug-in (normal)--0.0009;
Plug-in(Epanech.)--0.002.
this is the condition: N=1010 data sample is from normal distribution
N(0,0.0077^2). The grid points are taken to be [-0.05,0.05] and increment is
10. Bandwidth is taken the above h value r respectively and the kernel can
be Epanechnikov kernel or Gaussian kernel.
The following is my code:
#########################################################
# Define the Epanechnikov kernel function
kernel<-function(x){0.75*(1-x^2)*(abs(x)<=1)}
############################################################### # Define the
kernel density estimator
kernden=function(x,z,h,ker){
# parameters: x=variable; h=bandwidth; z=grid point; ker=kernel
nz<-length(z)
nx<-length(x)
x0=rep(1,nx*nz)
dim(x0)=c(nx,nz)
x1=t(x0)
x0=x*x0
x1=z*x1
x0=x0-t(x1)
if(ker==1){x1=kernel(x0/h)} # Epanechnikov kernel
if(ker==0){x1=dnorm(x0/h)} # normal kernel
f1=apply(x1,2,mean)/h
return(f1)
}
#################################################################################################################################
Simulation for different bandiwidths and different kernels
n=1010 # n=1010
ker=1 # ker=1=>Epan; ker=0=> Gaussian
h0=c(0.00452,0.001984) # set initial bandwidths
z=seq(-0.05,0.05,by=10) # grid points
nz=length(z) # number of grid points
x=rnorm(1010, mean=0, sd=0.0077) # simulate x-N(0,0.0077^2)
if(ker==1){h_o=2.34*n^{-0.2}} # bandwidth for Epanechnikov kernel
if(ker==0){h_o=1.06*n^{-0.2}} # bandwidth for normal kernel
f1=kernden(x,z,h0[1],ker)
f2=kernden(x,z,h0[2],ker)
f3=kernden(x,z,h0[3],ker)
f4=kernden(x,z,h0[4],ker)
text1=c("True","h=0.0025","h=0.00452","h=0.0009","h=0.002")
data=cbind(dnorm(z),f1,f2,f3,f4) # combine them as a matrix win.graph()
matplot(z,data,type="l",lty=1:5,col=1:5,xlab="",ylab="")
legend(-1,0.2,text1,lty=1:5,col=1:5)
################################################################
But the error message is "Error in strwidth(legend, units =
"user", cex cex, font = text.font) :
plot.new has not been called yet".
I know something is wrong in the code but don't know where.
Thanks
Regards
--
View this message in context:
http://r.789695.n4.nabble.com/about-different-bandwidths-in-one-graph-tp4636780.html
Sent from the R help mailing list archive at Nabble.com.
On 2012-07-17 11:15, chester123 wrote:> Thank you in advance. > > Now I want to make comparison of the different bandwidth h in a normal > distribution graph. > > This is the table of bandwidth h: thumb rule (normal)--0.00205; thumb > rule(Epanech.)--0.00452; Plug-in (normal)--0.0009; > Plug-in(Epanech.)--0.002. > > this is the condition: N=1010 data sample is from normal distribution > N(0,0.0077^2). The grid points are taken to be [-0.05,0.05] and increment is > 10. Bandwidth is taken the above h value r respectively and the kernel can > be Epanechnikov kernel or Gaussian kernel. > > The following is my code: > ######################################################### > # Define the Epanechnikov kernel function > kernel<-function(x){0.75*(1-x^2)*(abs(x)<=1)} > ############################################################### # Define the > kernel density estimator > kernden=function(x,z,h,ker){ > # parameters: x=variable; h=bandwidth; z=grid point; ker=kernel > nz<-length(z) > nx<-length(x) > x0=rep(1,nx*nz) > dim(x0)=c(nx,nz) > x1=t(x0) > x0=x*x0 > x1=z*x1 > x0=x0-t(x1) > if(ker==1){x1=kernel(x0/h)} # Epanechnikov kernel > if(ker==0){x1=dnorm(x0/h)} # normal kernel > f1=apply(x1,2,mean)/h > return(f1) > } > ################################################################################################################################# > Simulation for different bandiwidths and different kernels > n=1010 # n=1010 > ker=1 # ker=1=>Epan; ker=0=> Gaussian > h0=c(0.00452,0.001984) # set initial bandwidths > z=seq(-0.05,0.05,by=10) # grid points > nz=length(z) # number of grid points > x=rnorm(1010, mean=0, sd=0.0077) # simulate x-N(0,0.0077^2) > if(ker==1){h_o=2.34*n^{-0.2}} # bandwidth for Epanechnikov kernel > if(ker==0){h_o=1.06*n^{-0.2}} # bandwidth for normal kernel > f1=kernden(x,z,h0[1],ker) > f2=kernden(x,z,h0[2],ker) > f3=kernden(x,z,h0[3],ker) > f4=kernden(x,z,h0[4],ker) > text1=c("True","h=0.0025","h=0.00452","h=0.0009","h=0.002") > data=cbind(dnorm(z),f1,f2,f3,f4) # combine them as a matrix win.graph() > matplot(z,data,type="l",lty=1:5,col=1:5,xlab="",ylab="") > legend(-1,0.2,text1,lty=1:5,col=1:5) > ################################################################ > > But the error message is "Error in strwidth(legend, units = "user", cex > cex, font = text.font) : > plot.new has not been called yet". > > I know something is wrong in the code but don't know where.I would guess that you somehow de-activated the graphics device before issuing the legend() call. I can't comment on the rest of your code - haven't examined it. Peter Ehlers
Apparently Analagous Threads
- BreastCancer Dataset for Classification in kknn
- Estimate of baseline hazard in survival
- npudens(np) Error missing value where TRUE/FALSE needed
- density(kernel = "cosine") .. the `wrong cosine' ..
- Treated - KernSmooth pckg - dpik function gives numeric(0) for kernel="epanech"