[R] Dirichlet surface

From: Kehl Dániel <kehld_at_ktk.pte.hu>
Date: Tue, 29 Mar 2011 13:45:59 -0700


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



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. Received on Tue 29 Mar 2011 - 21:05:47 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Wed 30 Mar 2011 - 12:10:25 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.

list of date sections of archive