[R] problem with mode of marginal distriubtion of rdirichlet{gtools}

From: hong qin <alongway_at_gmail.com>
Date: Sat 21 Oct 2006 - 14:21:24 GMT


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. Received on Sun Oct 22 04:25:35 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Sun 22 Oct 2006 - 19:30:12 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.