LaplaceDeconv (version 1.0.4)

LagLaplDeconv: function LagLaplDeconv

Description

Main function of this package : computes the Laplace deconvolution with noisy discrete non-equally spaced observations on a finite time interval (see reference).

Usage

LagLaplDeconv(Y, g, times = 1:length(Y), sigma, cpen = 2, atab = -1, Mmax = 25, ncores.max = Inf, verbose = FALSE, withplot = FALSE)

Arguments

Y

numeric vector, the observed noisy observations of the Laplace convolution

g

numeric vector, the known kernel of the Laplace convolution

times

numeric vector, the observation times (default 1:length(Y))

sigma

numeric, the noise level

cpen

numeric, the penalization constant (default 2)

atab

numeric vector, an array of value for a (default -1, see Details)

Mmax

integer, the maximum degree (default 25)

ncores.max

numeric, number of used cores (default Inf)

verbose

boolean to control information output (default FALSE)

withplot

boolean to control plot output (default FALSE)

Value

a list containing:

• f.hat, numeric vector, the estimate at observation times
• q.hat, numeric vector, the reconstructed convolution of the kernel g and the estimate f.hat, computed at observation times
• f.coef, numeric vector, the coefficients of f.hat in the selected Laguerre function basis
• info, only for internal use --- not documented
• a.hat, numeric, the selected value of the parameter a
• M0, numeric, the dimension of the selected model
• sigma, numeric, the noise level
• cpen, numeric, the penalization constant used in the penalty

Details

atab defines the values of the scale parameter used to build the Laguerre functions basis. If atab is of length 1 and negative then atab<-seq(0.4,3,by=0.05/abs(atab))/sqrt(Tmax/10) where Tmax=max(times).

References

Laplace deconvolution on the basis of time domain data and its application to Dynamic Contrast Enhanced imaging by F. Comte, C-A. Cuenod, M. Pensky, Y. Rozenholc (ArXiv http://arxiv.org/abs/1405.7107)

Examples

## Not run: 
#  #### AN ARTICIAL EXAMPLE ####
# 
#  library(LaplaceDeconv)
#  par(mfrow=c(1,1))
#  set.seed(29102015)
# 
#  sigma=0.02
#  a = 1
#  t = seq(0,5,l=100)
#  g = 20*t^2*exp(-5*t)
#  f.coef = c(0.4,0.02,0.01)
# 
#  # compute the Laplace convolution from g, kernel computed at times t, and the function
#  # described by its decomposition in Laguerre function basis with scale a :
#  fg = LaguerreLaplaceConvolution(t,g,f.coef,a)
# 
#  # the noisy observations :
#  Y = fg+sigma*rnorm(length(fg))
# 
#  # estimation of f from the observation and the kernel :
#  L = LagLaplDeconv(Y,g,t,sigma)
#  matplot(t,cbind(g,MakeLaguerreMatrix(a,3)(t)%*%f.coef,fg,L$q.hat,L$f.hat,Y),lty=1,
#    type=c('b',rep('l',4),'p'),ylab='',pch='x')
# 
#  # display results of estimation
#  legend('topright',lty=c(rep(1,5),0),pch=c('x',rep('',4),'x'),
#    legend=c(
#      'g: partially observed kernel',
#      'f: unknown',
#      'q=fxg: unknown convolution',
#      expression(hat(q)*': plug-in convolution'),
#      expression(hat(f)*': estimation of f'),
#      'Y: observations'),
#    col=1:6)
#  ## End(Not run)

 ## Not run: 
#  #### A REAL EXAMPLE USING DCE-MRI DATA FROM A TUMOR ####
# 
#  library(LaplaceDeconv)
#  par(mfrow=c(1,2))
# 
#  # load data from patient before the treatment
#  data(EX_DCEMRI_t0)
# 
#  # display AIF and tumoral enhancements
#  matplot(ex_dcemri$times,
#    cbind(ex_dcemri$AIF,ex_dcemri$TUM_1,ex_dcemri$TUM_2,ex_dcemri$TUM_3),
#    ylab='',lty=1,type=c('b',rep('p',3)),pch='+',main='Observations')
#  legend('topright',pch='+',legend=c('AIF','TUM_1','TUM_2','TUM_3'),col=1:4)
# 
#  # estimation of the contrast agent survival functions
#  L1 = LagLaplDeconv(ex_dcemri$TUM_1,ex_dcemri$AIF,ex_dcemri$times,ex_dcemri$sigma)
#  L2 = LagLaplDeconv(ex_dcemri$TUM_2,ex_dcemri$AIF,ex_dcemri$times,ex_dcemri$sigma)
#  L3 = LagLaplDeconv(ex_dcemri$TUM_3,ex_dcemri$AIF,ex_dcemri$times,ex_dcemri$sigma)
# 
#  matlines(ex_dcemri$times,cbind(L1$q.hat,L2$q.hat,L3$q.hat),type='l',lty=1,col=2:4)
# 
#  # display results of estimation
#  matplot(ex_dcemri$times,cbind(L1$f.hat,L2$f.hat,L3$f.hat),type='l',lty=1,col=2:4,
#    ylab='survival',main='Contrast agent survival fcts')
#  legend('topright',lty=1,col=2:4,
#    legend=c(
#      paste0('TUM_1 - a.hat=',round(L1$a.hat,digits=2)),
#      paste0('TUM_2 - a.hat=',round(L2$a.hat,digits=2)),
#      paste0('TUM_3 - a.hat=',round(L3$a.hat,digits=2))
#      )
#    )
#  ## End(Not run)