N.I.M.R.O.D.  
Functions/Subroutines

optimization.f90 File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine optim (b, m, ni, v, rl, istop, conv, llikelihood, LCV_aOUT)
 Main optimization (RVS and marquardt)
subroutine LCVa (b, m, LCV_a, llikelihood)
 Quality of adjustement evaluation.
subroutine computeMeasurmentError (b, err)

Function Documentation

subroutine computeMeasurmentError ( double precision,dimension(npm),intent(in)  b,
double precision,dimension(npmcomp),intent(out)  err 
)

Definition at line 820 of file optimization.f90.

References WorkingSharedValues::b1, WorkingSharedValues::bayes_data, WorkingSharedValues::emp_bayes, WorkingSharedValues::numpat1, solution(), WorkingSharedValues::startsauv, and WorkingSharedValues::systeme.

Referenced by optim().

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine LCVa ( double precision,dimension(m),intent(in)  b,
integer,intent(in)  m,
double precision,intent(out)  LCV_a,
double precision,intent(inout)  llikelihood 
)

Quality of adjustement evaluation.

AUTHOR : Melanie Prague Daniel Commenges Julia Drylewicz Jeremy guedj Rodolphe Thiebaut

DESCRIPTION :

Commenges et al. 2006 (A newton- like algorithm for likelihood maximization : The robust-variance scoring algorithm. Arxiv math/0610402) proposed an approximate cross-validation criterion that can be used to choose estimators based on penalized likelihood. We propose to use this criterion for model choice:

\[ \rm{LCV_a} = -n^{-1}[L(\hat \theta)- {\rm Tr} (H^{-1}_{L^P}H_{_L})], \nonumber \]

where $H_{_L}$ and $H_{L^P}$ are the Hessians of minus the log-likelihood and the penalized log-likelihood respectively, taken in the optimized value of parameters.

MODIFICATION :

01/09/2012 - Prague - Refactoring

INFORMATIONS:

Parameters:
[in,out]bInitial starting points then optimized value.
[in]moptimized parameters vector lenth
[out]llikelihoodunpenalized log-likelihood
[out]LCV_aModel quality of fit assessment

Definition at line 757 of file optimization.f90.

References dsinv(), WorkingSharedValues::hessianNonPenalisee, WorkingSharedValues::hessianPenalisee, WorkingSharedValues::LLnonPenalisee, WorkingSharedValues::nbpatOK, and WorkingSharedValues::withexclusion.

Referenced by optim().

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine optim ( double precision,dimension(m),intent(inout)  b,
integer,intent(in)  m,
integer,intent(out)  ni,
double precision,dimension(m*(m+3)/2),intent(out)  v,
double precision,intent(out)  rl,
integer,intent(out)  istop,
double precision,dimension(3),intent(out)  conv,
double precision,intent(out)  llikelihood,
double precision,dimension(2),intent(out)  LCV_aOUT 
)

Main optimization (RVS and marquardt)

AUTHOR : Melanie Prague Daniel Commenges Julia Drylewicz Jeremy guedj Rodolphe Thiebaut

DESCRIPTION :

OPTIMIZATION METHODS : Under the conditions of the Bernstein-Von Mises theorem, the posterior tends to a normal, thus, when $n$ is large enough the penalized log-likelihood is close to a quadratic form in $\theta$. The use of the Newton-Raphson method is justified as this method is based on a quadratic approximation of the surface to maximize. This justify the use of the Robust variance Scoring algorithm (RVS). However with $n$ not very large, the surface may be close to a quadric only in a region near the maximum and may not be convex everywhere. We implemented a switch to a more robust methods : the Levenberg-Marquardt algorithm.

