Talk:Welch's t-test: Difference between revisions

Content added Content deleted
Line 853: Line 853:
For the record, here is a quick Stata translation of this betain implementation.
For the record, here is a quick Stata translation of this betain implementation.


The original C code can be simplified in several places using shortcuts (it could be simplified even more in C). Also, ifault is not used here: instead, the function returns a missing value if the arguments are invalid.
The original C code can be simplified in several places using shortcuts (it could be simplified even more in C). Also, ifault is not used here: instead, the function returns a missing value if the arguments are invalid. And log(beta(p,q)) is computed inside betain, rather than passed as an argument.


<lang stata>mata
<lang stata>mata
Line 860: Line 860:
Modified 28dec2017
Modified 28dec2017
Author:
Author:
Original Fortran 77 version by KL Majumder, GP Bhattacharjee.
Original FORTRAN77 version by KL Majumder, GP Bhattacharjee.
C version by John Burkardt.
C version by John Burkardt.
Stata version by Jean-Claude Arbaut
Stata version by Jean-Claude Arbaut
*/
*/
function betain(x,p,q,beta) {
function betain(x,p,q) {
acu = 1e-15

/* Check the input arguments. */
/* Check the input arguments. */
if (p<=0 | q<=0 | x<0 | x>1) return(.)
if (p<=0 | q<=0 | x<0 | x>1) return(.)

/* Special cases. */
/* Special cases. */
if (x==0 | x==1) {
if (x==0 | x==1) {
Line 875: Line 873:
}
}
acu = 1e-15
lnbeta = lngamma(p)+lngamma(q)-lngamma(p+q)

/* Change tail if necessary and determine S. */
/* Change tail if necessary and determine S. */
psq = p+q
psq = p+q
Line 906: Line 907:
if (temp<=acu & temp<=acu*value) {
if (temp<=acu & temp<=acu*value) {
value = value*exp(pp*log(xx)+(qq-1)*log(cx)-beta)/pp
value = value*exp(pp*log(xx)+(qq-1)*log(cx)-lnbeta)/pp
return(indx?1-value:value)
return(indx?1-value:value)
}
}