[R] Dirichlet surface

[Kehl Dániel]

Dear list members,

I want to draw surfaces of Dirichlet distributions with different 
parameter settings.
My code is the following:
#<begin code>
a1 <- a2 <- a3 <- 2
#a2 <- .5
#a3 <- .5
x1 <- x2 <- seq(0.01, .99, by=.01)

f <- function(x1, x2){
       term1 <- gamma(a1+a2+a3)/(gamma(a1)*gamma(a2)*gamma(a3))
       term2 <- x1^(a1-1)*x2^(a2-1)*(1-x1-x2)^(a3-1)
       term3 <- (x1 + x2 < 1)
       term1*term2*term3
       }

z <- outer(x1, x2, f)
z[z<=0] <- NA

persp(x1, x2, z,
       main = "Dirichlet Distribution",
       col = "lightblue",
       theta = 50,
       phi = 20,
       r = 50,
       d = 0.1,
       expand = 0.5,
       ltheta = 90,
       lphi = 180,
       shade = 0.75,
       ticktype = "detailed",
       nticks = 5)
#<end code>

It works fine (I guess), except for a1=a2=a3=1. In that case I get the 
error: Error in persp.default...  invalid 'z' limits.
The z matrix has only elements 2 and NA.

Any ideas are appreciated.

Thank you:
Daniel
University of Pécs