Last chance! 50% off unlimited learning
Sale ends in
Maximum likelihood estimation of the parameters of the Dirichlet distribution
dirichlet.mle(x, weights=NULL, eps=10^(-5), convcrit=1e-05, maxit=1000,
oldfac=.3, progress=FALSE)
Data frame with
Optional vector of frequency weights
Tolerance number which is added to prevent from logarithms of zero
Convergence criterion
Maximum number of iterations
Convergence acceleration factor. It must be a parameter between 0 and 1.
Display iteration progress?
A list with following entries
Vector of
The concentration parameter
Vector of proportions
Minka, T. P. (2012). Estimating a Dirichlet distribution. Technical Report.
For simulating Dirichlet vectors with matrix-wise
dirichlet.simul
.
For a variety of functions concerning the Dirichlet distribution see the DirichletReg package.
# NOT RUN {
#############################################################################
# EXAMPLE 1: Simulate and estimate Dirichlet distribution
#############################################################################
# (1) simulate data
set.seed(789)
N <- 200
probs <- c(.5, .3, .2 )
alpha0 <- .5
alpha <- alpha0*probs
alpha <- matrix( alpha, nrow=N, ncol=length(alpha), byrow=TRUE )
x <- sirt::dirichlet.simul( alpha )
# (2) estimate Dirichlet parameters
dirichlet.mle(x)
## $alpha
## [1] 0.24507708 0.14470944 0.09590745
## $alpha0
## [1] 0.485694
## $xsi
## [1] 0.5045916 0.2979437 0.1974648
# }
# NOT RUN {
#############################################################################
# EXAMPLE 2: Fitting Dirichlet distribution with frequency weights
#############################################################################
# define observed data
x <- scan( nlines=1)
1 0 0 1 .5 .5
x <- matrix( x, nrow=3, ncol=2, byrow=TRUE)
# transform observations x into (0,1)
eps <- .01
x <- ( x + eps ) / ( 1 + 2 * eps )
# compare results with likelihood fitting package maxLik
miceadds::library_install("maxLik")
# define likelihood function
dirichlet.ll <- function(param) {
ll <- sum( weights * log( ddirichlet( x, param ) ) )
ll
}
#*** weights 10-10-1
weights <- c(10, 10, 1 )
mod1a <- sirt::dirichlet.mle( x, weights=weights )
mod1a
# estimation in maxLik
mod1b <- maxLik::maxLik(loglik, start=c(.5,.5))
print( mod1b )
coef( mod1b )
#*** weights 10-10-10
weights <- c(10, 10, 10 )
mod2a <- sirt::dirichlet.mle( x, weights=weights )
mod2a
# estimation in maxLik
mod2b <- maxLik::maxLik(loglik, start=c(.5,.5))
print( mod2b )
coef( mod2b )
#*** weights 30-10-2
weights <- c(30, 10, 2 )
mod3a <- sirt::dirichlet.mle( x, weights=weights )
mod3a
# estimation in maxLik
mod3b <- maxLik::maxLik(loglik, start=c(.25,.25))
print( mod3b )
coef( mod3b )
# }
Run the code above in your browser using DataLab