Re: [Rd] Beta binomial and Beta negative binomial

From: Joan Maspons <j.maspons_at_creaf.uab.cat>
Date: Sat, 17 Mar 2012 19:38:41 +0100

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
Received on Sat 17 Mar 2012 - 18:59:34 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 - 02:40: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