[R] Quantile function for the generalized beta distribution of the 2nd kind

[Peter Dalgaard]

Florent Bresson <f_bresson at yahoo.fr> writes:

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

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. 

> Thanks for help
> 
> --- Prof Brian Ripley <ripley at stats.ox.ac.uk> a écrit
> :
> 
> > 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.
> > 
> > > but it doesn't work
> > 
> > Please see the posting guide, and tell us useful
> > information about what 
> > precisely happened.
> > 
> > On Sun, 11 Dec 2005, Florent Bresson wrote:
> > 
> > > 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 ?
> > >
> > > ______________________________________________
> > > R-help at stat.math.ethz.ch mailing list
> > > https://stat.ethz.ch/mailman/listinfo/r-help
> > > PLEASE do read the posting guide!
> > http://www.R-project.org/posting-guide.html
> > >
> > 
> > -- 
> > Brian D. Ripley,                 
> > ripley at stats.ox.ac.uk
> > Professor of Applied Statistics, 
> > http://www.stats.ox.ac.uk/~ripley/
> > University of Oxford,             Tel:  +44 1865
> > 272861 (self)
> > 1 South Parks Road,                     +44 1865
> > 272866 (PA)
> > Oxford OX1 3TG, UK                Fax:  +44 1865
> > 272595
> >
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
> 

-- 
   O__  ---- Peter Dalgaard             Øster Farimagsgade 5, Entr.B
  c/ /'_ --- Dept. of Biostatistics     PO Box 2099, 1014 Cph. K
 (*) \(*) -- University of Copenhagen   Denmark          Ph:  (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)                  FAX: (+45) 35327907