title: R multinomial distribution variance
body: <p>In the code below, <code>p_hat</code> contains the MLE's of the probabilities for X1, X2 and X3 in the given data sample. According to the <a href="http://en.wikipedia.org/wiki/Multinomial_distribution" rel="nofollow">multinomial distribution</a> page on Wikipedia, the covariance matrix for the estimated probabilities is calculated as below:</p>
  
  <pre><code>set.seed(102)
  X &lt;- rmultinom(n=1, size=100, prob =c(0.1,0.3,0.6))
  p_hat &lt;- X/sum(X)
  
  # print covariance matrix
  cov_matrix &lt;- matrix(0, nrow=length(p_hat), ncol=length(p_hat))
  rownames(cov_matrix) &lt;- c("X1","X2","X3"); colnames(cov_matrix) &lt;- 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] &lt;- p_hat[r] * (1-p_hat[r])}
      else{cov_matrix[r,c] &lt;- -p_hat[r] *p_hat[c]}
    }
  }
  </code></pre>
  
  <p>Is this implementation correct? </p>
  
  <p>Is there an R function that can produce this covariance matrix given <code>prob =c(0.1,0.3,0.6)</code> for a multinomial distribution?</p>