Adaptive weigths smoothing using patches

The function implements a version the propagation separation approach that uses patches instead of individuel voxels for comparisons in parameter space. Functionality is analog to function aws. Using patches allows for an improved handling of locally smooth functions and in 2D and 3D for improved smoothness of discontinuities at the expense of increased computing time.

paws(y, hmax = NULL, mask=NULL, onestep = FALSE, aws = TRUE, family = "Gaussian",
     lkern = "Triangle", aggkern = "Uniform", sigma2 = NULL, shape = NULL,
     scorr = 0, spmin = 0.25, ladjust = 1, wghts = NULL, u = NULL,
     graph = FALSE, demo = FALSE, patchsize = 1)

Arguments

y	array y containing the observe response (image intensity) data. dim(y) determines the dimensionality and extend of the grid design.
mask	logical array defining a mask. All computations are restricted to the mask.
hmax	hmax specifies the maximal bandwidth. Defaults to hmax=250, 12, 5 for 1D, 2D, 3D images, respectively. In case of lkern="Gaussian" the bandwidth is assumed to be given in full width half maximum (FWHM) units, i.e., 0.42466 times gridsize.
onestep	apply the last step only (use for test purposes only)
aws	logical: if TRUE structural adaptation (AWS) is used.
family	family specifies the probability distribution. Default is family="Gaussian", also implemented are "Bernoulli", "Poisson", "Exponential", "Volatility", "Variance" and "NCchi". family="Volatility" specifies a Gaussian distribution with expectation 0 and unknown variance. family="Volatility" specifies that p*y/theta is distributed as \(\chi^2\) with p=shape degrees of freedom. family="NCchi" uses a noncentral Chi distribution with p=shape degrees of freedom and noncentrality parameter theta
lkern	character: location kernel, either "Triangle", "Plateau", "Quadratic", "Cubic" or "Gaussian". The default "Triangle" is equivalent to using an Epanechnikov kernel, "Quadratic" and "Cubic" refer to a Bi-weight and Tri-weight kernel, see Fan and Gijbels (1996). "Gaussian" is a truncated (compact support) Gaussian kernel. This is included for comparisons only and should be avoided due to its large computational costs.
aggkern	character: kernel used in stagewise aggregation, either "Triangle" or "Uniform"
sigma2	sigma2 allows to specify the variance in case of family="Gaussian". Not used if family!="Gaussian". Defaults to NULL. In this case a homoskedastic variance estimate is generated. If length(sigma2)==length(y) then sigma2 is assumed to contain the pointwise variance of y and a heteroscedastic variance model is used.
shape	Allows to specify an additional shape parameter for certain family models. Currently only used for family="Variance", that is \(\chi\)-Square distributed observations with shape degrees of freedom.
scorr	The vector scorr allows to specify a first order correlations of the noise for each coordinate direction, defaults to 0 (no correlation).
spmin	Determines the form (size of the plateau) in the adaptation kernel. Not to be changed by the user.
ladjust	factor to increase the default value of lambda
wghts	wghts specifies the diagonal elements of a weight matrix to adjust for different distances between grid-points in different coordinate directions, i.e. allows to define a more appropriate metric in the design space.
u	a "true" value of the regression function, may be provided to report risks at each iteration. This can be used to test the propagation condition with u=0
graph	If graph=TRUE intermediate results are illustrated after each iteration step. Defaults to graph=FALSE.
demo	If demo=TRUE the function pauses after each iteration. Defaults to demo=FALSE.
patchsize	positive integer defining the size of patches. Number of grid points within the patch is (2*patchsize+1)^d with d denoting the dimensionality of the design.

Details

see aws. The procedure is supposed to produce superior results if the assumption of a local constant image is violated or if smooothness of discontinuities is desired.

Value

returns an object of class aws with slots

y = "numeric"

y

dy = "numeric"

dim(y)

x = "numeric"

numeric(0)

ni = "integer"

integer(0)

mask = "logical"

logical(0)

theta = "numeric"

Estimates of regression function, length: length(y)

hseq = "numeric"

sequence of bandwidths employed

mae = "numeric"

Mean absolute error for each iteration step if u was specified, numeric(0) else

psnr = "numeric"

Peak signal-to-noise ratio for each iteration step if u was specified, numeric(0) else

var = "numeric"

approx. variance of the estimates of the regression function. Please note that this does not reflect variability due to randomness of weights.

xmin = "numeric"

numeric(0)

xmax = "numeric"

numeric(0)

wghts = "numeric"

numeric(0), ratio of distances wghts[-1]/wghts[1]

degree = "integer"

0

hmax = "numeric"

effective hmax

sigma2 = "numeric"

provided or estimated error variance

scorr = "numeric"

scorr

family = "character"

family

shape = "numeric"

shape

lkern = "integer"

integer code for lkern, 1="Plateau", 2="Triangle", 3="Quadratic", 4="Cubic", 5="Gaussian"

lambda = "numeric"

effective value of lambda

ladjust = "numeric"

effective value of ladjust

aws = "logical"

aws

memory = "logical"

memory

homogen = "logical"

homogen

earlystop = "logical"

FALSE

varmodel = "character"

"Constant"

vcoef = "numeric"

numeric(0)

call = "function"

the arguments of the call to aws

References

J. Polzehl, K. Tabelow (2019). Magnetic Resonance Brain Imaging: Modeling and Data Analysis Using R. Springer, Use R! series. Appendix A. Doi:10.1007/978-3-030-29184-6.

J. Polzehl, K. Papafitsoros, K. Tabelow. Patch-wise adaptive weights smoothing, Preprint no. 2520, WIAS, Berlin, 2018, DOI 10.20347/WIAS.PREPRINT.2520. (to appear in Journal of Statistical Software).

Author

Joerg Polzehl, polzehl@wias-berlin.de, http://www.wias-berlin.de/people/polzehl/

Note

use setCores='number of threads' to enable parallel execution.

See also

See also aws, lpaws, vpaws,link{awsdata}

Examples

if (FALSE) {
setCores(2)
y <- array(rnorm(64^3),c(64,64,64))
yhat <- paws(y,hmax=6)
}