[R] Quantile function for the generalized beta distribution of the 2nd kind
Peter Dalgaard
p.dalgaard at biostat.ku.dk
Mon Dec 12 11:38:30 CET 2005
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
More information about the R-help
mailing list