R многодиапазонное распределение дисперсии

Question

В приведенном ниже коде p_hat содержит MLE вероятностей для X1, X2 и X3 в данном образце данных. Согласно multinomial distribution странице Википедии, ковариационная матрица для расчетных вероятностей рассчитывается по следующей формуле:

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]} 
    } 
} 

Является ли эта реализация правильно?

Есть ли функция R, которая может производить эту матрицу ковариации, заданную prob =c(0.1,0.3,0.6) для многомерного распределения?

Answer 1

Вы можете даже использовать outer и diag, чтобы получить тот же результат

> 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