Hi all,
I have a problem using rdirichlet{gtools}.
For Dir( a1, a2, ..., a_n), its mode can be found at $( a_i -1)/ (
\sum_{i}a_i - n)$;
The means are $a_i / (\sum_{i} a_i ) $;
I tried to study the above properties using rdirichlet from gtools. The code
are:
##############
library(gtools)
alpha = c(1,3,9) #totoal=13
mean.expect = c(1/13, 3/13, 9/13)
mode.expect = c(0, 2/10, 8/10) # this is for the overall mode.
mode1 = numeric(3);
mode2 = numeric(3);
theta = data.frame( rdirichlet( 10000, alpha) )
m1 = mean(theta)
for( i in 1:3) {
h = hist(theta[,i], breaks=20);
mode1[i]= h$mid[ h$density == max(h$density) ]
}
theta = data.frame( rdirichlet( 10000, alpha) )
m2 = mean(theta)
for( i in 1:3) {
h = hist(theta[,i], breaks=20);
mode2[i]= h$mid[ h$density == max(h$density) ]
}
rbind(m1,m2,mean.expect) #the means are consistent
rbind(mode1, mode2, mode.expect); #the marginal modes are quite
different from the global mode.
############
An example output is:
> rbind(m1,m2,mean.expect) #the means are consistent
X1 X2 X3
m1 0.07609384 0.2301840 0.6937222
m2 0.07716923 0.2300336 0.6927971
mean.expect 0.07692308 0.2307692 0.6923077
> rbind(mode1, mode2, mode.expect);
[,1] [,2] [,3]
mode1 0.025 0.175 0.725
mode2 0.010 0.175 0.725
mode.expect 0.000 0.200 0.800
So, what is the problem with the mode?
Many thanks,
HQ
[[alternative HTML version deleted]]
______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
