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.masp...@creaf.uab.cat>wrote: > Hello, > > > El 16 de març de 2012 20:34, Christophe Dutang <duta...@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.masp...@creaf.uab.cat > http://www.creaf.uab.cat > > ______________________________________________ > R-devel@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@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
