#############################################################################
# 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 <- 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 <- 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 <- 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 <- dirichlet.mle( x , weights= weights )
# mod3a
# # estimation in maxLik
# mod3b <- maxLik::maxLik(loglik, start=c(.25,.25))
# print( mod3b )
# coef( mod3b )
# ## End(Not run)
Run the code above in your browser using DataLab