Re: [R] Dirichlet surface

From: Gavin Simpson <gavin.simpson_at_ucl.ac.uk>
Date: Wed, 30 Mar 2011 18:31:17 +0100

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
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%

______________________________________________
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 Wed 30 Mar 2011 - 17:44:00 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 - 18:20:26 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