rdirichlet: Dirichlet Distribution

In SciencesPo: A Tool Set for Analyzing Political Behavior Data

Description

Density function and random number generation for the Dirichlet distribution

Usage

rdirichlet(n, alpha)

Arguments

n

number of random observations to draw.

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 observation), see “Details”. If alpha is a matrix, a complete set of α-parameters must be supplied for each observation. log returns the logarithm of the densities (therefore the log-likelihood) and sum.up returns the product or sum and thereby the likelihood or log-likelihood.

Value

The rdirichlet returns a matrix with n rows, each containing a single random number according to the supplied alpha vector or matrix.

Author(s)

Daniel Marcelino, [e-mail]

Examples

# 1) General usage:
rdirichlet(20, c(1,1,1) );
alphas <- cbind(1:10, 5, 10:1);
alphas;
rdirichlet(10, alphas );
alpha.0 <- sum( alphas );
test <- rdirichlet(10, alphas );
apply( test, 2, mean );
alphas / alpha.0;
apply( test, 2, var );
alphas * ( alpha.0 - alphas ) / ( alpha.0^2 * ( alpha.0 + 1 ) );

# 2) A practical example of usage:
# A Brazilian face-to-face poll by Datafolha conducted on Oct 03-04
# with 18,116 interviews asking for their vote preferences among the
# presidential candidates.

## First, draw a sample from the posterior
set.seed(1234);
n <- 18116;
poll <- c(40,24,22,5,5,4) / 100 * n; # The data
mcmc <- 100000;
sim <- rdirichlet(mcmc, alpha = poll + 1);

## Second, look at the margins of Aecio over Marina in the very last moment of the campaign:
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;

## Third, plot the posterior density
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);

SciencesPo documentation built on May 29, 2017, 9:28 p.m.