Re: How to derive Covariance of Dirichlet distirubtion from its variance

On Tue, Dec 19, 2023 at 12:38 PM marc nicole <[email protected]> wrote:

> Thanks, how about me exposing it as follows:
>
> import numpy as npdef covariance_dirichlet(X,alpha):
>   cov_list = []
>   sigma = np.sum(alpha)
>   for i in range(0,len(X)):
>     for j in range(i+1,len(X)):
>       if i == j:
>         iter_val_num = np.multiply(X[i],sigma)-np.power(X[i],2)#**2
>       else:
>         iter_val_num = -np.multiply(X[i],X[j])#a*b
>       iter_res = iter_val_num/(sigma**2)*(sigma+1)
>       cov_list.append(iter_res)
>   return np.reshape(np.array(cov_list),(-1,len(X)))print(covariance_dirichlet([[0.2, 0.2, 0.6,0.8],[0.2, 0.2, 0.6,0.8]],[0.4, 5, 15,11]))
>
>
> Could you provide opinion on this implementation? or how to improve it / refactor it ?
>
>
I'm not sure what you are doing with `X` there. `dirichlet` is
parameterized only by `alpha`.

>>> import numpy as np
>>> from scipy.stats import dirichlet

>>> alpha = np.array([0.4, 5, 15])

>>> X = dirichlet.rvs(alpha, size=1_000_000, random_state=1204044331)

>>> np.cov(X.T)
array([[ 0.00089554, -0.0002219 , -0.00067364],
       [-0.0002219 ,  0.00863528, -0.00841339],
       [-0.00067364, -0.00841339,  0.00908703]])

>>> dirichlet.var(alpha)
array([0.00089829, 0.00864603, 0.00909517])

>>> def cov_dirichlet(alpha):
...     alpha = np.asarray(alpha)
...     alpha0 = np.sum(alpha)
...     denom = alpha0 * alpha0 * (1 + alpha0)
...     cov = (np.diag(alpha * alpha0) - np.multiply.outer(alpha, alpha)) /
denom
...     return cov
...
>>> cov_dirichlet(alpha)
array([[ 0.00089829, -0.00022457, -0.00067372],
       [-0.00022457,  0.00864603, -0.00842146],
       [-0.00067372, -0.00842146,  0.00909517]])

-- 
Robert Kern

_______________________________________________
SciPy-Dev mailing list -- [email protected]
To unsubscribe send an email to [email protected]
https://mail.python.org/mailman3/lists/scipy-dev.python.org/
Member address: [email protected]