SciencesPo (version 1.02.12)

rDirichlet: Samples from Dirichlet distribution

Description

Generates random deviates from the Dirichlet distribution. This code was originally posted by Ben Bolker in the R-News on Fri Dec 15 2000. But Ben attributed the code to Ian Wilson.

Usage

rDirichlet(n, alpha)

Arguments

n

Number of random vectors to generate

alpha

Vector containing shape parameters

Value

• Returns a matrix with n rows, each containing a single Dirichlet random deviate.

encoding

UTF-8

Examples

# 1 - Simple usage
rDirichlet(20, c(1,1,1) )
# 2 -
alpha = c( 5.0, 1.0, 2.0 )
alpha.0 = sum( alpha )
test = rDirichlet( 100000, alpha )
apply( test, 2, mean )
alpha / alpha.0
apply( test, 2, var )
alpha * ( alpha.0 - alpha ) / ( alpha.0^2 * ( alpha.0 + 1 ) )
# 3 - Brazil poll, by Datafolha
## Face-to-face interviews conducted on Oct 03-04 with n = 18116
n <- 18116
poll <- c(40,24,22,5,5,4) / 100 * n # data
set.seed(1234)
mcmc <- 100000
sim <- rDirichlet(mcmc, alpha = poll + 1)
margin <- sim[,2] - sim[,3]
mn <- mean(margin) # Bayes estimate
mn
s <- sd(margin) # posterior standard deviation

qnts <- quantile(margin, probs = c(0.025, 0.975)) # 90\% credible interval
qnts
pr <- mean(margin > 0) # posterior probability of a positive margin
pr
hist(margin, prob = TRUE, # posterior distribution
    breaks = "FD", xlab = expression(p[2] - p[3]),
    main = expression(paste(bold("Posterior distribution of "), p[2] - p[3])))
abline(v=mn, col='red', lwd=3, lty=3)