laws {aws} R Documentation

## 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` `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 `y` determines if one, two or three-dimensional AWS is used. `x` `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 a Chi-square-distribution with number of parameters in the polynomial model as degrees of freedom. If `qlambda==NULL` a standard value depending on `model` and `symmetric` is choosen. `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"`. The Kernel operates on the squared distance, so `"Triangle"` corresponds to the use of an Epanechnikov kernel in kernel smoothing. `"Exponential"` requires larger values of `hmax` and therefore more iterations to reach comparable results. `model` `model` determines the distribution type of `y`. Currently implemented models are `"Poisson"`, `"Bernoulli"`, `"Gaussian"`, `"Exponential"`, `"Weibull"`, `"Volatility"` (Estimation of the scale parameter of a Gaussian distribution). Defaults to `"Poisson"`. `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.

### 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`):

Binary response
`model="Bernoulli"`
Poisson regression
`model="Poisson"` This model allows e.g. for density estimation or for the analysis of poisson count data on a grid (e.g. Positron emission tomography (PET)).
Exponential regression
`model="Exponential"` Applications of this model include e.g. test for constant (over time) failure rates and estimation of tail indicies.
Gaussian regression
`model="Gaussian"` This essentially coincides with the local constant regression model with additive sub-gaussian errors as described in Polzehl and Spokoiny (2000, 2003).
Weibull regression
`model="Weibull"` Applications in reliability analysis.
Volatility model
`model="Volatility"` .

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.

### Value

 `theta ` Parameter estimates, first dimension corresponds to parameter components `y` values provided in `y` `x` values provided in `x` `call` actual function call

### 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 `aws`, `awsdens`, `awstindex`

### 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)
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)
```

[Package aws version 1.3-0 Index]