Folders and files
| Name | Name | Last commit date | ||
|---|---|---|---|---|
Repository files navigation
ALTERNATIVES TO THE EM ALGORITHN FOR ML-ESTIMATION OF LOCATION, SCATTER MATRIX AND DEGREE OF FREEDOM OF THE STUDENT-T DISTRIBUTION This code belongs to the paper [1]. It is available at https://doi.org/10.1007/s11075-020-00959-w An Arxiv-Preprint is avialble at https://arxiv.org/abs/1910.06623 Please cite [1], if you use this code. The code contains an implementation of the accelerations of the EM algorithm from [1]. Further, the SQUAREM [4] and DAAREM [2] acceleration are applied to the algorithms and an application for noise level estimation is provided. For more details see [1]. For the noise level estimation, we use functions to detect homogenous areas. The code of these functions is part of the implementation of [3]. The original code for the homogenous area detection is availible at https://github.com/csutour/RNLF. Note that the functions might be modified. The code for the homogenous area detection can be found in the directory hom_area_detection and is protected by the CeCILL Licence, which is also located in the directory hom_area_detection. For questions and bug reports, please contact Johannes Hertrich (j.hertrich(at)math.tu-berlin.de). CONTENTS: 1. REQUIREMENTS 2. USAGE 3. CLASSES AND FUNCTIONS 4. EXAMPLES 5. REFERENCES --------------- 1. REQUIREMENTS ------------------------ This code was written and testet with Matlab R2019a. For some of the examples the Matlab coder and gcc are required. --------------- 2. USAGE ------------------------------- To compute the ML-estimator of samples X of size d in a dxn matrix, weights w in a 1xn matrix and a method out of 'EM', 'aEM', 'MMF', 'GMMF' run the command iterate_studentT(X,w,method,anz_steps); Here anz_steps is the maximal number of steps. To incorporate the SQUAREM acceleration use squarem_studentT(X,w,method,anz_steps); To incorporate the DAAREM acceleration use daarem_studentT(X,w,method,anz_steps); A detailed description of all functions and their corresponding input and output parameters can be found in Section 3. --------------- 3. FUNCTIONS --------------------------- In the following we describe all functions in the directory algorithms and specify the corresponding input and output parameters. The same description can be found in the headers of the functions in the code. - daarem_studentT(X,w,step_algorithm,anz_steps,stop,abs_criteria,regularize,save_obj): This function implements the algorithms - MMF - GMMF - EM algorithm - accelerated EM algorithm - ECME algorithm for Student-t distributions with DAAREM acceleration. Inputs: X - d x n array containing the samples w - 1 x n array containing the weights of the samples step_algorithm- Choice of the algorithm, which determines the nu, mu and sigma step. step_algorithm has to be of type string. For step_algorithm='MMF' the function performs the MMF. For step_algorithm='GMMF' the function performs the GMMF. For step_algorithm='aEM' the function performs the accelerated EM algorith, For step_algorithm='ECME' the function perfors the ECME algorithm. Alternatively the step_algorithm can be declared by a function handle. anz_steps - Maximum number of iterations stop - stop==1 means we apply the stopping criteria, stop==0 means we perform exactly anz_steps spteps. If nargin<5, stop is set to 0. abs_criteria - Replaces the relative stopping criteria by an absolute stopping criteria If nargin<6, abs_criteria is set to 0. regularize - we add in each step regularize*eye(d) to sigma to prevent that sigma becomes singular. If X contains at least n affine independent samples, this parameter should be set to 0. If nargin<7, regularize is set to 0. save_obj - save_obj==1 means we save in each step the negative log-likelihood value. Otherwise save_obj should be set to 0. Note that save_obj==1 slows down the performance of the algorithm. If nargin<8, stop is set to 0. Outputs [mu,nu,sigma,num_steps,time,objective]: mu,nu,sigma - estimated parameters num_steps - if stop==1 then num_steps contains the number of steps, until the stopping criteria is reached. Otherwise NaN is returned. time - execution time of the algorithm objective - if save_obj==1, objective(r) contains the negative log-likelihood value of nu_r,mu_r,Sigma_r. Otherwise objective(r)=0. - iterate_studentT(X,w,step_algorithm,anz_steps,stop,abs_criteria,regularize,save_obj): This function implements the algorithms - MMF - GMMF - EM algorithm - accelerated EM algorithm - ECME algorithm for Student-t distributions. Inputs: X - d x n array containing the samples w - 1 x n array containing the weights of the samples step_algorithm- Choice of the algorithm, which determines the nu, mu and sigma step. step_algorithm has to be of type string. For step_algorithm='MMF' the function performs the MMF. For step_algorithm='GMMF' the function performs the GMMF. For step_algorithm='aEM' the function performs the accelerated EM algorith, For step_algorithm='ECME' the function perfors the ECME algorithm. Alternatively the step_algorithm can be declared by a function handle. anz_steps - Maximum number of iterations stop - stop==1 means we apply the stopping criteria, stop==0 means we perform exactly anz_steps spteps. If nargin<5, stop is set to 0. abs_criteria - Replaces the relative stopping criteria by an absolute stopping criteria If nargin<6, abs_criteria is set to 0. regularize - we add in each step regularize*eye(d) to sigma to prevent that sigma becomes singular. If X contains at least n affine independent samples, this parameter should be set to 0. If nargin<7, regularize is set to 0. save_obj - save_obj==1 means we save in each step the negative log-likelihood value. Otherwise save_obj should be set to 0. Note that save_obj==1 slows down the performance of the algorithm. If nargin<8, stop is set to 0. Outputs [mu,nu,sigma,num_steps,time,objective]: mu,nu,sigma - estimated parameters num_steps - if stop==1 then num_steps contains the number of steps, until the stopping criteria is reached. Otherwise NaN is returned. time - execution time of the algorithm objective - if save_obj==1, objective(r) contains the negative log-likelihood value of nu_r,mu_r,Sigma_r. Otherwise objective(r)=0. - newton(start,f,der_f,max_steps,tol) This function implements newtons method for determining the zeros in the nu-steps. Inputs: start - initial value f - function handle of the function der_f - function handle of the derivative max_steps - maximum number of steps If nargin<4, max_steps is set to 1000 tol - stopping criteria. Stop if abs(f(val))<tol. If nargin<5 tol is set to 1e-5 Outputs [zero,evals]: zero - root of f evals - number of evaluations of f and its derivative. - nu_step_em(nu_r,delta_r,d,n): This function implements the nu-step of the EM algorithm and the accelerated EM algorithm Inputs: nu_r - current nu delta_r - 1 x n array, where delta_r(i)=(x_i-mu_r)'*sigma_r^(-1)*(x_i-mu_r) d - dimension of the data n - number of samples Outputs [nu_r_plus_one,evals]: nu_r_plus_one - next nu evals - number of evaluations of the subfunctions A and B - nu_step_gmmf(nu_r,delta_r,d,n): This function implements the nu-step of the GMMF Inputs: nu_r - current nu delta_r - 1 x n array, where delta_r(i)=(x_i-mu_r)'*sigma_r^(-1)*(x_i-mu_r) d - dimension of the data n - number of samples Outputs [nu_r_plus_one,evals]: nu_r_plus_one - next nu evals - number of evaluations of the subfunctions A and B - nu_step_mmf(nu_r,delta_r,d,n): This function implements the nu-step of the MMF Inputs: nu_r - current nu delta_r - 1 x n array, where delta_r(i)=(x_i-mu_r)'*sigma_r^(-1)*(x_i-mu_r) d - dimension of the data n - number of samples Outputs [nu_r_plus_one,evals]: nu_r_plus_one - next nu evals - number of evaluations of the subfunctions A and B - sample(n,mu,nu,sigma): This function generates samples of a multivariate Student-t distribution. Inputs: n - int. Number of samples generated. mu - dx1 matrix. Location of the samples. nu - double. Degree of freedom of the samples. sigma - dxd matrix. Scatter matrix Outputs: X - dxn matrix. Samples. - squarem_studentT(X,w,step_algorithm,anz_steps,stop,abs_criteria,regularize,save_obj): This function implements the algorithms - MMF - GMMF - EM algorithm - accelerated EM algorithm - ECME algorithm for Student-t distributions with SQUAREM acceleration. Inputs: X - d x n array containing the samples w - 1 x n array containing the weights of the samples step_algorithm- Choice of the algorithm, which determines the nu, mu and sigma step. step_algorithm has to be of type string. For step_algorithm='MMF' the function performs the MMF. For step_algorithm='GMMF' the function performs the GMMF. For step_algorithm='aEM' the function performs the accelerated EM algorith, For step_algorithm='ECME' the function perfors the ECME algorithm. Alternatively the step_algorithm can be declared by a function handle. anz_steps - Maximum number of iterations stop - stop==1 means we apply the stopping criteria, stop==0 means we perform exactly anz_steps spteps. If nargin<5, stop is set to 0. abs_criteria - Replaces the relative stopping criteria by an absolute stopping criteria If nargin<6, abs_criteria is set to 0. regularize - we add in each step regularize*eye(d) to sigma to prevent that sigma becomes singular. If X contains at least n affine independent samples, this parameter should be set to 0. If nargin<7, regularize is set to 0. save_obj - save_obj==1 means we save in each step the negative log-likelihood value. Otherwise save_obj should be set to 0. Note that save_obj==1 slows down the performance of the algorithm. If nargin<8, stop is set to 0. Outputs [mu,nu,sigma,num_steps,time,objective]: mu,nu,sigma - estimated parameters num_steps - if stop==1 then num_steps contains the number of steps, until the stopping criteria is reached. Otherwise NaN is returned. time - execution time of the algorithm objective - if save_obj==1, objective(r) contains the negative log-likelihood value of nu_r,mu_r,Sigma_r. Otherwise objective(r)=0. - studentT_step(X,w,nu_r,mu_r,sigma_r,regularize,nu_step,delta_r) This function performs one step of the MMF, GMMF or accelerated EM algorithm for estimating the parameters of a Student-t distribution. Inputs: X - d x n array containing the samples w - 1 x n array containing the weights of the samples nu_r,mu_r,sigma_r - current estimates of nu, mu and sigma regularize - we add in each step regularize*eye(d) to sigma to prevent that sigma becomes singular. If X contains at least n affine independent samples, this parameter should be set to 0. nu_step - declares the nu-step of the algorithm. Choose nu_step=@nu_step_em for the accelerated EM algorithm. Choose nu_step=@nu_step_gmmf for the GMMF. Choose nu_step=@nu_step_mmf for the MMF. delta_r - 1 x n array containing the values of delta_r. Set delta_r=nan to compute the values within this function. If nargin<7 delta_r is computed wighin this function Outputs [nu_r_plus_one,mu_r_plus_one,sigma_r_plus_one,delta_r_plus_one]: nu_r_plus_one - updated estimate of nu. mu_r_plus_one - updated estimate of mu. sigma_r_plus_one - updated estimate of sigma. delta_r_plus_one - updated values of delta. - studentT_step_em(X,w,nu_r,mu_r,sigma_r,regularize,nu_step,delta_r) This function performs one step of the EM algortihm or the ECME algorithm for estimating the parameters of a Student-t distribution. Inputs: X - d x n array containing the samples w - 1 x n array containing the weights of the samples nu_r,mu_r,sigma_r - current estimates of nu, mu and sigma regularize - we add in each step regularize*eye(d) to sigma to prevent that sigma becomes singular. If X contains at least n affine independent samples, this parameter should be set to 0. nu_step - declares the nu-step of the algorithm. Choose nu_step=@nu_step_em for the EM algorithm. Choose nu_step=@nu_step_gmmf for the ECME algorithm. delta_r - 1 x n array containing the values of delta_r. Set delta_r=nan to compute the values within this function. If nargin<7 delta_r is computed wighin this function Outputs [nu_r_plus_one,mu_r_plus_one,sigma_r_plus_one,delta_r_plus_one]: nu_r_plus_one - updated estimate of nu. mu_r_plus_one - updated estimate of mu. sigma_r_plus_one - updated estimate of sigma. delta_r_plus_one - returns NaN. --------------- 4. EXAMPLES ---------------------------- In the following we give a short description of the scripts in the main directory. Note that you can find a detailed description in [1]. - historgrams_nu: This script draws n=1000 samples from a student t distribution for a fixed nu out of [1,2,5,10,100,200], with identity scatter matrix and zero mean in two dimensions. Then it uses the GMMF to compute the ML-Estimator. This procedure is repeated 1000 times and the histograms of the outcoming nu estimates is plotted. The resulting histograms are given in Figure 3 in [1]. - Likelihood_plot: The script draws n=1000 samples from a student t distribution with fixed nu out of [1,2,5,10,100,200] and computes the ML-estimator using the EM, aEM, MMF and GMMF algorithm. Then it generates plots of the Likelihood value versus the number of steps of the corresponding algorithms. Some exemplary resulting plots are given in Figure 2 in [1]. - nu_estimation_camera: This script implements the noise level estimation as described in Section 6.3 in [1]. First the script detects homogenous areas in the noisy input image using kendalls tau test. For this we use the code of [3], see the header of the readme. On the homogenous areas we estimate the scatter and degree of freedom of the noise on the image. Then histograms of results and the homogenous areas are plotted. The results are given in Figure 4 in [1]. - nu_estimation_muehle: Same as nu_estimation_camera with another image. The results are given in Figure 5 in [1]. - simulation_study: Runtime comparison of the EM, aEM, MMF and GMMF algorithm using different choices of nu and sigma. The script generates for fixed nu out of [1,2,5,10,100] and Sigma out of [0.1 Id, Id,10 Id, [2,-1;-1,2]] n=1000 samples- Then it runs the algorithm and records execution time and number of steps of the algorithms using a relative stopping criteria, see [1]. This procedure is repeated 10000 times and the mean and standard deviation of the execution times and number of steps are computed. The results are printed into the file 'results'. - simulation_study_accelerations: Runtime comparison of the EM, aEM, MMF, GMMF and ECME algorithms with no acceleration, SQUAREM acceleration and DAAREM acceleration. The script generates for nu out of [1,2,5,10,100] n=1000 samples of a student-t distribution with zero mean and 0.1 Id scatter. Then it runs the algorithms and records execution time and number of steps of the algorithms using a relative stopping criteria, see [1]. This procedure is repeated 1000 times and the mean and standard deviation of the execution times and number of steps are computed. The results are printed into the file 'results_with_accelerations'. Exemplary results are given in Table 3 and 4 in [1]. - simulation_study_accelerations_absolute: Similar as simulations_study_acceleration with the only difference, that an absolute stopping criterion is used instead of a relative one. The results are printed into the file 'results_with_accelerations_absolute'. Exemplary results are given in Table 5 and 6 in [1]. --------------- 5. REFERENCES -------------------------- [1] M. Hasannasab, J. Hertrich, F. Laus, and G. Steidl. Alternatives to the EM algorithm for ML-estimation of location, scatter matrix and degree of freedom of the student-t distribution. Numerical Algorithms, 2020. DOI: https://doi.org/10.1007/s11075-020-00959-w ArXiv preprint available: arXiv:1910.06623. [2] N. C. Henderson and R. Varadhan. Damped Anderson acceleration with restarts and monotonicity control for accelerating EM and EM-like algorithms. Journal of Computational and Graphical Statistics, 28(4):834–846, 2019 [3] C. Sutour, C.A. Deledalle and J.F. Aujol. Estimation of the noise level function based on a nonparametric detection of homogeneous image regions. SIAM Journal on Imaging Sciences, 8(4):2622–2661, 2015. [4] R. Varadhan and C. Roland. Simple and globally convergent methods for accelerating the convergence of any EM algorithm. Scandinavian Journal of Statistics. Theory and Applications, 35(2):335–353, 2008.