CODEWITHSUNDEEP
rcovariancevariancemultinomial
Zhubarb
Zhubarb

Reputation: 11905

R multinomial distribution variance

In the code below, p_hat contains the MLE's of the probabilities for X1, X2 and X3 in the given data sample. According to the multinomial distribution page on Wikipedia, the covariance matrix for the estimated probabilities is calculated as below:

set.seed(102)
X <- rmultinom(n=1, size=100, prob =c(0.1,0.3,0.6))
p_hat <- X/sum(X)

# print covariance matrix
cov_matrix <- matrix(0, nrow=length(p_hat), ncol=length(p_hat))
rownames(cov_matrix) <- c("X1","X2","X3"); colnames(cov_matrix) <- c("X1","X2","X3");
for (r in 1: length(p_hat)){
  for (c in 1: length(p_hat)){
    if(r==c){cov_matrix[r,c] <- p_hat[r] * (1-p_hat[r])}
    else{cov_matrix[r,c] <- -p_hat[r] *p_hat[c]}
  }
}

Is this implementation correct?

Is there an R function that can produce this covariance matrix given prob =c(0.1,0.3,0.6) for a multinomial distribution?

Upvotes: 1

Views: 1770

Answers (1)

Jilber Urbina
Jilber Urbina

Reputation: 61204

You can even use outer and diag to get the same result

> p <- drop(p_hat)
> variance         <-  p*(1-p)
> covariance       <- -outer(p, p)
> diag(covariance) <-  variance
> covariance
       [,1]    [,2]    [,3]
[1,]  0.090 -0.0290 -0.0610
[2,] -0.029  0.2059 -0.1769
[3,] -0.061 -0.1769  0.2379

Upvotes: 4

Related Questions