sfMM-rhoLvfree: Estimation of the local regularization parameters
In bozenne/MRIaggr: Management, Display, and Processing of Medical Imaging Data

Description Usage Arguments Details Value References See Also Examples

Estimation of the local regularization parameters on external data using an efficient MCMC procedure.

rhoLvfree(Y, W_SR, rho_max = 10, rho_init = "PML", sd_rho_init = "PML", site_order = NULL,
     dprior = function(x){stats::dunif(x,0,rho_max)}, epsilon = 0.001, update_epsilon = 1,
	 n.sample = 2500, iter_max = 3, burn_in = round(0.25 * n.sample), thin = 1,
	 n.chain = 1, verbose = 3, cpus = 1, export.coda = TRUE)

`Y`	a matrix containing the observations (by rows) for the various groups (by columns). REQUIRED.
`W_SR`	the local neighbourhood matrix. dgCMatrix. Should be normalized by row (i.e. `rowSums(W_SR)=1`). REQUIRED.
`rho_max`	maximum possible rho value, minimum is 0. double.
`rho_init`	the initial rho value. double. If `NULL`, it is sample in a uniform law between 0 and `rho_max`
`sd_rho_init`	the standard deviation of the rho value for the proposal distribution. double. It is updated at the end of the burn in.
`site_order`	a specific order to go all over the sites. Can be `NULL` or an integer vector.
`dprior`	the prior distribution. function. Must be a density function between 0 and `rho_max`.
`epsilon`	the tolerance parameter regulating the acceptance of the parameter. double.
`update_epsilon`	if over 1 the value by which the epsilon parameter is devided every 5% of the total iterations, if under one the acceptance rate to reach. double.
`n.sample`	the number of iterations of the Gibbs sampler. integer.
`iter_max`	the number of Gibbs move for each simulation. integer.
`burn_in`	the number of iteration of the burn in phase. integer.
`thin`	the thinning interval between consecutive observations. integer.
`n.chain`	the number of chain to use. integer.
`verbose`	how should the execution of the function should be traced ? 1 traces display the initialization and the final result, `2` traces every 5% of the total iterations and `3` also displays the update performed every 5%.
`cpus`	the number of CPU to use. integer.
`export.coda`	should the results be convert to the mcmc format of the coda package ? logical.

FUNCTION:
This function uses a likelihood-free Metropolis-Hastings method (proposed by Pereyra et al., 2013) to estimate the local regularization parameter. It should give a less biased estimation of the parameter but require a higher computational cost.

A good initialisation (rho_init) may be obtained with the rhoMF function. It is automatically done if rho_init is set to "PML".

Setting sd_rho_init to "PML" leads to adjustement of the initial variance according to the rho value and the sample size.

M. Pereyra, N. Dobigeon, H. Batatia, and J.Y. Tourneret. Estimation the granularity coefficient of a Potts-Markov random field within an MCMC algorithm. IEEE Trans. Image Porcessing, 22(6):2385-2397, 2013.

calcW to compute the neighbourhood matrix,
simulPotts to simulate from a Potts model. rhoMF to estimate the regularization parameters using mean field approximation. calcPottsParameter general interface for estimating the regularization parameters.

# spatial field
## Not run: 
n <- 50

## End(Not run)

G <- 3
coords <- which(matrix(0, nrow = n * G, ncol = n * G) == 0, arr.ind = TRUE)

# neighbourhood matrix
resW <- calcW(as.data.frame(coords), range = sqrt(2), row.norm = TRUE, calcBlockW = TRUE)
W_SR <- resW$W
site_order <- unlist(resW$blocks$ls_groups)-1

# initialisation
set.seed(10)
sample <- simulPotts(W_SR, G = 3, rho = 3.5, iter_max = 500, 
                     site_order = site_order)$simulation

multiplot(as.data.frame(coords), sample,palette = "rgb")

# estimation

## Not run: 
rho <- rhoLvfree(Y = sample, W_SR = W_SR, site_order = site_order)

## End(Not run)