R help - Mar 2011 - Dirichlet surface

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

On Mar 29, 2011, at 4:45 PM, Kehl D?niel wrote:
> 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.
Might be a buglet in persp. If you set the zlim argument to c(0,2.1),  
you get the expected constant plane at z=2 for those
arguments.
>
> Any ideas are appreciated.
>
> Thank you:
> Daniel
> University of P?cs
>
> ______________________________________________
> R-help at r-project.org 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.
David Winsemius, MD
Heritage Laboratories
West Hartford, CT