[R] Dirichlet surface

Kehl Dániel kehld at ktk.pte.hu
Wed Mar 30 18:55:03 CEST 2011


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?)

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



More information about the R-help mailing list