Input parameters in Dirichlet density function

Question

I'm trying to feed the Dirichlet density function with some sample data and estimated alpha vector using R in Python (library rpy2), but I don't understand the error that I got:

>>> q['dirichlet2b'] = DirichletReg.ddirichlet(np.asarray(mySample),np.asarray(listOfAlphas), log=False, sum_up=True)
Error in rowSums(alpha) : 'x' must be an array of at least two dimensions
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python2.7/dist-packages/rpy2/robjects/functions.py", line 82, in __call__
    return super(SignatureTranslatedFunction, self).__call__(*args, **kwargs)
  File "/usr/lib/python2.7/dist-packages/rpy2/robjects/functions.py", line 34, in __call__
    res = super(Function, self).__call__(*new_args, **new_kwargs)
rpy2.rinterface.RRuntimeError: Error in rowSums(alpha) : 'x' must be an array of at least two dimensions

The function I'm dealing with is declared as: ddirichlet(x, alpha, log = FALSE, sum.up = FALSE) . Its implementation is here.

From the documentation:

alpha the Dirichlet distribution’s parameters. Can be a vector (one set of parameters for all observations) or a matrix (a different set of parameters for each observa- tion), see details Details:

Usually, alpha is a vector thus the same parameters will be used for all observations. If alpha is a
matrix, a complete set of alpha-parameters must be supplied for each observation.

So In a test I am running, these are the values of the two variables:

>>> mySample
[0.23947368421052631, 0.29122807017543861, 0.13596491228070176, 0.24473684210526317]
>>> listOfAlphas
[0.96321625726816873, 0.010664397021223898, 0.0073408705340313662, 0.011574312890446362]

Why should I need a matrix of alphas for my observations? What exactly is needed to feed a Dirichlet density function?

Answer 1

The error appears to be caused by the fact that x is not a matrix. The reason why x needs to be a matrix, in principle is that the Dirichlet is a multivariate distribution; each of your $N$ observations has, let us say, $K$ elements, leading naturally to a representation of x as an $N \times K$ matrix. Even if you only have one observation, the function as implemented in DirichletReg is looking for a $1 \times K$ matrix (but see below for an alternative.) Here's the relevant line of code from the implementation:

  if(is.null(dim(x))) stop("x must be a matrix")

alpha does not need to be a matrix; however, if it is a matrix, it needs to be $N \times K$ as well. The latter representation allows for a Dirichlet with each observation having a different parameterization.

Try constructing x as a matrix and see if that fixes the problem.

The example below uses the ddirichlet function in the gtools package, which has the same calling structure as that in DirichletReg (for some directory-related reason I can't get DirichletReg to load in python, but I can get gtools to load):

>>> xval = robjects.FloatVector([0.3,0.5,0.4,0.3,0.3,0.2])
>>> x = robjects.r['matrix'](xval, nrow=2)
>>> alpha = robjects.FloatVector([4,4,2])
>>> print(x)
     [,1] [,2] [,3]
[1,]  0.3  0.4  0.3
[2,]  0.5  0.3  0.2

>>> print(alpha)
[1] 4 4 2

>>> gtools.ddirichlet(x,alpha)
<FloatVector - Python:0x3e5e830 / R:0x2ec7748>
[5.225472, 6.804000]

Also note that gtools's implementation of ddirichlet treats x differently than DirichletReg's implementation of ddirichlet:

if (!is.matrix(x)) 
    if (is.data.frame(x)) 
        x <- as.matrix(x)
    else x <- t(x)

It looks like you can just load a vector in, as you wished, and it will be converted to a row vector. Let's try it:

>>> x2 = robjects.FloatVector([0.3,0.4,0.3])
>>> alpha = robjects.FloatVector([3,4,3])
>>> gtools.ddirichlet(x2,alpha)
<FloatVector - Python:0x2c63a28 / R:0x266ee58>
[7.838208]

So that's another fix to your problem.