I apologise for not including a reproducible example with this query but I
hope that I can make things clear without one.
I am fitting some finite mixture models to data. Each mixture component
has p parameters (p=29 in my application) and there are q components to
the mixture. The number of data points is n ~ 1500.
I need to select a good q and I have been considering model selection
methods suggested in Chapter 6 of
@BOOK{mp01,
author = {McLachlan, G. J. and Peel, D.},
title = {Finite Mixture Models},
publisher = {Wiley},
address = {New York},
year = {2001}
}
One of these methods involves an "empirical information matrix" which is
the matrix of products of parameter scores at the observation level
evaluated at the MLE and summed over all observations. For example a
six-component mixture will have 6 - 1 + 29*6 = 179 parameters. So for
observation i I form the 179 by 179 matrix of products of scores and sum
these up over all 1500-odd observations.
Actually it is the log of the determinant of the resultant matrix that I
really need rather than the matrix itself. I am seeking advice on what may
be the best way to evaluate this log(det()). I have been encountering
problems using
determinant(SS,logarithm=TRUE)
and eigen(SS,only.values = TRUE)$values
shows some negative eigenvalues.
Advice will be gratefully received!
Murray Jorgensen
Calculating large determinants
4 messages · Murray Jorgensen, Martin Maechler, Ravi Varadhan
"MJ" == maj <maj at stats.waikato.ac.nz>
on Wed, 5 Dec 2007 14:18:23 +1300 (NZDT) writes:
MJ> I apologise for not including a reproducible example
MJ> with this query but I hope that I can make things clear
MJ> without one.
MJ> I am fitting some finite mixture models to data. Each
MJ> mixture component has p parameters (p=29 in my
MJ> application) and there are q components to the
MJ> mixture. The number of data points is n ~ 1500.
MJ> I need to select a good q and I have been considering model selection
MJ> methods suggested in Chapter 6 of
MJ> @BOOK{mp01,
MJ> author = {McLachlan, G. J. and Peel, D.},
MJ> title = {Finite Mixture Models},
MJ> publisher = {Wiley},
MJ> address = {New York},
MJ> year = {2001}
MJ> }
MJ> One of these methods involves an "empirical information
MJ> matrix" which is the matrix of products of parameter
MJ> scores at the observation level evaluated at the MLE and
MJ> summed over all observations. For example a six-component
MJ> mixture will have 6 - 1 + 29*6 = 179 parameters. So for
MJ> observation i I form the 179 by 179 matrix of products of
MJ> scores and sum these up over all 1500-odd observations.
MJ> Actually it is the log of the determinant of the resultant matrix that I
MJ> really need rather than the matrix itself. I am seeking advice on what may
MJ> be the best way to evaluate this log(det()). I have been encountering
MJ> problems using
MJ> determinant(SS,logarithm=TRUE)
MJ> and eigen(SS,only.values = TRUE)$values
MJ> shows some negative eigenvalues.
which is a problem?
In that case I guess your problem is that you want to estimate a
positive definite matrix S but your estimate S^ is not quite
positive definite.
Function posdefify() in CRAN package "sfsmisc" provides an old
cheap solution to this problem,
where nearPD() in package 'Matrix' (based on a donation from
Jens Oehlschlaegel) provides a more sophisticated algorithm for
this problem.
If you really only need the eigenvalues of the "corrected"
matrix, you might want to abbreviate the nearPD() function by
just returning the final 'd' vector of eigenvalues.
MJ> Advice will be gratefully received!
I'll be glad to hear if and how you'd use one of these two functions.
Martin Maechler, ETH Zurich
MJ> Murray Jorgensen
Hi Murray, A likely reason for the observed information matrix not to be positive definite is the inaccuracies in the numerical estimation of the scores. You might want to try more accurate methods (e.g. Richardson extrapolation) for approximating numerical derivatives, such as available in the package "numDeriv". This would likely alleviate your problem. Please disregard my advice, if you are using analytical, closed-form expressions for the scores. In this case, you might try the approaches suggested by Martin. Ravi. ---------------------------------------------------------------------------- ------- Ravi Varadhan, Ph.D. Assistant Professor, The Center on Aging and Health Division of Geriatric Medicine and Gerontology Johns Hopkins University Ph: (410) 502-2619 Fax: (410) 614-9625 Email: rvaradhan at jhmi.edu Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html ---------------------------------------------------------------------------- -------- -----Original Message----- From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf Of Martin Maechler Sent: Wednesday, December 05, 2007 4:42 AM To: maj at stats.waikato.ac.nz Cc: r-help at stat.math.ethz.ch Subject: Re: [R] Calculating large determinants
"MJ" == maj <maj at stats.waikato.ac.nz>
on Wed, 5 Dec 2007 14:18:23 +1300 (NZDT) writes:
MJ> I apologise for not including a reproducible example
MJ> with this query but I hope that I can make things clear
MJ> without one.
MJ> I am fitting some finite mixture models to data. Each
MJ> mixture component has p parameters (p=29 in my
MJ> application) and there are q components to the
MJ> mixture. The number of data points is n ~ 1500.
MJ> I need to select a good q and I have been considering model
selection
MJ> methods suggested in Chapter 6 of
MJ> @BOOK{mp01,
MJ> author = {McLachlan, G. J. and Peel, D.},
MJ> title = {Finite Mixture Models},
MJ> publisher = {Wiley},
MJ> address = {New York},
MJ> year = {2001}
MJ> }
MJ> One of these methods involves an "empirical information
MJ> matrix" which is the matrix of products of parameter
MJ> scores at the observation level evaluated at the MLE and
MJ> summed over all observations. For example a six-component
MJ> mixture will have 6 - 1 + 29*6 = 179 parameters. So for
MJ> observation i I form the 179 by 179 matrix of products of
MJ> scores and sum these up over all 1500-odd observations.
MJ> Actually it is the log of the determinant of the resultant matrix
that I
MJ> really need rather than the matrix itself. I am seeking advice on
what may
MJ> be the best way to evaluate this log(det()). I have been
encountering
MJ> problems using
MJ> determinant(SS,logarithm=TRUE)
MJ> and eigen(SS,only.values = TRUE)$values
MJ> shows some negative eigenvalues.
which is a problem?
In that case I guess your problem is that you want to estimate a
positive definite matrix S but your estimate S^ is not quite
positive definite.
Function posdefify() in CRAN package "sfsmisc" provides an old
cheap solution to this problem,
where nearPD() in package 'Matrix' (based on a donation from
Jens Oehlschlaegel) provides a more sophisticated algorithm for
this problem.
If you really only need the eigenvalues of the "corrected"
matrix, you might want to abbreviate the nearPD() function by
just returning the final 'd' vector of eigenvalues.
MJ> Advice will be gratefully received!
I'll be glad to hear if and how you'd use one of these two functions.
Martin Maechler, ETH Zurich
MJ> Murray Jorgensen
______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Hi Murray, A likely reason for the observed information matrix not to be positive definite is the inaccuracies in the numerical estimation of the scores. You might want to try more accurate methods (e.g. Richardson extrapolation) for approximating numerical derivatives, such as available in the package "numDeriv". This would likely alleviate your problem. Please disregard my advice, if you are using analytical, closed-form expressions for the scores. In this case, you might try the approaches suggested by Martin. It is also likely that your likelihood is quite flat around the MLE, which, of course, means that you don't have enough information in your data to separate the mixtures and estimate 179 parameters. Have you tried estimating with fewer mixture components and then working your way up? Also, with this many parameters and relatively few data points, can you ensure that you have a "global" maximum? Ravi. ---------------------------------------------------------------------------- ------- Ravi Varadhan, Ph.D. Assistant Professor, The Center on Aging and Health Division of Geriatric Medicine and Gerontology Johns Hopkins University Ph: (410) 502-2619 Fax: (410) 614-9625 Email: rvaradhan at jhmi.edu Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html ---------------------------------------------------------------------------- -------- -----Original Message----- From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf Of Martin Maechler Sent: Wednesday, December 05, 2007 4:42 AM To: maj at stats.waikato.ac.nz Cc: r-help at stat.math.ethz.ch Subject: Re: [R] Calculating large determinants
"MJ" == maj <maj at stats.waikato.ac.nz>
on Wed, 5 Dec 2007 14:18:23 +1300 (NZDT) writes:
MJ> I apologise for not including a reproducible example
MJ> with this query but I hope that I can make things clear
MJ> without one.
MJ> I am fitting some finite mixture models to data. Each
MJ> mixture component has p parameters (p=29 in my
MJ> application) and there are q components to the
MJ> mixture. The number of data points is n ~ 1500.
MJ> I need to select a good q and I have been considering model
selection
MJ> methods suggested in Chapter 6 of
MJ> @BOOK{mp01,
MJ> author = {McLachlan, G. J. and Peel, D.},
MJ> title = {Finite Mixture Models},
MJ> publisher = {Wiley},
MJ> address = {New York},
MJ> year = {2001}
MJ> }
MJ> One of these methods involves an "empirical information
MJ> matrix" which is the matrix of products of parameter
MJ> scores at the observation level evaluated at the MLE and
MJ> summed over all observations. For example a six-component
MJ> mixture will have 6 - 1 + 29*6 = 179 parameters. So for
MJ> observation i I form the 179 by 179 matrix of products of
MJ> scores and sum these up over all 1500-odd observations.
MJ> Actually it is the log of the determinant of the resultant matrix
that I
MJ> really need rather than the matrix itself. I am seeking advice on
what may
MJ> be the best way to evaluate this log(det()). I have been
encountering
MJ> problems using
MJ> determinant(SS,logarithm=TRUE)
MJ> and eigen(SS,only.values = TRUE)$values
MJ> shows some negative eigenvalues.
which is a problem?
In that case I guess your problem is that you want to estimate a
positive definite matrix S but your estimate S^ is not quite
positive definite.
Function posdefify() in CRAN package "sfsmisc" provides an old
cheap solution to this problem,
where nearPD() in package 'Matrix' (based on a donation from
Jens Oehlschlaegel) provides a more sophisticated algorithm for
this problem.
If you really only need the eigenvalues of the "corrected"
matrix, you might want to abbreviate the nearPD() function by
just returning the final 'd' vector of eigenvalues.
MJ> Advice will be gratefully received!
I'll be glad to hear if and how you'd use one of these two functions.
Martin Maechler, ETH Zurich
MJ> Murray Jorgensen
______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.