CONVERGENCE CRITERIA : Any iterative algorithm must have a stopping rule; We implemented two commonly used criteria to assess the stability of the algorithm concerning the displacement in the parameter space $||\theta_{k+1}-\theta_k||$ (ca) and the variation of the objective function, here: $L^P(\theta_{k+1})-L^P(\theta_k)$ (cb). If both of these criteria take very low values, this means that the algorithm does not move any more. The thrid and most worth used criterium is the Relative Distance to maximum (rdm), it is characterized by $U^P(\hat{\theta})=0$. It was initially proposed in Commenges et al. 2006 (A newton- like algorithm for likelihood maximization : The robust-variance scoring algorithm. Arxiv math/0610402) :

\[ \rm{RDM}(\theta_k)=\frac{U^P(\theta_k)^TG^{-1}(\theta_k)U^P(\theta_k)}{m}\nonumber \]

In the case where $G$ (the hessian approximation) is not inversible, RDM can not be computed and the program will stop only when the maximum number of iterations is reached; this case is a failure of convergence. Nevertheless an evaluation fo the gradients steepness is provided. Interpretation the optimization results REALLY depends on the exit status istop: istop=1 All convergence criteria are met istop=2 Maximum number of iteration reached istop=3 Failure, the hessian is non inversible. istop=4 "no convergence" program is stuck (restart is propably needed)

MODEL QUALITY ASSESSMENT : Commenges et al. 2006 (A newton- like algorithm for likelihood maximization : The robust-variance scoring algorithm. Arxiv math/0610402) proposed an approximate cross-validation criterion that can be used to choose estimators based on penalized likelihood. We propose to use this criterion for model choice:

\[ \rm{LCV_a} = -n^{-1}[L(\hat \theta)- {\rm Tr} (H^{-1}_{L^P}H_{_L})], \nonumber \]

where $H_{_L}$ and $H_{L^P}$ are the Hessians of minus the log-likelihood and the penalized log-likelihood respectively, taken in the optimized value of parameters. The criterion is implemented in NIMROD, the lower the better.

PARAMETER SIGNIFICANCE : We implemented a Bayesian p-value as an analogue to the frequentist one. For the one-sided test $\beta \leq 0$, it is simply $P(\beta \leq 0|Y)$. For the two-sided Bayesian p-value testing $\beta=0$, called $\alpha_{min}$, it is such that the largest HPD credible set not containing zero has probability $1-\alpha_{min}$. This can be computed as $2~\rm{min}[P(\beta \leq 0|Y),P(\beta \geq 0|Y))]$. If all priors are flat, this Bayesian p-value is consistent with the frequentist p-value obtained from a two-sided Wald test. Moreover, as credible sets and confidence intervals of the same coverage are asymptotically equivalent, this Bayesian p-value is also asymptotically consistent with the frequentist p-value for any prior.

MODIFICATION :

01/09/2012 - Prague - Refactoring

INFORMATIONS:

Parameters:
[in,out]bInitial starting points then optimized value.
[in]moptimized parameters vector lenth
[out]niNumber of iteration during optimization
[out]vHessian and gradient at maximum
[out]rlFunction (penalized likelihood) value at maximum
[out]istopConvergence status
[out]convConvergence criteria value (difference in parameter norm, difference in function value and RDM)
[out]llikelihoodunpenalized log-likelihood
[out]LCV_aModel quality of fit assessment

Definition at line 82 of file optimization.f90.

References WorkingSharedValues::abserrfuncpa, computeMeasurmentError(), WorkingSharedValues::controleCA, WorkingSharedValues::controleCB, WorkingSharedValues::controleRDM, dchole(), derivMARC(), derivRVS(), dmaxt(), dsinv(), WorkingSharedValues::firstFuncpa, LCVa(), WorkingSharedValues::loglike, WorkingSharedValues::marqONLY, WorkingSharedValues::maxiteration, WorkingSharedValues::nbpatOK, WorkingSharedValues::nomfileOut, mpimod::numproc, WorkingSharedValues::OutputFolder, searpas(), WorkingSharedValues::SwitchMarquartd, and WorkingSharedValues::withexclusion.

Referenced by nimrod().

Here is the call graph for this function:

Here is the caller graph for this function: