R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.1.1  (2005-06-20), ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for a HTML browser interface to help.
Type 'q()' to quit R.

> ### * <HEADER>
> ###
> attach(NULL, name = "CheckExEnv")
> assign(".CheckExEnv", as.environment(2), pos = length(search())) # base
> ## add some hooks to label plot pages for base and grid graphics
> setHook("plot.new", ".newplot.hook")
> setHook("persp", ".newplot.hook")
> setHook("grid.newpage", ".gridplot.hook")
> 
> assign("cleanEx",
+        function(env = .GlobalEnv) {
+ 	   rm(list = ls(envir = env, all.names = TRUE), envir = env)
+            RNGkind("default", "default")
+ 	   set.seed(1)
+    	   options(warn = 1)
+ 	   delayedAssign("T", stop("T used instead of TRUE"),
+ 		  assign.env = .CheckExEnv)
+ 	   delayedAssign("F", stop("F used instead of FALSE"),
+ 		  assign.env = .CheckExEnv)
+ 	   sch <- search()
+ 	   newitems <- sch[! sch %in% .oldSearch]
+ 	   for(item in rev(newitems))
+                eval(substitute(detach(item), list(item=item)))
+ 	   missitems <- .oldSearch[! .oldSearch %in% sch]
+ 	   if(length(missitems))
+ 	       warning("items ", paste(missitems, collapse=", "),
+ 		       " have been removed from the search path")
+        },
+        env = .CheckExEnv)
> assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now
> assign("ptime", proc.time(), env = .CheckExEnv)
> grDevices::postscript("aws-Examples.ps")
> assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv)
> options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly"))
> options(warn = 1)    
> library('aws')
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "aws"
> 
> ### * aws
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: aws
> ### Title: Local polynomial Adaptive Weights Smoothing for regression with
> ###   additive errors
> ### Aliases: aws
> ### Keywords: smooth nonparametric regression models
> 
> ### ** Examples
> 
> #######################################################
> #                                                     
> #                  univariate examples                
> #                                                     
> ####################################################### 
> #  
> #    Blocks data (Example 6 from Fan & Gijbels (1996)
> #
>      mofx6 <- function(x){
+      xj <- c(10,13,15,23,25,40,44,65,76,78,81)/100
+      hj <- c(40,-50,30,-40,50,-42,21,43,-31,21,-42)*.37
+      Kern <- function(x) (1-sign(x))/2
+      apply(Kern(outer(xj,x,"-"))*hj,2,sum)
+      }
> #
> #    sigma==1   step by step with graphics
> #
>      fx6 <- mofx6(seq(0,1,1/2047))
>      y <- rnorm(fx6,fx6,1)
>      tmp <- aws(y,p=0,hmax=100,graph=TRUE)
>      par(mfrow=c(1,1),mar=c(3,3,2.5,.5),mgp=c(2,1,0))
>      plot(seq(0,1,1/2047),y)
>      lines(seq(0,1,1/2047),tmp$theta,col=3,lwd=2)
>      lines(seq(0,1,1/2047),fx6,col=2)
> #
> #    sigma==3   without graphics (much faster)
> #
>      y <- rnorm(fx6,fx6,3)
>      tmp <- aws(y,hmax=500)
>      par(mfrow=c(1,1),mar=c(3,3,2.5,.5),mgp=c(2,1,0))
>      plot(seq(0,1,1/2047),y)
>      lines(seq(0,1,1/2047),tmp$theta,col=3,lwd=2)
>      lines(seq(0,1,1/2047),fx6,col=2)
>      rm(mofx6,fx6,y,tmp)
> #  
> #    second example from Polzehl and Spokoiny (2002)
> #
>      f2 <- function(x) sin(2*pi*1.2/(x+.2))
>      n <- 1000
>      x <- seq(0,1,length=n)
>      fx2 <- f2(x)
>      set.seed(1)
>      sigma <- .25
>      y <- rnorm(x,fx2,sigma)
> #   increase hmax to 2 for good results
>      ex1p0 <- aws(y,x,p=0,hmax=.1)$theta[1,]
>      ex1p1 <- aws(y,x,p=1,hmax=.1)$theta[1,]
>      ex1p2 <- aws(y,x,p=2,hmax=.1)$theta[1,]
>      ex1p3 <- aws(y,x,p=3,hmax=.1)$theta[1,]
>      par(mfrow=c(2,2),mar=c(2.5,2.5,2.5,.5),mgp=c(1.5,.5,0))
>      plot(x,y)
>      lines(x,fx2,col=2)
>      lines(x,ex1p0,col=3,lwd=2)
>      title("local constant AWS")
>      plot(x,y)
>      lines(x,fx2,col=2)
>      lines(x,ex1p1,col=3,lwd=2)
>      title("local linear AWS")
>      plot(x,y)
>      lines(x,fx2,col=2)
>      lines(x,ex1p2,col=3,lwd=2)
>      title("local quadratic AWS")
>      plot(x,y)
>      lines(x,fx2,col=2)
>      title("local cubic AWS")
>      lines(x,ex1p3,col=3,lwd=2)
>      rm(f2,fx2,x,sigma,y,ex1p0,ex1p1,ex1p2,ex1p3)
> #######################################################
> #                                                     
> #                  bivariate examples                 
> #                                                     
> ####################################################### 
> #  
> #  Local constant image from Polzehl and Spokoiny (2000)
> #
>      xy <- rbind(rep(0:255,256),rep(0:255,rep(256,256)))
>      indw <- c(1:12,29:48,73:100,133:168,209:256)
>      w0 <- matrix(rep(FALSE,256*256),ncol=256)
>      w0[indw,] <- TRUE
>      w0[,indw] <- !w0[,indw]
>      w0 <- w0-.5
>      
>      w0[((xy[1,]-129)^2+(xy[2,]-129)^2)<=10000&((xy[1,]-129)^2+(xy[2,]-129)^2)>=4900] <- 0
>      w0[abs(xy[1,]-xy[2,])<=20&((xy[1,]-129)^2+(xy[2,]-129)^2)<4900] <- 0
>      w0[((xy[1,]-225)^2+2*(xy[2,]-30)^2)-(xy[1,]-225)*(xy[2,]-30)<=625] <- 0
>      
>      w0[((xy[1,]-225)^2+2*(xy[2,]-30)^2)-(xy[1,]-225)*(xy[2,]-30)<=625&xy[2,]>27&xy[2,]<31] <- -.5
>      
>      w0[((xy[1,]-225)^2+2*(xy[2,]-30)^2)-(xy[1,]-225)*(xy[2,]-30)<=625&xy[1,]>223&xy[1,]<227] <- .5
>      w0[((xy[2,]-225)^2+2*(xy[1,]-30)^2)+(xy[2,]-225)*(xy[1,]-30)<=625] <- 0
>      
>      w0[((xy[2,]-225)^2+2*(xy[1,]-30)^2)+(xy[2,]-225)*(xy[1,]-30)<=625&xy[1,]>27&xy[1,]<31] <- -.5
>      
>      w0[((xy[2,]-225)^2+2*(xy[1,]-30)^2)+(xy[2,]-225)*(xy[1,]-30)<=625&xy[2,]>223&xy[2,]<227] <- .5
>      w0[((xy[2,]-225)^2+(xy[1,]-225)^2)+1*(xy[2,]-225)*(xy[1,]-225)<=400] <- 0
>      w0[((xy[2,]-30)^2+(xy[1,]-30)^2)<=256] <- 0
>      rm(xy,indw)
>      set.seed(1)
>      sigma <- .5
>      y <- w0+rnorm(w0,0,sigma)
> #    run some steps interactively  (increase hmax)
>      tmp <- aws(y,graph=TRUE,hmax=2,demo=TRUE)
Press return#    now without graphics and larger hmax
Press return#   increase hmax for better results
Press return     tmp <- aws(y,hmax=5)
Press return     par(mfrow=c(1,3))
Press return     image(y,col=gray((0:255)/255))
Press return     image(tmp$theta,zlim=range(y),col=gray((0:255)/255))
Press return     image(w0,zlim=range(y),col=gray((0:255)/255))
>      rm(y,w0,tmp,sigma)
> #
> #   piecewise smooth image from Polzehl and Spokoiny (2003)
> #
>      x <- (0:99)/99
>      fbi <- function(x,y) (x^2+y^3)*sign(x^2-1*x*y-.75*y^3)
>      z0 <- outer(2*(x-.5),2*(x-.25),FUN=fbi)
>      z <- z0+rnorm(z0,0,.5)
>      par(mfrow=c(1,3),mar=c(3,3,3,.5),mgp=c(2,1,0))
>      set.seed(1)
>      persp(x,x,z0,phi=15,theta=150,d=5,r=2,expand=1.5,col="white")
>      title("True function")
>      persp(x,x,z,phi=15,theta=150,d=5,r=2,expand=1.5,col="white")
>      title("Observed values")
>      image(x,x,z,col=gray((0:255)/255))
>      title("Observed values")
> #
> #   local constant
> #
>      ex1p0 <- aws(z,hmax=3,symmetric=FALSE)$theta # use hmax=5 or larger
> #
> #   local linear
> #
>      ex1p1 <- aws(z,p=1,hmax=3)$theta[1,,]  # use hmax=12 or larger
> #
> #   local quadratic
> #
>      ex1p2 <- aws(z,p=2,hmax=3)$theta[1,,]  # use hmax=20 or larger
> #
> #   display results
> #
>      par(mfrow=c(2,2),mar=c(0.25,.5,.5,.5),mgp=c(.5,.5,0))
>      persp(x,x,z,phi=15,theta=150,d=5,r=2,expand=1,xlab="",ylab="",zlab="",box=FALSE)
>      title("original data",line=-3)
>      persp(x,x,ex1p0,phi=15,theta=150,d=5,r=2,expand=1,xlab="",ylab="",zlab="",box=FALSE)
>      title("local constant AWS",line=-3)
>      persp(x,x,ex1p1,phi=15,theta=150,d=5,r=2,expand=1,xlab="",ylab="",zlab="",box=FALSE)
>      title("local linear AWS (small hmax)",line=-3)
>      persp(x,x,ex1p2,phi=15,theta=150,d=5,r=2,expand=1,xlab="",ylab="",zlab="",box=FALSE)
>      title("local quadratic AWS (small hmax)",line=-3)
>      rm(x,fbi,z0,ex1p0,ex1p1,ex1p2)
> 
> #######################################################
> #
> #                  three-variate examples
> #
> #######################################################
> #
> #  Local constant image from Polzehl and Spokoiny (2000)
> #
>      xy <- rbind(rep(0:30,31),rep(0:30,rep(31,31)))
>      w3 <- array(0,c(31,31,31))
>      w3[4:28,4:28,4:28] <- 1
>      dim(w3) <- c(961,31)
>      w3[((xy[1,]-15)^2+(xy[2,]-15)^2)<=144,16] <- 0
>      for(i in 1:12) {
+         r2 <- 144-i*i
+         w3[((xy[1,]-15)^2+(xy[2,]-15)^2)<=r2,16+c(-i,i)] <- 0
+      }
>      dim(w3) <- c(31,31,31)
>      w3[10:22,10:22,10:22] <- 1
>      dim(w3) <- c(961,31)
>      w3[((xy[1,]-15)^2+(xy[2,]-15)^2)<=36,16] <- 0
>      for(i in 1:6) {
+         r2 <- 36-i*i
+         w3[((xy[1,]-15)^2+(xy[2,]-15)^2)<=r2,16+c(-i,i)] <- 0
+      }
>      dim(w3) <- c(31,31,31)
>      rm(xy)
>      sigma <- .5
>      set.seed(1)
>      y <- w3+rnorm(w3,0,sigma)
> #   increase hmax for reasonable results (hmax >=5)
>      tmp <- aws(y,hmax=2)
>      par(mfrow=c(1,3))
>      for(i in 14:18){
+      image(y[,,i],col=gray((0:255)/255))
+      image(tmp$theta[,,i],zlim=range(y),col=gray((0:255)/255))
+      image(w3[,,i],zlim=range(y),col=gray((0:255)/255))
+ #     readline()
+      }
>      rm(w3,y,tmp,sigma)
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> cleanEx(); ..nameEx <- "awsbi"
> 
> ### * awsbi
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: awsbi
> ### Title: Two-dimensional Adaptive Weights Smoothing
> ### Aliases: awsbi
> ### Keywords: regression nonparametric smooth
> 
> ### ** Examples
> 
> xy<-rbind(rep(0:255,256),rep(0:255,rep(256,256)))
> indw<-c(1:12,29:48,73:100,133:168,209:256)
> w0<-matrix(rep(0,256*256),ncol=256)
> w0[indw,]<-1
> w0[,indw]<-!w0[,indw]
> w0<-w0-.5
> w0[((xy[1,]-129)^2+(xy[2,]-129)^2)<=10000&((xy[1,]-129)^2+(xy[2,]-129)^2)>=4900]<- 0
> w0[abs(xy[1,]-xy[2,])<=20&((xy[1,]-129)^2+(xy[2,]-129)^2)<4900]<- 0
> w0[((xy[1,]-225)^2+2*(xy[2,]-30)^2)-(xy[1,]-225)*(xy[2,]-30)<=625]<- 0
> w0[((xy[1,]-225)^2+2*(xy[2,]-30)^2)-(xy[1,]-225)*(xy[2,]-30)<=625&xy[2,]>27&xy[2,]<31]<- -.5
> w0[((xy[1,]-225)^2+2*(xy[2,]-30)^2)-(xy[1,]-225)*(xy[2,]-30)<=625&xy[1,]>223&xy[1,]<227]<- .5
> w0[((xy[2,]-225)^2+2*(xy[1,]-30)^2)+(xy[2,]-225)*(xy[1,]-30)<=625]<- 0
> w0[((xy[2,]-225)^2+2*(xy[1,]-30)^2)+(xy[2,]-225)*(xy[1,]-30)<=625&xy[1,]>27&xy[1,]<31]<- -.5
> w0[((xy[2,]-225)^2+2*(xy[1,]-30)^2)+(xy[2,]-225)*(xy[1,]-30)<=625&xy[2,]>223&xy[2,]<227]<- .5
> w0[((xy[2,]-225)^2+(xy[1,]-225)^2)+1*(xy[2,]-225)*(xy[1,]-225)<=400]<- 0
> w0[((xy[2,]-30)^2+(xy[1,]-30)^2)<=256]<-0
> sigma<-.25
> y<-w0+rnorm(w0,0,sigma)
> #  increase rmax for better results
> yhat<-awsbi(y,rmax=3)
Control sceme:  1 1 0 1 0 0 
Estimated variance: 0.06823233 
> par(mfrow=c(1,3))
> image(y,col=gray((0:255)/255))
> title("Noisy image")
> image(yhat$yhat,zlim=range(y),col=gray((0:255)/255))
> title("AWS reconstruction")
> image(w0,zlim=range(y),col=gray((0:255)/255))
> title("Original image")
> rm(y,w0,yhat,xy)
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> cleanEx(); ..nameEx <- "awsdens"
> 
> ### * awsdens
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: awsdens
> ### Title: local constant density estimation
> ### Aliases: awsdens
> ### Keywords: smooth nonparametric models
> 
> ### ** Examples
> 
> ###
> ###   univariate density estimation
> ###
> f0 <- function(x) 1.5*(x>=0)-(x>=.25)+(x>=.75)-1.5*(x>1)
> set.seed(1)
> m <- 1000
> x1 <- runif(m)
> ind <- (runif(m)<f0(x1)/1.5)
> n <- 200
> y <- x1[ind][1:n]
> tmp0 <- awsdens(y,440,20,hmax=2000)
> plot(tmp0$xgrid[[1]],tmp0$dens,type="l",lwd=2,
+      ylim =range(0,1.5,tmp0$dens),xlab="x",ylab="Density")
> lines(tmp0$xgrid[[1]],f0(tmp0$xgrid[[1]]),lty=3,col=2,lwd=2)
> title("Estimated and true density (n=200)")
> ###
> ###   bivariate density estimation
> ###
> f1 <- function(x,y) 7.5*pmax(x*(1-x^2-y^2),0)*(x>0)*(y>0)
> set.seed(1)
> m <- 10000
> x1 <- runif(m)
> x2 <- runif(m)
> fx12 <- f1(x1,x2)
> ind <- (runif(m)<f1(x1,x2)/.385/7.5)
> n <- 2500
> y <- rbind(x1[ind],x2[ind])[,1:n]
> tmp <- awsdens(y,120,10,hmax=5)
> image(tmp$xgrid[[1]],tmp$xgrid[[2]],tmp$dens,xlab="x1",ylab="x2",col=gray(.5+(255:0)/511),lwd=2)
> #   thats the estimated density on the grid
> lines(seq(0,1,.01),sqrt(1-seq(0,1,.01)^2),col=2,lty=2,lwd=2)
> lines(c(0,1),c(0,0),col=2,lty=2,lwd=2)
> lines(c(0,0),c(0,1),col=2,lty=2,lwd=2)
> #    thats the boundary of the support
> title("Estimated density (n=2500)")
> #    now add contour lines of the estimated density
> contour(tmp$xgrid[[1]],tmp$xgrid[[2]],tmp$dens,nlevels=20,add=TRUE,col=1,lty=1,labcex=.1)
> rm(f0,m,x1,ind,n,y,tmp0,f1,x2,fx12,tmp)
> 
> 
> 
> cleanEx(); ..nameEx <- "awsh"
> 
> ### * awsh
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: awsh
> ### Title: Univariate local polynomial Adaptive Weights Smoothing for
> ###   regression with heteroscedastic additive errors
> ### Aliases: awsh
> ### Keywords: smooth nonparametric regression models
> 
> ### ** Examples
> 
> ##---- Should be DIRECTLY executable !! ----
> ##-- ==>  Define data, use random,
> ##--    or do  help(data=index)  for the standard data sets.
> 
> 
> 
> cleanEx(); ..nameEx <- "awstindex"
> 
> ### * awstindex
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: awstindex
> ### Title: Tail index estimation
> ### Aliases: awstindex
> ### Keywords: smooth nonparametric models
> 
> ### ** Examples
> 
> ###
> ###   Estimate the tail-index of a cauchy distribution
> ###   absolute values can be used because of the symmetry of centered cauchy
> ###
> set.seed(1)
> n <- 500
> x <- rcauchy(n)
> tmp <- awstindex(abs(x),hmax=n)
> tmp$tindex
[1] 1.243385
> ###
> ###   now show the segmentation generated by AWS 
> ###
> plot(tmp$intensity[1:250],type="l")
> 
> 
> 
> cleanEx(); ..nameEx <- "awstri"
> 
> ### * awstri
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: awstri
> ### Title: Three-dimensional Adaptive Weights Smoothing
> ### Aliases: awstri
> ### Keywords: regression nonparametric smooth
> 
> ### ** Examples
> 
> xy <- rbind(rep(0:30,31),rep(0:30,rep(31,31)))
> w3 <- array(0,c(31,31,31))
> w3[4:28,4:28,4:28] <- 1
> dim(w3) <- c(961,31)
> w3[((xy[1,]-15)^2+(xy[2,]-15)^2)<=144,16] <- 0
> for(i in 1:12) {
+    r2 <- 144-i*i
+    w3[((xy[1,]-15)^2+(xy[2,]-15)^2)<=r2,16+c(-i,i)] <- 0
+ }
> dim(w3) <- c(31,31,31)
> w3[10:22,10:22,10:22] <- 1
> dim(w3) <- c(961,31)
> w3[((xy[1,]-15)^2+(xy[2,]-15)^2)<=36,16] <- 0
> for(i in 1:6) {
+    r2 <- 36-i*i
+    w3[((xy[1,]-15)^2+(xy[2,]-15)^2)<=r2,16+c(-i,i)] <- 0
+ }
> dim(w3) <- c(31,31,31)
> sigma <- .4
> y <- w3+rnorm(w3,0,sigma)
> #  increase rmax for better results
> yhat <- awstri(y,rmax=2)
Control sceme:  1 1 0 1 
> rm(y,yhat,w3,xy)
> 
> 
> 
> cleanEx(); ..nameEx <- "awsuni"
> 
> ### * awsuni
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: awsuni
> ### Title: One-dimensional Adaptive Weights Smoothing
> ### Aliases: awsuni
> ### Keywords: regression nonparametric smooth
> 
> ### ** Examples
> 
> #  Blocks data (from Donoho, Johnstone, Kerkyacharian and Picard (1995))
> mofx6 <- function(x){
+ xj <- c(10,13,15,23,25,40,44,65,76,78,81)/100
+ hj <- c(40,-50,30,-40,50,-42,21,43,-31,21,-42)*.37
+ Kern <- function(x) (1-sign(x))/2
+ apply(Kern(outer(xj,x,"-"))*hj,2,sum)
+ }
> x <- seq(0,1,1/2047)
> fx6 <- mofx6(x)
> #    sigma==3
> y <- rnorm(fx6,fx6,3)
> tmp <- awsuni(y)
Control sceme:  1 1 0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
> par(mfrow=c(1,1))
> plot(x,y)
> lines(x,tmp$yhat,col=2)
> lines(x,fx6,col=3)
> title(expression(paste("AWS Reconstruction of blocks data  ",sigma==3)))
> rm(x,y,fx6,mofx6,tmp)
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> cleanEx(); ..nameEx <- "getnubi"
> 
> ### * getnubi
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: getnubi
> ### Title: Calculate number of grid points in a bivariate elliptical
> ###   neighborhood
> ### Aliases: getnubi
> ### Keywords: misc
> 
> ### ** Examples
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "getnutri"
> 
> ### * getnutri
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: getnutri
> ### Title: Calculate number of grid points in a 3D ellipsoidal neighborhood
> ### Aliases: getnutri
> ### Keywords: misc
> 
> ### ** Examples
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "laws"
> 
> ### * laws
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: laws
> ### Title: Likelihood based Adaptive Weights Smoothing
> ### Aliases: laws
> ### Keywords: smooth nonparametric regression models
> 
> ### ** Examples
> 
> ###
> ###    Artificial PET data
> ###
> x <- 1:128
> f12 <- function(x,y){
+ 2*((1.5*(x-64)^2+(y-64)^2<3025)) +
+ ((x-64)^2+(y-104)^2<16)-((x-92)^2+(y-84)^2<25)+
+ ((x-78)^2+(y-84)^2<25)-((x-64)^2+(y-84)^2<36)+
+ ((x-50)^2+(y-84)^2<36)-((x-36)^2+(y-84)^2<25)+
+ ((x-92)^2+(y-64)^2<25)-((x-78)^2+(y-64)^2<16)+
+ ((x-64)^2+(y-64)^2<16)-((x-50)^2+(y-64)^2<25)+
+ ((x-36)^2+(y-64)^2<25)-((x-92)^2+(y-44)^2<36)+
+ ((x-78)^2+(y-44)^2<36)-((x-64)^2+(y-44)^2<25)+
+ ((x-50)^2+(y-44)^2<25)-((x-36)^2+(y-44)^2<16)+
+ ((x-64)^2+(y-24)^2<16)
+ }
> u0 <- 2*outer(x,x,"f12")
> set.seed(1)
> y <- matrix(rpois(u0,u0),128,128)
> # use larger hmax for good results
> yhat <- laws(y,model="Poisson",hmax=4)$theta
> par(mfrow=c(1,3),mar=c(3,3,3,.5),mgp=c(2,1,0))
> image(y,col=gray((0:255)/255))
> title("Observed image")
> image(yhat,col=gray((0:255)/255))
> title("AWS-Reconstruction")
> image(u0,col=gray((0:255)/255))
> title("True image")
> rm(u0,y,yhat,x)
> ###
> ###   VOLATITILTY ESTIMATION
> ###
> ###   artificial example
> ###
> sigma <- c(rep(1,125),rep(2,125),rep(.5,125),rep(1,125))
> set.seed(1)
> x <- rnorm(sigma,0,sigma)
> tmpa <- laws(x,model="Volatility",u=sigma,graph=TRUE,hmax=250)
bandwidth:  1    MSE:  0.2706934    MAE:  0.3883171 
bandwidth:  1.25    MSE:  0.2462917    MAE:  0.3493625 
bandwidth:  1.56    MSE:  0.2458382    MAE:  0.3464986 
bandwidth:  1.95    MSE:  0.2508973    MAE:  0.3594376 
bandwidth:  2.44    MSE:  0.2069912    MAE:  0.3254337 
bandwidth:  3.05    MSE:  0.1748197    MAE:  0.3005164 
bandwidth:  3.81    MSE:  0.1391698    MAE:  0.269995 
bandwidth:  4.77    MSE:  0.1100888    MAE:  0.2393491 
bandwidth:  5.96    MSE:  0.08903676    MAE:  0.2132892 
bandwidth:  7.45    MSE:  0.0727766    MAE:  0.1919125 
bandwidth:  9.31    MSE:  0.06079471    MAE:  0.1758755 
bandwidth:  11.6    MSE:  0.05184914    MAE:  0.1615062 
bandwidth:  14.6    MSE:  0.04492663    MAE:  0.1471510 
bandwidth:  18.2    MSE:  0.03926773    MAE:  0.1356497 
bandwidth:  22.7    MSE:  0.03520471    MAE:  0.1254111 
bandwidth:  28.4    MSE:  0.03348095    MAE:  0.1183937 
bandwidth:  35.5    MSE:  0.03362099    MAE:  0.1160992 
bandwidth:  44.4    MSE:  0.03477255    MAE:  0.1140405 
bandwidth:  55.5    MSE:  0.03699834    MAE:  0.1109075 
bandwidth:  69.4    MSE:  0.03959886    MAE:  0.1076274 
bandwidth:  86.7    MSE:  0.04195846    MAE:  0.1053840 
bandwidth:  108    MSE:  0.04432575    MAE:  0.1037195 
bandwidth:  136    MSE:  0.04587928    MAE:  0.1019359 
bandwidth:  169    MSE:  0.04682174    MAE:  0.1005104 
bandwidth:  212    MSE:  0.04771838    MAE:  0.09935786 
> tmps <- laws(x,model="Volatility",u=sigma,hmax=250,symmetric=TRUE)
bandwidth:  1    MSE:  0.2706934    MAE:  0.3883171 
> par(mfrow=c(1,1),mar=c(3,3,3,1))
> plot(abs(x),col=3,xlab="time t",ylab=expression(abs( R[t] )))
> lines(tmpa$theta,col=1,lwd=2)
> lines(tmps$theta,col=1,lwd=2,lty=2)
> lines(sigma,col=2,lwd=2,lty=3)
> legend(350,5.5,c("asymmetric AWS","symmetric AWS","true sigma"),
+                  lwd=c(2,2,2),lty=c(1,2,3),col=c(1,1,2))
> title(expression(paste("Estimates of  ",sigma(t))))
> rm(tmpa,tmps,sigma,x)
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  33.69 3.01 37.51 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')