[R] Dirichlet ...
ben@zoo.ufl.edu
ben at zoo.ufl.edu
Fri Dec 15 15:12:50 CET 2000
Some code (originally contributed by Ian Wilson
<i.wilson at maths.abdn.ac.uk>
## functions for the "Dirichlet function", the multidimensional ##
generalization of the beta distribution: it's the Bayesian canonical ##
distribution for the parameter estimates of a multinomial ## distribution.
"pdirichlet" and "qdirichlet" (distribution function and quantiles) would
be more difficult because you'd first have to decide how to define the
distribution function for a multivariate distribution ... I'm sure this
could be done but I don't know how
logD <- function(a) {
sum(lgamma(a)) - lgamma(sum(a))
}
######################################################
ddirichlet<-function(x,alpha)
## probability density for the Dirichlet function, where x=vector of
## probabilities
## and (alpha-1)=vector of observed samples of each type
## ddirichlet(c(p,1-p),c(x1,x2)) == dbeta(p,x1,x2)
{
s<-sum((alpha-1)*log(x))
exp(sum(s)-logD(alpha))
}
######################################################
rdirichlet<-function(n,a)
## pick n random deviates from the Dirichlet function with shape
## parameters a
{
l<-length(a);
x<-matrix(rgamma(l*n,a),ncol=l,byrow=TRUE);
sm<-x%*%rep(1,l);
x/as.vector(sm);
}
On Thu, 14 Dec 2000 j.logsdon at lancaster.ac.uk wrote:
> Is there an implementation of the Dirichlet distribution somewhere in one
> of the libraries? Or for S-plus somewhere?
>
> John
>
>
> -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
> r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
> Send "info", "help", or "[un]subscribe"
> (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
> _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
>
--
318 Carr Hall bolker at zoo.ufl.edu
Zoology Department, University of Florida http://www.zoo.ufl.edu/bolker
Box 118525 (ph) 352-392-5697
Gainesville, FL 32611-8525 (fax) 352-392-3704
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list