Example R programs and commands

	      Visualize the Dirichlet pdf on p1+p2+p3=1
				   

# All lines preceded by the "#" character are my comments.
# All other left-justified lines are my input.
# All other indented lines are the R program output.


###### VISUALIZE DIRICHLET

# The Dirichlet density is not implemented in base R.
# It is necessary to get a contributed package.
# Several are available, here is one:

install.packages('gtools')
require('gtools')

# ...provides functions 
# ddirichlet(x,alpha)   # Density function at (x1,...,xk)
# rdirichlet(n,alpha)   # n random samples

###### Dirichlet density on 3 variables x=(x1,x2,x3), 
######  with x1+x2+x3=1, and shape parameters a=(a1,a2,a3)
#
# Generate a grid of x1,x2 values:
x1 <- seq(0,1, by=0.01)
x2 <- seq(0,1, by=0.01)

# Initialize a matrix to hold the pdf value:
z <- matrix(0, nrow=length(x1), ncol=length(x2))

# Choose some reasonable shape parameters:
alpha <- c( 3,5,12 )

# Fill z by looping over all valid (x1,x2) pairs, putting in
# x3=1-x1-x2 as the third variable in ddirichlet():
for(i in 1:length(x1)) {
   for(j in 1:length(x2) ) {
     if( x1[i]+x2[j] < 1) {
       x <- c(x1[i], x2[j], 1-x1[i]-x2[j])  # so x1+x2+x3=1     
       z[i,j] <- ddirichlet(x,alpha);
     } else { 
       z[i,j] <- 0
     }
   }
}

# View the results as a perpective plot:
persp(x1,x2,z)

# View the results as a contour plot:
contour(x1,x2,z)


########  MULTIPLE DIRICHLET PLOTS as a function of shape parameters

dirichlet3persp <- function(alpha1=2, alpha2=2, alpha3=2) {
  alpha<-c(alpha1,alpha2,alpha3)
  x1 <- seq(0,1, by=0.01)
  x2 <- seq(0,1, by=0.01)
  z <- matrix(0, nrow=length(x1), ncol=length(x2))
  # Fill z by looping over all valid (x1,x2) pairs, putting in
  # x3=1-x1-x2 as the third variable in ddirichlet():
  for(i in 1:length(x1)) {
     for(j in 1:length(x2) ) {
       if( x1[i]+x2[j] < 1) {
	 x <- c(x1[i], x2[j], 1-x1[i]-x2[j])  # so x1+x2+x3=1     
	 z[i,j] <- ddirichlet(x,alpha);
       } else { 
	 z[i,j] <- 0
       }
     }
  }
  persp(x1,x2,z)
}

dirichlet3contour <- function(alpha1=2, alpha2=2, alpha3=2) {
  alpha<-c(alpha1,alpha2,alpha3)
  x1 <- seq(0,1, by=0.01)
  x2 <- seq(0,1, by=0.01)
  z <- matrix(0, nrow=length(x1), ncol=length(x2))
  # Fill z by looping over all valid (x1,x2) pairs, putting in
  # x3=1-x1-x2 as the third variable in ddirichlet():
  for(i in 1:length(x1)) {
     for(j in 1:length(x2) ) {
       if( x1[i]+x2[j] < 1) {
	 x <- c(x1[i], x2[j], 1-x1[i]-x2[j])  # so x1+x2+x3=1     
	 z[i,j] <- ddirichlet(x,alpha);
       } else { 
	 z[i,j] <- 0
       }
     }
  }
  contour(x1,x2,z)
}

### 
# Examples of use:
dirichlet3persp(12,5,32)

dirichlet3contour(12,5,32)