[Rd] Beta binomial and Beta negative binomial
Joan Maspons
j.maspons at creaf.uab.cat
Sun Mar 18 12:36:07 CET 2012
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
More information about the R-devel
mailing list