On Jun 20, 2012, at 06:17 , <nata...@orchidpharma.com> <nata...@orchidpharma.com> wrote:
> Dear Peter, > > Thanks for your reply, as per my last post; my current problem is not on > positive definite related but to compare the objective function narrated by > Kjetil for optimization of parameters Vs Objective function that I have > posted in my last post. > > Though I appreciate the Kjetil approach that he uses only the covariance > matrix information to optimize the parameters of the matrix but in my case I > have other values to be used in the objective function to optimize my matrix > parameters and I can not figure out how the Kjetil approach extendable for my > specific problem and in this regard my inexperience in R as well as lack of > knowledge of linear algebra is to be blamed but I have to get out of the > problem earlier, is my current situation :-). Isn't it just a composite function issue? Got: objective(Sigma) makeSigma(param,...) ## param=c(Vs, Vn) Wanted: objective2(param, ...) Solution: objective2 <- function(param, ...) objective(makeSigma(param,...)) > > Regards, > B.Nataraj > > > > > > > -----Original Message----- > From: peter dalgaard [mailto:pda...@gmail.com] > Sent: Tuesday, June 19, 2012 8:07 PM > To: Nataraj B (ORLL-Biotech) > Cc: kjetilbrinchmannhalvor...@gmail.com; r-help@r-project.org > Subject: Re: [R] Cholesky decomposition error > > > On Jun 19, 2012, at 10:31 , <nata...@orchidpharma.com> > <nata...@orchidpharma.com> wrote: > >> Dear Kjetil, thanks for your time to detail the code. Sorry, I am still not >> able to grasp why you are creating random normal distribution for the >> covariance matrix and also why I should optimize the whole matrix value >> while I only need to optimize the two variables ?. >> >> This is my 15 x 15 covariance matrix >>> M_cov >> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] >> [,15] >> [1,] 4 3 3 2 2 1 2 2 3 2 2 2 3 >> 3 1 >> [2,] 3 5 3 2 2 1 2 2 2 2 2 2 2 >> 2 1 >> [3,] 3 3 3 2 2 1 2 2 2 2 2 2 2 >> 2 1 >> [4,] 2 2 2 3 2 1 2 3 2 1 1 1 1 >> 1 1 >> [5,] 2 2 2 2 2 1 2 2 2 1 1 1 1 >> 1 1 >> [6,] 1 1 1 1 1 1 1 1 1 1 0 1 1 >> 0 0 >> [7,] 2 2 2 2 2 1 3 2 2 1 1 1 1 >> 1 1 >> [8,] 2 2 2 3 2 1 2 4 2 1 1 1 1 >> 1 1 >> [9,] 3 2 2 2 2 1 2 2 3 1 1 1 2 >> 2 1 >> [10,] 2 2 2 1 1 1 1 1 1 2 1 2 2 >> 1 0 >> [11,] 2 2 2 1 1 0 1 1 1 1 2 1 1 >> 2 1 >> [12,] 2 2 2 1 1 1 1 1 1 2 1 3 2 >> 1 0 >> [13,] 3 2 2 1 1 1 1 1 2 2 1 2 3 >> 2 0 >> [14,] 3 2 2 1 1 0 1 1 2 1 2 1 2 >> 3 1 >> [15,] 1 1 1 1 1 0 1 1 1 0 1 0 0 >> 1 1 >> > > > Notice that that matrix has rank 9 since eigenvalues no. 10 and up are > effectively zero: > >> eigen(M_cov) > $values > [1] 2.539980e+01 5.456591e+00 3.694394e+00 2.757612e+00 1.746257e+00 > [6] 1.385264e+00 6.232218e-01 5.853408e-01 3.515233e-01 1.480189e-16 > [11] 9.492860e-17 -5.918493e-17 -1.142008e-16 -2.348974e-16 -7.429798e-16 > > $vectors > ... > > >> >> And the function to check the posdef condition is >> >>> isPosDef <- function(M) >> { if ( all(M == t(M) ) ) { # first test symmetric-ity >> if ( all(eigen(M)$values > 0) ) {TRUE} >> else {FALSE} } # else {FALSE} # not symmetric >> } >> >> And my call for optimization is >> >>> optim(c(rep(1,2)),obj,method="BFGS",hessian=TRUE,F1=F1,mu1=mu1,M_cov=M_cov) >> >> And I pass the vector F1 and mu1 used in my objective function along with >> M_cov (the covariance matrix) >> >> My objective function is >> >> >>> obj<-function(theta,F1,mu1,M_cov) { >> Vs<-theta[1] >> Vn<-theta[2] >> k<-length(F1) >> for(i in 1:15) { >> for(j in 1:15) { >> >> # to make the covariance matrix with the variables to optimize >> if(i != j) >> { >> M_cov[i,j]<-M_cov[i,j]*Vs >> } >> else >> M_cov[i,j]<-M_cov[i,j]*Vs+Vn >> } >> } >> > > ...alias > > M_cov <- Vs * M_cov > diag(M_cov) <- diag(M_cov) + Vn > > Right? > > Anyways, immediate observations are that > > (1) Vn had better be positive > (2) Vs could be negative but 25.993*Vs + Vn must be positive. > > A little thought indicates that these are also sufficient conditions for > M_cov to be pos.def. So you can replace the test for positive definiteness > with bounds on your parameters. My instinct tells me that you might not be > able to interpret the Vs < 0 case, in which case, the easy way is to just > check that Vs and Vn are positive. People often ensure this by using a > log-parametrization, Vs = exp(lVs), etc. > > -pd > > >> # here I check the posdef condition of the modified covariance matrix in >> #earlier step >> checkPD<-isPosDef(M_cov) >> if(checkPD) >> { >> e<-t(F1-mu1)%*%solve(M_cov)%*%(F1-mu1) >> logl<- -k/2*log(2*pi)-0.5*log(det(M_cov))-0.5*e >> } >> >> # negative value since the function is for maximization >> return(-logl) >> } >> >> Is there anything missing in my approach of function optimization for >> variable estimation of a multivariate normal distribution? I looking forward >> your valuable comments. >> >> Regards, >> B.Nataraj >> >> >> >> -----Original Message----- >> From: Kjetil Halvorsen [mailto:kjetilbrinchmannhalvor...@gmail.com] >> Sent: Monday, June 18, 2012 10:07 PM >> To: Nataraj B (ORLL-Biotech) >> Cc: r-help@r-project.org >> Subject: Re: [R] Cholesky decomposition error >> >> see inline! >> >> On Mon, Jun 18, 2012 at 12:38 AM, <nata...@orchidpharma.com> wrote: >>> Thanks Kjetil for your detailed code for log-Cholesky but my situation is >>> like optimization of the variables inside the matrix. >> >> But you must of course call my function makemat inside the call >> to optim() or whatever function you use for optimization! >> >> >> Here is an extended example, assuming the functions I have defined earlier: >> >> library(MASS) # for mvrnorm >> library(mvtnorm) # for dmvnorm >> >> Let us simulate a multinorma sample with expectation zero (to simplify >> the example) >> and covariance matrix Sigma: >> >>> Sigma >> [,1] [,2] [,3] >> [1,] 3 1 1 >> [2,] 1 3 1 >> [3,] 1 1 3 >>> >> >>> testdat <- mvrnorm(n=100, rep(0, 3), Sigma) >> >> We can calculate the usual maxlik estimator of Sigma: >> >>> (t(testdat) %*% testdat) /100 >> [,1] [,2] [,3] >> [1,] 3.2424618 0.9670133 1.673897 >> [2,] 0.9670133 2.9447654 1.289127 >> [3,] 1.6738970 1.2891268 3.084106 >> >>> objective <- function(Sigma) { >> sum(apply(testdat, 1, FUN=function(x) dmvnorm(x, rep(0, 3), >> Sigma, log=TRUE))) >> } >>> objective(Sigma) >> >>> start >> $diag >> [1] 0.5493061 0.4904146 0.4581454 >> >> $upper >> [1] 0.5773503 0.5773503 0.4082483 >> >>> startSigma <- c(start$diag, start$upper) >>> startSigma >> [1] 0.5493061 0.4904146 0.4581454 0.5773503 0.5773503 0.4082483 >> >>> objective(Sigma) >> [1] -571.5945 >>> objective(Sigma/3) >> [1] -699.0552 >>> objective(3*Sigma) >> [1] -638.9688 >>> >> >> >> Finally: >> >> optim(startSigma, fn=function(parvec) { >> R <- makemat(parvec[1:3], parvec[4:6]) >> S <- t(R) %*% R >> print(S) # remove this when works! >> obj <- -objective(S) >> obj >> }, control=list(REPORT=1)) >> >> # with method="BFGS" dopes not work well, because takes too long >> # steps! default Nelder-Mead is slow, but works. >> >> >> Giving the following output: >> < much removed ...> >> $par >> [1] 0.5811796 0.4969062 0.3485920 0.5469560 0.9248746 0.4896758 >> >> $value >> [1] 567.648 >> >> $counts >> function gradient >> 501 NA >> >> $convergence >> [1] 1 >> >> $message >> NULL >> >>>> res <- .Last.value >>> res$par >> [1] 0.5811796 0.4969062 0.3485920 0.5469560 0.9248746 0.4896758 >>> makemat(res$par[1:3], res$par[4:6]) >> [,1] [,2] [,3] >> [1,] 1.788146 0.546956 0.9248746 >> [2,] 0.000000 1.643628 0.4896758 >> [3,] 0.000000 0.000000 1.4170710 >>> R <-makemat(res$par[1:3], res$par[4:6]) >>> t(R) %*% R >> [,1] [,2] [,3] >> [1,] 3.1974676 0.9780374 1.653811 >> [2,] 0.9780374 3.0006752 1.310711 >> [3,] 1.6538112 1.3107108 3.103266 >>> Sigma >> [,1] [,2] [,3] >> [1,] 3 1 1 >> [2,] 1 3 1 >> [3,] 1 1 3 >> >> >> >> ¿Easy? >> >> >>> >>> Let me change your own example matrix to explain the problem. >>> >>> A martrix is [,1] [,2] [,3] >>> [1,] 3*Vs+Vn 1*Vs 1*Vs >>> [2,] 1*Vs 3*Vs+Vn 1*Vs >>> [3,] 1*Vs 1*Vs 3*Vs+Vn >> >> So make your own makemat function specific to this problem! >> which you then use inside the call to optim. >> >> kjetil >>> >>> So I have to optimize two variables "Vs" and "Vn" using the maximum >>> likelihood function, which need the determinant and inverse of the matrix A >>> to compute the value. >>> >>> Isn't it logically correct that I have to seed some initial values for the >>> Vs and Vn and use it for predicting the value of the matrix and check for >>> its positive definite condition, >> >> If you construct the matrix by above method, it will automatically be >> posdef, so no need to check! >> (except as check on programming ... ) >> >> if the condition fulfill, then I have to compute the determinant value >> and inverse of the matrix in order to use them in the maximum >> likelihood function and the optimization iteration to be carried out >> further until parameters converges. >>> >>> I sincerely believe that you have points to implements log-Cholesky in my >>> situation, but I am simply not able to see where is the Choleksy >>> decomposition to be implemented in this my workflow. >>> >>> Regards, >>> B.Nataraj >>> >>> >>> >>> >>> -----Original Message----- >>> From: Kjetil Halvorsen [mailto:kjetilbrinchmannhalvor...@gmail.com] >>> Sent: Sunday, June 17, 2012 4:10 AM >>> To: Nataraj B (ORLL-Biotech) >>> Cc: r-help@r-project.org >>> Subject: Re: [R] Cholesky decomposition error >>> >>> see below. >>> >>> On Fri, Jun 15, 2012 at 11:53 PM, <nata...@orchidpharma.com> wrote: >>>> Dear Mr.Kjetil, >>>> >>>> Thanks for your comment. You have already pointed me the article in reply >>>> to one of my earlier post to this list and I am following the paper. Now I >>>> am checking for condition for positive definiteness for original matrix >>>> using a simple script (got from earlier posting in the list) and if the >>>> condition passed then the optimization function perform decomposition of >>>> the matrix using chol() function. >>>> >>>> I believe all 5 parameterization techniques listed in the paper weigh each >>>> other based on its performance and not on accuracy, please correct me if I >>>> am wrong. Since I am novice, I just want to start with the easiest method >>>> "Cholesky" for my purpose. >>>> >>>> Are you recommending log-Cholesky function , if so could you please tell >>>> me the name of the function which does log-Cholesky decomposition of a >>>> matrix. >>> >>> I dont think there is a prewritten function for log cholesy, but it is >>> very easy to write! >>> >>>> makelogchol >>> function(A) {#A should be a square posdef matrix >>> m <- dim(A)[1] >>> uind <- upper.tri(A) >>> R <- chol(A) >>> diag <- log(diag(R)) >>> upper <- R[uind] >>> return(list(diag=diag, upper=upper)) >>> } >>> >>> and its "inverse" --- to get the cholesky factor back: >>> >>>> makemat >>> function(diag, upper) { >>> m <- length(diag) >>> mu <- m*(m-1)/2 >>> if (length(upper)!=mu) stop("incompatible lengths") >>> A <- matrix(0.0, m, m) >>> ind <- upper.tri(A) >>> A[ind] <- upper >>> diag(A) <- exp(diag) >>> A >>> } >>> >>> and an example of its use: >>> >>> Lets say A is the posdef matrix where you will start your >>> optimization, so you need >>> the log cholesy parameters to feed to optim: >>> >>>> A >>> [,1] [,2] [,3] >>> [1,] 3 1 1 >>> [2,] 1 3 1 >>> [3,] 1 1 3 >>>> start <- makelogchol(A) >>>> start >>> $diag >>> [1] 0.5493061 0.4904146 0.4581454 >>> >>> $upper >>> [1] 0.5773503 0.5773503 0.4082483 >>> >>> Then to get the cholesky factors back you call makemat: >>> >>>> R <- makemat(start$diag, start$upper) >>>> R >>> [,1] [,2] [,3] >>> [1,] 1.732051 0.5773503 0.5773503 >>> [2,] 0.000000 1.6329932 0.4082483 >>> [3,] 0.000000 0.0000000 1.5811388 >>>> t(R) %*% R >>> [,1] [,2] [,3] >>> [1,] 3 1 1 >>> [2,] 1 3 1 >>> [3,] 1 1 3 >>>> >>> >>> Kjetil >>> >>> >>> >>> >>> >>> >>> >>>> >>>> Regards, >>>> B.Nataraj >>>> >>>> >>>> >>>> -----Original Message----- >>>> From: Kjetil Halvorsen [mailto:kjetilbrinchmannhalvor...@gmail.com] >>>> Sent: Friday, June 15, 2012 11:09 PM >>>> To: Nataraj B (ORLL-Biotech) >>>> Cc: gunter.ber...@gene.com; r-help@r-project.org >>>> Subject: Re: [R] Cholesky decomposition error >>>> >>>> see inline. >>>> >>>> On Fri, Jun 15, 2012 at 4:33 AM, <nata...@orchidpharma.com> wrote: >>>>> Thanks for your reply. I am sorry and I am bit hurried up to say before >>>>> doing a proper due diligence, I have found out that during the >>>>> optimization the variables tend to vary the values of the matrix , the >>>>> function report error at some point (in particular iteration step) when >>>>> the matrix become non-decomposable due to not a positive definiteness. >>>> >>>> ¿What are you trying to optimize? If the objective function depends on >>>> a matrix argument >>>> which has to be a positive definite function, you must parametrize the >>>> matrix such that the matrix >>>> inside the optimizer always is positive definite. So if your positive >>>> definite matrix is A, then, for example, represent it as >>>> its cholesky decomposition A= L L^T where L is lower triangular with >>>> positive diagonal. Here the stricly >>>> upper diagonal part varies freely, but the diagonal not, so represent >>>> the diagonal as exp( l_i) >>>> where now the l_i varies freely. This is called the log-Cholesky >>>> parametrization. For other ideas along this lines, see the paper by >>>> Douglas Bates: "Unconstrained Parameterizations for >>>> Variance-Covariance Matrices " >>>> which you can find by googling. >>>> >>>> Kjetil >>>> >>>> >>>> This I observed when I change the maximum iteration of the optim >>>> function set to 1 and upto iteration no. 3 it runs , it stuck at >>>> iteration 4 and above. >>>>> >>>>> Now, I am trying to find ways to escalate such a condition inside the >>>>> function during the iteration process and if possible please help me to >>>>> do that. >>>>> >>>>> Regards, >>>>> B.Nataraj >>>>> >>>>> >>>>> -----Original Message----- >>>>> From: Bert Gunter [mailto:gunter.ber...@gene.com] >>>>> Sent: Friday, June 15, 2012 1:51 PM >>>>> To: Nataraj B (ORLL-Biotech) >>>>> Cc: r-help@r-project.org >>>>> Subject: Re: [R] Cholesky decomposition error >>>>> >>>>> Follow the posting guide,please: I believe at this point we need >>>>> reproducible code and your data to provide you help. See ?dput to post >>>>> your matrix. >>>>> >>>>> -- Bert >>>>> >>>>> On Thu, Jun 14, 2012 at 11:30 PM, <nata...@orchidpharma.com> wrote: >>>>>> >>>>>> Thanks for your reply. To my surprise I can find one more strange >>>>>> behavior of my 15X15 matrix "A", that is if I call the function >>>>>> chol(A) in the terminal it decompose the matrix fine without any errors >>>>>> or warnings. >>>>>> But if I call the function chol() within a function, which I have >>>>>> written in order to call the function (contains formula) for >>>>>> optimization routine "optim()" and also supplied with the same matrix >>>>>> "A" as argument, the error mentioned >>>>>> >>>>>>> Error in chol.default(M_cov) : >>>>>>> the leading minor of order 10 is not positive definite >>>>>> >>>>>> is surfaced during the function call by optim. >>>>>> >>>>>> Why the matrix fulfill the symmetric and positive definite for chol() in >>>>>> one case but fails in other case when the function chol() is called in >>>>>> other function ? >>>>>> >>>>>> I played around parameters of "optim" function but nothing seems to be >>>>>> working and I am confused and I am looking for some hints to introspect >>>>>> the problem further. >>>>>> >>>>>> Regards, >>>>>> B.Nataraj >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> -----Original Message----- >>>>>> From: Bert Gunter [mailto:gunter.ber...@gene.com] >>>>>> Sent: Thursday, June 14, 2012 6:18 PM >>>>>> To: Nataraj B (ORLL-Biotech) >>>>>> Cc: r-help@r-project.org >>>>>> Subject: Re: [R] Cholesky decomposition error >>>>>> >>>>>> Your matrix is not symmetric, positive definite. If you don't know >>>>>> what this means, you shouldn't be using chol() >>>>>> >>>>>> This may be because it isn't to begin with, or due to numerical error, >>>>>> it doesn't behave as one in the decomposition. My relative ignorance >>>>>> of numeric methods for linear algebra prevents me from saying more >>>>>> than that. >>>>>> >>>>>> -- Bert >>>>>> >>>>>> On Thu, Jun 14, 2012 at 4:23 AM, <nata...@orchidpharma.com> wrote: >>>>>>> Dear friends, >>>>>>> >>>>>>> When I do Cholesky decomposition for a 15x15 matrix using the function >>>>>>> chol(), I get the following error for which I do not understand the >>>>>>> meaning of the error >>>>>>> >>>>>>> Error in chol.default(M_cov) : >>>>>>> the leading minor of order 10 is not positive definite >>>>>>> >>>>>>> When I searched online for similar error reported earlier I could get >>>>>>> few hits but not of much help to resolve my error and one post >>>>>>> suggested to use different function called sechol() from accuracy >>>>>>> package but that did not work and it leads to different errors. So I >>>>>>> want to stick to function chol() itself. >>>>>>> >>>>>>> Could you please help me to find where things are going wrong in my >>>>>>> matrix? >>>>>>> >>>>>>> >>>>>>> Thanks and regards, >>>>>>> B.Natarj >>>>>>> >>>>>>> ______________________________________________ >>>>>>> R-help@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. >>>>>> >>>>>> >>>>>> >>>>>> -- >>>>>> >>>>>> Bert Gunter >>>>>> Genentech Nonclinical Biostatistics >>>>>> >>>>>> Internal Contact Info: >>>>>> Phone: 467-7374 >>>>>> Website: >>>>>> http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm >>>>>> >>>>>> >>>>> >>>>> >>>>> >>>>> -- >>>>> >>>>> Bert Gunter >>>>> Genentech Nonclinical Biostatistics >>>>> >>>>> Internal Contact Info: >>>>> Phone: 467-7374 >>>>> Website: >>>>> http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm >>>>> >>>>> ______________________________________________ >>>>> R-help@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. >>>> >>>> >>> >>> >> >> >> [[alternative HTML version deleted]] >> >> ______________________________________________ >> R-help@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. > > -- > Peter Dalgaard, Professor > Center for Statistics, Copenhagen Business School > Solbjerg Plads 3, 2000 Frederiksberg, Denmark > Phone: (+45)38153501 > Email: pd....@cbs.dk Priv: pda...@gmail.com > > -- Peter Dalgaard, Professor, Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd....@cbs.dk Priv: pda...@gmail.com ______________________________________________ R-help@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.
