Deux formules de calcul de la densité (pdf) d'une distribution normale bivariée retournant différents résultats

Question

Avec le code je calcule la densité d'une distribution normale bivariée. Ici j'utilise deux formules qui devraient retourner le même résultat.

La première formule utilise le dmvnorm du package mvtnorm et la deuxième formule utilise la formule de Wikipedia (https://en.wikipedia.org/wiki/Multivariate_normal_distribution). Lorsque l'écart-type des deux distributions est égal à un (la matrice de covariance n'en comporte qu'une sur la diagonale primaire), les résultats sont les mêmes. Mais lorsque vous faites varier les deux entrées de la matrice de covariance à deux ou à un tiers ... les résultats ne sont pas tous les deux identiques.

(j'espère) J'ai lu l'aide correctement et aussi ce document (https://cran.r-project.org/web/packages/mvtnorm/vignettes/MVT_Rnews.pdf).

Ici, sur stackoverflow (How to calculate multivariate normal distribution function in R), j'ai trouvé ceci parce que ma matrice de covariance est peut-être mal définie.

Mais jusqu'à présent, je ne pouvais pas trouver une réponse ...

Ma question: Pourquoi mon code retour de différents résultats lorsque l'écart-type ne correspond à un?

J'espère avoir donné assez d'informations ... mais quand quelque chose manque, veuillez commenter. Je vais éditer ma question.

Merci beaucoup d'avance!

Et maintenant mon code:

library(mvtnorm) # for loading the package if necessary 

    mu=c(0,0) 
    rho=0 
    sigma=c(1,1) # the standard deviation which should be changed to two or one third or… to see the different results 
    S=matrix(c(sigma[1],0,0,sigma[2]),ncol=2,byrow=TRUE) 

    x=rmvnorm(n=100,mean=mu,sigma=S) 
    dim(x) # for control 
    x[1:5,] # for visualization 

    # defining a function 
    Comparison=function(Points=x,mean=mu,sigma=S,quantity=4) { 
    for (i in 1:quantity) { 
      print(paste0("The ",i," random point")) 
      print(Points[i,]) 
      print("The following two results should be the same") 
      print("Result from the function 'dmvnorm' out of package 'mvtnorm'") 
      print(dmvnorm(Points[i,],mean=mu,sigma=sigma,log=FALSE)) 
      print("Result from equation out of wikipedia") 
      print(1/(2*pi*S[1,1]*S[2,2]*(1-rho^2)^(1/2))*exp((-1)/(2*(1-rho^2))*(Points[i,1]^2/S[1,1]^2+Points[i,2]^2/S[2,2]^2-(2*rho*Points[i,1]*Points[i,2])/(S[1,1]*S[2,2])))) 
      print("----") 
      print("----") 
    } # end for-loop  
    } # end function 

    # execute the function and compare the results 
    Comparison(Points=x,mean=mu,sigma=S,quantity=4) 

Answer 1

Rappelez-vous que S est la matrice de variance-covariance. La formule que vous utilisez de Wikipédia utilise l'écart-type et non la variance. Par conséquent, vous devez insérer la racine carrée des entrées diagonales dans la formule. C'est aussi la raison pour laquelle cela fonctionne quand vous choisissez 1 comme les entrées diagonales (à la fois la variance et le SD est 1).

Voir votre code modifié ci-dessous:

library(mvtnorm) # for loading the package if necessary 

mu=c(0,0) 
rho=0 
sigma=c(2,1) # the standard deviation which should be changed to two or one  third or… to see the different results 
S=matrix(c(sigma[1],0,0,sigma[2]),ncol=2,byrow=TRUE) 

x=rmvnorm(n=100,mean=mu,sigma=S) 
dim(x) # for control 
x[1:5,] # for visualization 

# defining a function 
Comparison=function(Points=x,mean=mu,sigma=S,quantity=4) { 
    for (i in 1:quantity) { 
    print(paste0("The ",i," random point")) 
    print(Points[i,]) 
    print("The following two results should be the same") 
    print("Result from the function 'dmvnorm' out of package 'mvtnorm'") 
    print(dmvnorm(Points[i,],mean=mu,sigma=sigma,log=FALSE)) 
    print("Result from equation out of wikipedia") 
    SS <- sqrt(S) 
    print(1/(2*pi*SS[1,1]*SS[2,2]*(1-rho^2)^(1/2))*exp((-1)/(2*(1-rho^2))*(Points[i,1]^2/SS[1,1]^2+Points[i,2]^2/SS[2,2]^2-(2*rho*Points[i,1]*Points[i,2])/(SS[1,1]*SS[2,2])))) 
    print("----") 
    print("----") 
    } # end for-loop  
} # end function 

# execute the function and compare the results 
Comparison(Points=x,mean=mu,sigma=S,quantity=4) 

Donc, votre commentaire lorsque vous définissez sigma est incorrect. Dans votre code, sigma est les variances, pas les écarts-types si vous jugez par la façon dont vous construisez S.