From j.maspons at creaf.uab.cat  Sun Mar 18 12:36:07 2012
From: j.maspons at creaf.uab.cat (Joan Maspons)
Date: Sun, 18 Mar 2012 12:36:07 +0100
Subject: [Rd] Beta binomial and Beta negative binomial
In-Reply-To: <CAC+N9BWHxC71r2kdt6=AgexzmvAUwgRGshK0Nd4wn9wQ5k7FKQ@mail.gmail.com>
References: <CANF72pP7nw_nuoSz+P7j-jt5qtm2ohqmyX=tnWdNL=n3LG7=2Q@mail.gmail.com>
	<76F908AE-7BB0-40B8-A198-0BEA6EB18F53@gmail.com>
	<CANF72pOnJSrAgTcu0FfHS=qD0zAwDHgGBh_qoAsinLv_4QHNaw@mail.gmail.com>
	<CAC+N9BWHxC71r2kdt6=AgexzmvAUwgRGshK0Nd4wn9wQ5k7FKQ@mail.gmail.com>
Message-ID: <CANF72pOsVL6rw554zb5dpVuenpUq4ftz_u+bEYeC3Uv5XXekCQ@mail.gmail.com>

El 18 de mar? de 2012 2:46, Tim Triche, Jr. <tim.triche at gmail.com> ha escrit:
> use the gsl package for Kummer's hypergeometric and others.

I looks nice but I'm a little bit lost. Gsl have 10 hypergeometric functions:

hyperg_0F1(c, x, give=FALSE, strict=TRUE)
*hyperg_1F1_int(m, n, x, give=FALSE, strict=TRUE)
*hyperg_1F1(a, b, x, give=FALSE, strict=TRUE)
**hyperg_U_int(m, n, x, give=FALSE, strict=TRUE)
*hyperg_U(a, b, x, give=FALSE, strict=TRUE)
hyperg_2F1(a, b, c, x, give=FALSE, strict=TRUE)
hyperg_2F1_conj(aR, aI, c, x, give=FALSE, strict=TRUE)
hyperg_2F1_renorm(a, b, c, x, give=FALSE, strict=TRUE)
hyperg_2F1_conj_renorm(aR, aI, c, x, give=FALSE, strict=TRUE)
*hyperg_2F0(a, b, x, give=FALSE, strict=TRUE)

* functions with the same number of parameters
** functions with the swame number of parameters and types (a,b,c,x<-
integer, m,n<- real)
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)
genhypergeo(c(int, int, real), c(int,real), 1)

I picked the hyperg_U_int(c(1,2,2.5),c(2,3.1),1) but this give a
vector of three numbers and a warning. I'm not mathematician and this
seems to much for me. Which function is the equivalent to
Mathematica's HypergeometricPQF [1,2]?

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

I would like to but with the dependences I don't think it will fit
R-base. Is there any package who want to include these distributions
(BB and BNB)?

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

[1] http://mathworld.wolfram.com/GeneralizedHypergeometricFunction.html
[2] http://reference.wolfram.com/mathematica/ref/BetaNegativeBinomialDistribution.html

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



[Rd] Beta binomial and Beta negative binomial
Joan Maspons