Skip to content

Quantile function for the generalized beta distribution of the 2nd kind

4 messages · Brian Ripley, Florent Bresson, Peter Dalgaard

Original
Florent Bresson
# 
I have succeded in defining the cdf of the generalized
beta of the second kind, eg.

pgbeta2 <-  function(quint,b,a,p1,p2) {
integrate(function(x)
{exp(log(a)+(a*p1-1)*log(x)-(a*p1)*log(b)-log(beta(p1,p2))-(p1+p2)*log(1+(x/b)^a))},0,quint)$value
}

but I'm facing problems with the quantile function. I
tried something like 

qgbeta2 <-  function(proba,b,a,p1,p2) {
optimize(function(z)
{(proba-pgbeta2(z,b,a,p1,p2))^2},lower=0,
upper=10^200) }

but it doesn't work. I tried with other non linear
optimization command like optim but it is apparently
not the solution.
Any idea ?
Brian Ripley
# 
Don't you want to use uniroot() to find quantiles?  It is the usual way.

Note that if you use integrate(), the result is not guaranteed to be 
smooth function of the parameters.  I may well help to decrease the 
tolerances.
Please see the posting guide, and tell us useful information about what 
precisely happened.
On Sun, 11 Dec 2005, Florent Bresson wrote:

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Florent Bresson
    

    

      
      #
    

  


  

      
The problem was with the use of the integrate command
for the definition of the cdf of the generalized beta
of the second kind. I solved the problem with a
transformation of the function qbeta. I then used an
optimize command but using uniroot is maybe  nicer. So
my function is :

qgbeta2  <-  function(proba,b,a,p1,p2) { val
<-qbeta(proba,p1,p2)   
                b*(uniroot(function(z) {z/(z+1)-val},
lower=0, upper=10000000,
tol=.Machine$double.eps^20)$root)^(1/a) }

The code is maybe not pretty but it works perfectly. I
just regret that it is not possible to fix the upper
limit of uniroot to Inf.

Thanks for help

--- Prof Brian Ripley <ripley at stats.ox.ac.uk> a ??crit
:

            

      
      
      
    

        
{exp(log(a)+(a*p1-1)*log(x)-(a*p1)*log(b)-log(beta(p1,p2))-(p1+p2)*log(1+(x/b)^a))},0,quint)$value

            

      
      
      
    

        

      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Peter Dalgaard
    

    

      
      #
    

  


  

      
Florent Bresson <f_bresson at yahoo.fr> writes:

            

      
      
      
    

        
Eh? You don't need uniroot to solve z/(z+1) == y 

Just set

 z <- y/(1-y)

In general, uniroot works by binary search, so is unhappy with
infinite-length intervals, but you can usually transform to a finite
interval.