See in line below.
On 08/06/13 00:00, Hertzog Gladys wrote:
Dear all,
I’m new in R and I’m trying to estimate the covariance matrix of a bivariate
normal distribution by maximizing the log-likelihood. The maximization really
has to be performed with the non-linear minimization routine (nlm).
Why?
The 2 means of the distribution are known and equal to 1.
But you simulate your data using means 2.4 and 1.3.
I’ve already tried 2 different ways to compute this covariance but for each of
them I obtained a lot of warnings and illogical values for the covariance
matrix.
In the first one, I defined the bivariate normal distribution with the command dmvnorm:
x<-rnorm(6000, 2.4, 0.6)
x <- matrix(c(x), ncol=1)
Why do you use matrix() and c() here? Totally unnecessary
(and confusing).
y<-rlnorm(6000, 1.3,0.1)
y <- matrix(c(y), ncol=1)
XY <- cbind(x,y)
To estimate the covariance matrix underlying the data XY you could
simply do:
M <- var(XY).
If you want to impose the constraint that the true mean vector is c(1,1)
you could do:
mu <- c(1,1)
xbar <- apply(XY,2,mean)
M <- (5999/6000)*var(XY) + (xbar-mu)%*%t(xbar-mu)
It would of course make more sense in terms of your simulated data
to use:
mu <- c(2.4,1.3)
I haven't looked at your code below any further. Too much of a mess.
cheers,
Rolf Turner
L <- function(par,x,y) {
return (-sum(log(par[4]*dmvnorm(XY, mean=c(1,1), sigma= matrix(c(par[1],
par[1]*par[2]*par[3],par[1]*par[2]*par[3], par[2] ),nrow=2, ncol=2))
)))
}
par.start<- c(0.5, 0.5 ,0.5 ,0.5)
result<-nlm(L,par.start,y=y,x=x, hessian=TRUE)
par.hat <- result$estimate
par.hatIl y a eu 32 avis (utilisez warnings() pour les visionner)
par.hat <- result$estimate
par.hat
[1] 5.149919e+01 2.520721e+02 8.734212e-03 3.996771e+02
warnings()
Messages d'avis :
1: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
production de NaN
2: In nlm(L, par.start, y = y, x = x, hessian = TRUE) :
NA/Inf replaced by maximum positive value
3: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
production de NaN
4: In nlm(L, par.start, y = y, x = x, hessian = TRUE) :
NA/Inf replaced by maximum positive value
5: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
production de NaN
6: In nlm(L, par.start, y = y, x = x, hessian = TRUE) :
NA/Inf replaced by maximum positive value
7: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
production de NaN
8: In nlm(L, par.start, y = y, x = x, hessian = TRUE) :
NA/Inf replaced by maximum positive value
9: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
production de NaN
10: In nlm(L, par.start, y = y, x = x, hessian = TRUE) :
NA/Inf replaced by maximum positive value
11: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
production de NaN
12: In nlm(L, par.start, y = y, x = x, hessian = TRUE) :
NA/Inf replaced by maximum positive value
13: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
production de NaN
In the second one, I wrote step by step the bivariate normal distribution in order to have each parameter separately (not in a matrix) but it didn’t work as well:
x<-rnorm(6000, 2.4, 0.6)
y<-rlnorm(6000, 1.3,0.1)
L <- function(par,x,y) {
return (-sum(log((1-par[4])*( (1/(2*pi*par[1]*par[2]*sqrt(1-par[3])))*exp(
(-1/2*(1-par[3]^2))* ((y-1)/par[2])^2 +((x-1)/par[1])^2 -
2*(y-1)*(x-1)/(par[2]*par[1]) )) )))
}
#par [1]= sigma_x , par [2]= sigma_y par [3]= rho_xy par[4] is a mixing
parameter. The final step of my calculation will be to have a mixture of
bivariate normal and log-normal distributions.
par.start<- c(0.5, 0.5 ,0.5 ,0.5)
result<-nlm(L,par.start,y=y,x=x, hessian=T)
par.hat <- result$estimate
par.ha
When I run this script, I get always 50 advices like those below:
Messages d'avis :
1: In sqrt(1 - par[3]) : production de NaN
2: In nlm(L, par.start, y = y, x = x, hessian = T) :
NA/Inf replaced by maximum positive value
3: In sqrt(1 - par[3]) : production de NaN
4: In nlm(L, par.start, y = y, x = x, hessian = T) :
NA/Inf replaced by maximum positive value
5: In sqrt(1 - par[3]) : production de NaN
6: In nlm(L, par.start, y = y, x = x, hessian = T) :
NA/Inf replaced by maximum positive value
7: In log((1 - par[4]) * ((1/(2 * pi * par[1] * par[2] * ... : production de
NaN
8: In nlm(L, par.start, y = y, x = x, hessian = T) :
NA/Inf replaced by maximum positive value
9: In log((1 - par[4]) * ((1/(2 * pi * par[1] * par[2] * ... : production de
NaN
10: In nlm(L, par.start, y = y, x = x, hessian = T) :
NA/Inf replaced by maximum positive value
11: In log((1 - par[4]) * ((1/(2 * pi * par[1] * par[2] * ... : production de
NaN
12: In nlm(L, par.start, y = y, x = x, hessian = T) :
NA/Inf replaced by maximum positive value
……
Does one of you know how to use the nlm method to estimate the covariance
matrix (and mixing parameter ) of a bivariate normal distribution?
Thank you in advance for your help and answers.
______________________________________________
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.
