Re: [Rd] Beta binomial and Beta negative binomial

From: Tim Triche, Jr. <tim.triche_at_gmail.com>
Date: Sat, 17 Mar 2012 18:46:21 -0700

use the gsl package for Kummer's hypergeometric and others.

you might find implementing the distributions in C or C++ worthwhile for speed.

thanks for doing this, by the way.

On Sat, Mar 17, 2012 at 11:38 AM, Joan Maspons <j.maspons_at_creaf.uab.cat>wrote:

> Hello,
>
>
> El 16 de marš de 2012 20:34, Christophe Dutang <dutangc_at_gmail.com> ha
> escrit:
> > Hi,
> >
> > Please look at the distribution task view (
> http://cran.r-project.org/web/views/Distributions.html) and the package
> gamlss.dist.
>
> Thanks for the tip. There are Beta binomial functions but they don't
> have the number of trials parameter so I supose it's a Beta Bernoulli
> distribution.
>
> >
> > Regards
> >
> > Christophe
> >
> > --
> > Christophe Dutang
> > Ph.D. student at ISFA, Lyon, France
> > website: http://dutangc.free.fr
> >
> > Le 16 mars 2012 Ó 18:41, Joan Maspons a Úcrit :
> >
> >> Hi,
> >> I need Beta binomial and Beta negative binomial functions ...
> >>
> >> Can I implement these new functions inside stats
> >> package following the
> >> same patterns as other probability distributions?
> >>
> >> Yours,
> >> --
> >> Joan Maspons
>
> I have implemented a prototype of the beta negative binomial:
>
> FindParamBetaDist<- function(mu, sigma){ #
> return(data.frame(a=shape1,b=shape2))
> # mu<- a/(a+b) [mean]
> # sigma<- ab/((a+b)^2 (a+b+1)) [variance]
> # Maxima: solve([mu= a/(a+b) , sigma= a*b/((a+b)^2 * (a+b+1))], [a,b]);
> a<- -(mu * sigma + mu^3 - mu^2) / sigma
> b<- ((mu-1) * sigma + mu^3 - 2 * mu^2 + mu) / sigma
> if (a <= 0 | b <= 0) return (NA)
> return (data.frame(a,b))
> }
>
> #Rmpfr::pochMpfr()?
> pochhammer<- function (x, n){
> return (gamma(x+n)/gamma(x))
> }
>
> # PMF:
> # P (X = x) = ((alpha)_n (n)_x (beta)_x)/(x! (alpha+beta)_n
> (n+alpha+beta)_x) | for | x>=0
> # (a)_b Pochhammer symbol
> dbetanbinom<- function(x, size, mu, sigma){
> param<- FindParamBetaDist(mu, sigma)
> if (is.na(sum(param))) return (NA) #invalid Beta parameters
> if (length(which(x<0))) res<- 0
> else
> res<- (pochhammer(param$a, size) * pochhammer(size, x) *
> pochhammer(param$b, x)
> / (factorial(x) * pochhammer(param$a + param$b, size)
> * pochhammer(size + param$a + param$b, x)))
> return (res)
> }
>
> curve(dbetanbinom(x, size=12, mu=0.75, sigma=.1), from=0, to=24, n=25,
> type="p")
>
> # CDF:
> # P (X<=x) = 1-(Gamma(n+floor(x)+1) beta(n+alpha, beta+floor(x)+1)
> # genhypergeo(1, n+floor(x)+1, beta+floor(x)+1;floor(x)+2,
> n+alpha+beta+floor(x)+1;1))
> # /(Gamma(n) beta(alpha, beta) Gamma(floor(x)+2)) | for | x>=0
> pbetanbinom<- function(q, size, mu, sigma){
> require(hypergeo)
> param<- FindParamBetaDist(mu, sigma)
> if (is.na(sum(param))) return (NA) #invalid Beta parameters
> res<- numeric(length(q))
> for (i in 1:length(q)){
> if (q[i]<0) res[i]<- 0
> else res[i]<- (1-(gamma(size+floor(q[i])+1) *
> beta(size+param$a, param$b+floor(q[i])+1)
> * genhypergeo(c(1, 1+size+floor(q[i]), 1+param$b+floor(q[i])),
> c(2+floor(q[i]),1+size+param$a+param$b+floor(q[i])), 1))
> / (beta(param$a, param$b) * gamma(size) * gamma(2+floor(q[i]))))
> }
> return (res)
> }
>
> ## genhypergeo not converge. Increase iterations or tolerance?
> pbetanbinom(0:10x, size=20, mu=0.75, sigma=0.03)
>
> I have to investigate
> http://mathworld.wolfram.com/GeneralizedHypergeometricFunction.html
> Any tip on how to solve the problem?
>
>
> --
> Joan Maspons
> CREAF (Centre de Recerca Ecol˛gica i Aplicacions Forestals)
> Universitat Aut˛noma de Barcelona, 08193 Bellaterra (Barcelona), Catalonia
> Tel +34 93 581 2915 j.maspons_at_creaf.uab.cat
> http://www.creaf.uab.cat
>
> ______________________________________________
> R-devel_at_r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

-- 
*A model is a lie that helps you see the truth.*
*
*
Howard Skipper<http://cancerres.aacrjournals.org/content/31/9/1173.full.pdf>

	[[alternative HTML version deleted]]


______________________________________________ R-devel_at_r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel

Received on Sun 18 Mar 2012 - 01:58:04 GMT

This quarter's messages: by month, or sorted: [ by date ] [ by thread ] [ by subject ] [ by author ]

All messages

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 Sun 18 Mar 2012 - 13:30:30 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel. Please read the posting guide before posting to the list.

list of date sections of archive