R help - Oct 2006 - problem of mode in marginal distributions of dirichlet distribution using rdirichlet{gtools}

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]]