[R] Dirichlet surface

[David Winsemius]

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