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

How you manage to write  cov_dirichlet ?
Based on which definition of the dirichlet function ? (np.diag and
np.multiply.outer got me confused) could you explain more?

also the function is supposed to take two random variables Xi and Xj to
measure their covariance not just alpha (what is alpha anyways?)

Also could alpha be multi dimensional ?

This is the update of my version without the X


import numpy as npdef covariance_dirichlet(alpha):
  cov_list = []
  sigma = np.sum(alpha)
  for i in range(0,len(alpha)):
    for j in range(i+1,len(alpha)):
      if i == j:
        iter_val_num = np.multiply(alpha[i],sigma)-np.power(alpha[i],2)#**2
      else:
        iter_val_num = -np.multiply(alpha[i],alpha[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(alpha)))print(covariance_dirichlet([[0.2,
0.2, 0.6,0.8],[0.2, 0.2, 0.6,0.8]]))


Le mar. 19 déc. 2023 à 18:53, Robert Kern <[email protected]> a écrit :

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