laws: Likelihood based Adaptive Weights Smoothing

Description

This function implements a Adaptive Weights Smoothing procedure for local constant Poisson, Bernoulli, Exponential, Weibull, Volatility and Gaussian models as described in Polzehl & Spokoiny (2002).

Usage

laws(y, x = NULL, qlambda = NULL, eta = 0.5, lkern = "Triangle", 
     model = "Poisson", shape = NULL, hinit = NULL, hincr = NULL, 
     hmax = 10, NN = FALSE, u = NULL, graph = FALSE, demo = FALSE, 
     symmetric = FALSE, wghts=NULL)

Arguments

y contains the observed values at location x. In case of x==NULL (second parameter) y is assumed to be observed on a one, two or three-dimensional grid. The dimension of

x is either NULL, in this case y is assumed to be observed on a grid, or is a matrix, with rows corresponding to variables, containing the design points where y is observed.

qlambda

qlambda determines the scale parameter qlambda for the stochastic penalty. The scaling parameter in the stochastic penalty lambda is choosen as the qlambda-quantile of

eta

eta is a memory parameter used to stabilize the procedure. eta has to be between 0 and 1, with eta==0.5 being the default.

lkern

lkern determines the location kernel to be used. Options are "Uniform", "Triangle", "Quadratic", "Cubic" and "Exponential". Default is "Triangle"

model

model determines the distribution type of y. 
           Currently implemented models are "Poisson", "Bernoulli",
           "Gaussian", "Exponential", "Weibull", <

shape

used for additional parameters of the specified distribution if needed,
           i.e. variance if model=="Gaussian"

hinit

hinit Initial bandwidth for the location penalty. 
           Appropriate value is choosen in case of hinit==NULL

hincr

hincr hincr^(1/d), with d the 
           dimensionality of the design, is used as a factor to increase the 
           bandwidth between iterations. Defauts to hincr==1.2

hmax

hmax Maximal bandwidth to be used. Determines the 
           number of iterations and is used as the stopping rule.

NN

If NN==TRUE use nearest neighbor-rules instead of 
           distances in the location term.

u

u used to supply values of the true regression function 
           for test purposes to calculate  Mean Squared Error (MSE) and 
           Mean Absolute Error (MAE)

graph

graph if TRUE results are displayed after each 
           iteration step.

demo

demo if TRUE after each iteration step results 
           are displayed and the process waits for user interaction.

symmetric

If symmetric==TRUE the stochastic penalty is
           symmetrized, i.e. (sij + sji)/lambda is used instead of 
           sij/lambda. See references for details.

wghts

Specifies wghts for distance evaluation on a bi- or trivariate grid.
                Allows for anisotropic local neighborhoods. If wghts=NULL
                isotropic neighborhoods are used.

`Value`

thetaParameter estimates, first dimension corresponds to parameter 
components
yvalues provided in y
xvalues provided in x
callactual function call

`Details`

This function implements an adaptive weights 
   smoothing (AWS) procedure for a several classes of distributions for the dependent
   variable in local constant regression models. The approach generalizes the original 
   AWS procedure from Polzehl and Spokoiny (2000).
   
   Adaptive weights smoothing is an iterative data adaptive smoothing technique that
   is designed  for smoothing in regression problems with discontinuous regression
   function. The basic assumption is that the regression function can be approximated 
   by a simple local,  e.g. local constant or local polynomial, model. 
   The estimate of the regression function, i.e. the conditional expectation of y
   given x is computed as a weighted maximum likelihood estimate, with weights choosen
   in a completely data adaptive way. The procedure is edge preserving. If the assumed local
   model is globally valid, almost all weights used will be 1, i.e. the resulting estimate
   almost is the global estimate.
    
      
   Currently implemented are the following models (specified by parameter model):
   [object Object],[object Object],[object Object],[object Object],[object Object],[object Object]
   
   The essential parameter in the procedure is qlambda. This parameter has an 
   interpretation as a significance level of a test for equivalence of two local
   parameter estimates. Optimal values mainly depend on the choosen model
   and the value of symmetric (determines the use of an asymmetric or a symmetrized
   test). The optimal values only slightly depend on the model parameters, i.e. the
   default parameters should work in most situations. Larger values of qlambda
   may lead to oversmoothing, small values of qlambda lead to a random segmentation
   of homogeneous regions. A good value of qlambda can be obtained by the propagation
   condition, requiring that in case of global validity of the local model the 
   estimate for large hmax should be equal to the global estimate. 
   
   The numerical complexity of the procedure is mainly determined by hmax. The number
   of iterations is d*log(hmax/hinit)/log(hincr) with d being the dimension 
   of y. Comlexity in each iteration step is Const*hakt*n with hakt
   being the actual bandwith in the iteration step and n the number of design points.
   hmax determines the maximal possible variance reduction.

`References`

{Polzehl, J. and Spokoiny, V. (2003). Varying coefficient 
regression modeling 
     by adaptive weights smoothing, WIAS-Preprint 818.}
{Polzehl, J. and Spokoiny, V. (2002). Local likelihood 
modelling by adaptive weights smoothing, WIAS-Preprint 787.}
{ Polzehl, J. and Spokoiny, V. (2000). Adaptive Weights Smoothing
     with applications to image restoration, J.R.Statist.Soc. B, 62,
     Part 2, pp. 335-354}

`See Also`

SEE ALSO aws, awsdens, awstindex

`Examples`

Run this code###
###    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)
tmps <- laws(x,model="Volatility",u=sigma,hmax=250,symmetric=TRUE)
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)
Run the code above in your browser using DataLab