On Oct 18, 2013, at 08:37 , David Winsemius wrote:
>
> On Oct 17, 2013, at 9:11 PM, Steven LeBlanc wrote:
>
>> Greets,
>>
>> I'm using nlminb() to estimate the parameters of a multivariate normal
>> random sample with missing values and ran into an unexpected result from my
>> call to dmvnorm()
>
> There are at least 5 different version of dmvnorm. None of them are in the
> default packages.
>
>> within the likelihood function. Particular details are provided below.
>
> Complete? Except for the name of the package that has `dmvnorm`.
>
More importantly, it gives no clue as to the connection between sigma and the
data set. It is not the covariance matrix:
> s <- scan(what=list("",0,0))
1: [1,] 0.84761637 3.994261
2: [2,] 0.91487059 4.952595
3: [3,] 0.84527267 4.521837
....
40: [40,] 0.65938218 5.209301
41:
Read 40 records
> cor(s[[2]],s[[3]])
[1] 0.8812403
> colMeans(cbind(s[[2]],s[[3]]))
[1] 1.252108 5.540686
> var(cbind(s[[2]],s[[3]]))
[,1] [,2]
[1,] 1.284475 2.536627
[2,] 2.536627 6.450582
These are not the u and sigma stated.
Furthermore the matrix given as sigma is not a covariance matrix. Try working
out the correlation coefficient:
> 2.2289513/sqrt(0.6464647*5.697834)
[1] 1.161377
That should be enough to make any version of dmvnorm complain...
>> It appears that dmvnorm() makes a call to log(eigen(sigma)). Whereas
>> eigen(sigma) is returning a negative number, I understand log()'s complaint.
>> However, it is a mystery to me why this data set should produce such a
>> result.
>>
>> Any suggestions?
>>
>> Best Regards,
>> Steven
>>
>>> complete
>> [,1] [,2]
>> [1,] 0.84761637 3.994261
>> [2,] 0.91487059 4.952595
>> [3,] 0.84527267 4.521837
>> [4,] 2.53821358 8.374880
>> [5,] 1.16646209 6.255022
>> [6,] 0.94706527 4.169510
>> [7,] 0.48813564 3.349230
>> [8,] 3.71828469 9.441518
>> [9,] 0.08953357 1.651497
>> [10,] 0.68530515 5.498403
>> [11,] 1.52771645 8.484671
>> [12,] 1.55710697 5.231272
>> [13,] 1.89091603 4.152658
>> [14,] 1.08483541 5.401544
>> [15,] 0.58125385 5.340141
>> [16,] 0.24473250 2.965046
>> [17,] 1.59954401 8.095561
>> [18,] 1.57656436 5.335744
>> [19,] 2.73976992 8.572871
>> [20,] 0.87720252 6.067468
>> [21,] 1.18403087 3.526790
>> [22,] -1.03145244 1.776478
>> [23,] 2.88197343 7.720838
>> [24,] 0.60705218 4.406073
>> [25,] 0.58083464 3.374075
>> [26,] 0.87913427 5.247637
>> [27,] 1.10832692 3.534508
>> [28,] 2.92698371 8.682130
>> [29,] 4.04115277 11.827360
>> [30,] -0.57913297 1.476586
>> [31,] 0.84804365 7.009075
>> [32,] 0.79497940 3.671164
>> [33,] 1.58837762 5.535409
>> [34,] 0.63412821 3.932767
>> [35,] 3.14032433 9.271014
>> [36,] -0.18183869 1.666647
>> [37,] 0.57535770 6.881830
>> [38,] 3.21417723 10.901636
>> [39,] 0.29207932 4.120408
>> [40,] 0.65938218 5.209301
>>> u
>> [1] 1.267198 5.475045
>>> sigma
>> [,1] [,2]
>> [1,] 0.6461647 2.228951
>> [2,] 2.2289513 5.697834
>>> dmvnorm(x=complete,mean=u,sigma=sigma)
>
>> [1] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
>> NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
>> [30] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
>> Warning message:
>> In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
>> NaNs produced
>> ______________________________________________
>> 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.
>
> David Winsemius
> Alameda, CA, USA
>
> ______________________________________________
> 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
______________________________________________
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.
