0
$\begingroup$

I am trying to generate data from a beta-binomial distribution by inverting its cdf in R.

The code I have written to calculate the cdf seems to be working fine for most cases, but gives me values like 1.46 for the right-most support when alpha and beta are small. ( a and b are the alpha and beta parameters of the beta-binomial distribution):

cbetabin <- function(a,b,n,k){

    sum <- 0

    for (i in 0:k){
            b <- choose(n,i)*(beta(i+a,n-i+b)/beta(a,b))
            sum <- sum + b
    }

    return(sum)

}

Here's some sample output:

> res <- cbetabin(0.2,0.25,10,10)
> res
[1] 1.466977

I can't seem to figure out what's wrong here.

If it helps I also printed out the values of the pdf calculated by the function.

[1] 3.435193e-01 8.512258e-02 2.619279e-02 8.476199e-03 2.779507e-03
[6] 9.455368e-04 3.518647e-04 1.551412e-04 9.244767e-05 1.004164e-04
[11] 9.992411e-01

Thanks.

$\endgroup$

closed as off-topic by Xi'an, gung, Silverfish, Nick Cox, Sven Hohenstein Apr 2 '16 at 15:16

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question appears to be off-topic because EITHER it is not about statistics, machine learning, data analysis, data mining, or data visualization, OR it focuses on programming, debugging, or performing routine operations within a statistical computing platform. If the latter, you could try the support links we maintain." – Xi'an, Silverfish, Nick Cox, Sven Hohenstein
If this question can be reworded to fit the rules in the help center, please edit the question.

  • 3
    $\begingroup$ This appears to be a pure coding issue, which would generally be regarded as off topic here. Please see the help center in particular the section related to programming; you might be able to make it on topic at stackoverflow perhaps (see their own on-topic help), in which case edit this as needed and then flag for migration. Alternatively, if you can rephrase it to a more general question about the approach to solving the underlying statistics problem, it might be more clearly on topic here. This will likely be closed is left as it stands. $\endgroup$ – Glen_b Apr 2 '16 at 2:26
2
$\begingroup$

You have a simple implementation error. You are overwriting the input parameter b within your for loop. The following ought to do the trick:

cbetabin <- function(a,b,n,k){

    sum <- 0

    for (i in 0:k){
            summand <- choose(n,i)*(beta(i+a,n-i+b)/beta(a,b))
            sum <- sum + summand
    }

    return(sum)

}
$\endgroup$
0
$\begingroup$
b <- choose(n,i)*(beta(i+a,n-i+b)/beta(a,b))

Every time it generate a new 'b', it will use the new value of b to calculate beta() in the next loop. Don't use b, change it into any other letter, such as:

...    
t <- choose(n,i)*(beta(i+a,n-i+b)/beta(a,b))
sum <- sum + t
...
$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.