dirichlet.RdThe Dirichlet distribution in likelihood (for p) form, including the generalized Dirichlet distribution due to Connor and Mosimann
dirichlet(powers, alpha)
dirichlet3(powers, lambda=NULL)
GD(alpha, beta, beta0=0)
GD_wong(alpha, beta)
rdirichlet(n,H)
is.dirichlet(H)
rp_unif(n,H)
multinom(x)In functions dirichlet() and multinom(),
a (named) vector of powers
A vector of parameters for the Dirichlet or generalized Dirichlet distribution
In function GD(), an arbitrary parameter
Object of class hyper2
Vector of weights in dirichlet3()
Number of observations
These functions are really convenience functions.
Function rdirichlet() returns random samples drawn from a
Dirichlet distribution using the gamma distribution. If second
argument H is a hyper2 object, it is tested [with
is.dirichlet()] for being a Dirichlet distribution. If so,
samples from it are returned. If not, (e.g. icons), an error
is given. If H is not a hyper2 object, it is
interpreted as a (possibly named) vector of parameters
\(\alpha\) [not a vector of powers].
Function rp_unif() returns uniformly distributed vectors,
effectively using H*0; but note that this uses Dirichlet
sampling which is much faster and better than the Metropolis-Hastings
functionality documented at rp.Rd.
Functions GD() and GD_wong() return a likelihood
function corresponding to the Generalized Dirichlet distribution as
presented by Connor and Mosimann, and Wong, respectively. In
GD_wong(), alpha and beta must be named vectors;
the names of alpha give the names of
\(x_1,\ldots,x_k\) and the last element of beta
gives the name of \(x_{k+1}\).
Function dirichlet3() returns a hyper3 object with
weights lambda. If lambda is length less than that of
powers, it is padded with 1s [so default NULL
corresponds to unit weights, that is, a hyper2 object]. A
use-case is given in inst/rock_paper_scissors_monster.Rmd.
Experimental function multinom() also takes a named vector and
interprets it in almost the same way as dirichlet(), but it
allows for repeated names. It is interpreted as the support function
for multinomial trials. If x is a named vector with no
repeats, dirichlet(x) == multinom(x). But multinom()
uses frab formalism to deal with multiple entries with the
same name, which are summed. This is the Right Way (tm) to do this.
Negative numbers, having no natural interpretation in this context,
are not allowed. The numbers do not need to be integers. One
use-case is given in the examples (which should be faster than looping
over the named vector and accumulating a hyper2 object, element
by element; but I have no timings on this as yet).
R. J. Connor and J. E. Mosimann 1969. “Concepts of independence for proportions with a generalization of the Dirichlet distribution”. Journal of the American Statistical Association, 64:194–206
T.-T. Wong 1998. “Generalized Dirichlet distribution in Bayesian Analysis”. Applied Mathematics and Computation, 97:165–181
A dirichlet distribution can have a term with zero power. But
this poses problems for hyper2 objects as zero power brackets
are dropped.
Function dirichlet3() is a replacement for now removed function
pair3().
Function rdirichlet() commits a very mild (but necessary in the
absence of a working dismat package) violation of disordR
discipline, as the columns of the returned matrix have the same order
as pnames(H)
x1 <- dirichlet(c(a=1,b=2,c=3))
x2 <- dirichlet(c(c=3,d=4))
x1+x2
#> log( a * (a + b + c)^-6 * b^2 * c^6 * (c + d)^-7 * d^4)
multinom(c(x=2,a=3,b=3,x=1)) # NB repeated 'x'
#> log(a^3 * (a + b + x)^-9 * b^3 * x^3)
x <- setNames(rpois(40,2.2), sample(letters[1:4],40,replace=TRUE))
multinom(x)
#> log(a^21 * (a + b + c + d)^-102 * b^26 * c^27 * d^28)
H <- dirichlet(c(a=1,b=2,c=3,d=4))
rdirichlet(10,H)
#> a b c d
#> [1,] 0.01489053 0.2773005 0.17545514 0.5323539
#> [2,] 0.16490358 0.4473492 0.21954889 0.1681984
#> [3,] 0.10789406 0.5322661 0.08322359 0.2766163
#> [4,] 0.16966098 0.1607731 0.24824807 0.4213179
#> [5,] 0.19018396 0.2876375 0.13460494 0.3875736
#> [6,] 0.10645204 0.4885826 0.23210971 0.1728557
#> [7,] 0.06030589 0.3329259 0.21613456 0.3906336
#> [8,] 0.14609874 0.1888472 0.25023030 0.4148238
#> [9,] 0.07240817 0.4974437 0.08988788 0.3402603
#> [10,] 0.24799704 0.3797639 0.25643394 0.1158052
colMeans(rdirichlet(1e4,H))
#> a b c d
#> 0.1443673 0.2119955 0.2860204 0.3576168
dirichlet3(c(fish=3,chips=2),lambda=1.8)
#> log( (chips=1)^2 * (chips=1, fish=1.8)^-5 * (fish=1.8)^3)
dirichlet3(c(x=6,y=5,z=2),1:3)
#> log( (x=1)^6 * (x=1, y=2, z=3)^-13 * (y=2)^5 * (z=3)^2)