[R] Dirichlet surface

[Gavin Simpson]

On Wed, 2011-03-30 at 09:55 -0700, Kehl Dániel wrote:
> Dear David,
> 
> I think that is a small bug too, maybe because the function is constant?
> is there a nice way to put the c(0,2.1) argument optionally, only if all 
> the parameters are 1?
> Should I post the problem somewhere else (developers maybe?)

I don't think this is a bug really; the code is having to compute limits
of the z axis and you supplied it one bit of information: a 2. If your
data are so degenerate then it is not unreasonable to expect some user
intervention. Admittedly, persp doesn't seem to work like other R
plotting functions.

You could do something like:

persp(x1, x2, z, 
      zlim = if(length(na.omit(unique(as.vector(z)))) < 2){ c(0,2.1) }
else { NULL})

in your call to persp so it only uses user-defined limits if the number
of numeric values in z is less than 2.

HTH

G

> thanks:
> Daniel
> 
> 2011-03-30 04:42 keltezéssel, David Winsemius írta:
> >
> > 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
> >
> >
> >
> 
> ______________________________________________
> 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.

-- 
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
 Dr. Gavin Simpson             [t] +44 (0)20 7679 0522
 ECRC, UCL Geography,          [f] +44 (0)20 7679 0565
 Pearson Building,             [e] gavin.simpsonATNOSPAMucl.ac.uk
 Gower Street, London          [w] http://www.ucl.ac.uk/~ucfagls/
 UK. WC1E 6BT.                 [w] http://www.freshwaters.org.uk
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%