Probability density function for, and random sampling from, the hyperdirichlet distribution

Share:

Description

Probability density function for the hyperdirichlet distribution in terms of either p or e; and random sampling using Metropolis-Hastings

Usage

1
2
3
dhyperdirichlet_e(e, HD, include.Jacobian = TRUE)
dhyperdirichlet(p, HD, include.NC = FALSE, TINY = 1e-10, log = FALSE)
rhyperdirichlet(n, HD, start=NULL, sigma=NULL)

Arguments

HD

Object of class hyperdirichlet, or coerced thereto

p

Vector of length dim(HD), notionally summing to one

e

Vector of length dim(HD) giving the point in e-space

include.Jacobian

In function dhyperdiriclet_e(), Boolean with default TRUE meaning to include the Jabobian of the transform from e to p

include.NC

In function dhyperdirichlet_e(), Boolean with TRUE meaning to include the normalization factor and default FALSE meaning not to include it (it is expensive to calculate). Note that if the normalizing factor is not known, the function will return NA

TINY

In function dhyperdirichlet_p(), numeric, specifying minimum size for elements of p via p <- pmax(p , TINY)

log

In function dhyperdirichlet_p(), Boolean with default FALSE meaning to return the probability density and TRUE meaning to return its logarithm

n,start,sigma

In function rhyperdirichlet(), n is the number of observations to take, start is the start-point for the random walk (with default NULL meaning to use the neutral point), and sigmais the standard deviation for the (Gaussian) kernel, with default NULL meaning to use 1/d

Details

Function dhyperdirichlet() gives the density as a function of the p_1, ..., p_d.

Function dhyperdirichlet_e() gives the density as a function of the e_i. This is useful when integrating as the simplex (in p-space) transforms to a hypercube in e-space.

Value

Functions dhyperdirichlet() and dhyperdirichlet_e() return a scalar; function rhyperdirichlet() returns a matrix whose rows are k-tuples

Note

Function rhyperdirichlet() uses a Metropolis-Hastings algorithm to construct a Markov chain. Note that successive observations are not independent of one another. The details of this nonindependence are poorly understood (by me, at any rate). The examples section below shows how to generate a single random observation from the hyperdirichlet distribution.

Function dhyperdirichlet() silently normalizes p by p <- p/sum(p). Currently, no check for the elements being positive is made (negative elements are nullified by TINY).

The relationship between e and p is given in e_to_p.Rd.

Author(s)

Robin K. S. Hankin

See Also

maximum.likelihood,e_to_p

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
dhyperdirichlet(c(1,4,3,2)/10, dirichlet(1:4))

rhyperdirichlet(20, dirichlet(1:3))

diff(c(0,sort(runif(9)),1))  # random sample drawn from dirichlet(rep(1,10))

# how to draw a single observation from the hyperdirichlet:
f <- function(HD,n,...){rhyperdirichlet(n,HD=HD,...)[n,,drop=TRUE]}
f(dirichlet(1:3),n=100)
f(dirichlet(1:3),n=100)
f(dirichlet(1:3),n=100)
# Note, n=100 might not be enough burn-in.

Want to suggest features or report bugs for rdrr.io? Use the GitHub issue tracker.