Equivalent of beta functions

Hi!

I'd like to use the Matlab functions betapdf and betacdf, but I do not have access to the statistical toolbox of Matlab. Do you know if someone has developed an equivalent of these function, that can achieve calculation for high values of alpha and beta rather robustly?

Thanks!

Martin

Answer 1

[Youssef Khmou]

hi,

Concerning the function beta, you can use the gamma function, do you have access to the gamma function? anyway:

 BETA(x,y)=Gamma(x)*gamma(y)/gamma(x+y)

Write your own function :

 function y=myBeta(x,y)
 y=gamma(x)*gamma(y)./(gamma(y+x))

If you do not have the gamma function you can quickly write your own defined gamma using factorial or prod functions which a restriction that it will be correct only for INTEGER inputs :

 gamma(n+1)=n!               n in N

Answer 2

[Roger Stafford]

The 'betainc' function in the regular matlab package computes the same function as 'betacdf' in the Statistics Toolbox. As for 'betapdf', you can use the 'beta' function together with your own expression:

 y = 1/beta(a,b)*x^(a-1)*(1-x)^(b-1)

to achieve the equivalent of 'betapdf'.

In case it is of interest, you can even obtain the equivalent of 'betarnd' of the Statistics Toolbox using the regular package's 'betaincinv' along with the 'rand' function:

 betaincinv(rand,a,b)

[Martin]

Thank you, everyting works now! All this is a bit weird by the way, as betainc and betacdf seem to do the same thing.

Thanks again!

[Martin]

Hi Roger!

Another question, I have troubles with your formula when alpha or beta are > 1000, which I do not have with betapdf. (beta(1000,1000) = 0 but betapdf ~= 0). Would you have an idea to improve its accuracy still without using the toolbox? For the momemt I take the derivative of betainc.

Thanks again,

Martin

Answer 3

[Roger Stafford]

Use this formula for betapdf(x,a,b) to handle a larger range of a and b values:

 y = exp(-betaln(a,b)+(a-1)*log(x)+(b-1)*log(1-x));

The function 'betaln' is in the regular matlab package.

[Martin]

Thank you, it works perfectly. You know everything! :)