Spreadshit addiction [General Statistics]
Dear HS and ElMaestro, thanks for ElMaestro's sharing of that great paper and HS's useful links for R.
To improve the precision of Dr. Russel Lenth's
I have tested the
I found a solution to improve the precision of Excel's
How to quickly leave a blank at the head of each row in BEBAC? My code is with blank at several rows, but after I submit the post, the blanks disappeared. *) See note. [HS]
Using the improved BetaInverse(), I get:
It can be seen that my BetaInverse(0.001,2,4) is much closer to R's qbeta(0.001,2,4) 0.01010178788373775407572 than Excel's BetaInv(0.001,2,4).
Thus, the new results for ElMaestro's df 5.123 example were refreshed as:
Excel's calculation is performed with improved NctInv() containing improved BetaInverse()
Now, the bias has been greatly decreased. May be acceptable in most cases.
To improve the precision of Dr. Russel Lenth's
NctInv()
in Excel, I have to find a solution to improve the Excel's BetaInv()
function instead of direct using of BetaInv()
in MS Excel.I have tested the
BetaInv(0.001,2,4)
in MS Excel 2007, it returns 0.010101318359375
which is identical to that in Excel 2003. How disappointing with MS's so called improvment. Actually, they did nothing! I will test it in MS Excel 2010 if possible.I found a solution to improve the precision of Excel's
BetaInv()
function to a level of MS's so called 10^(-15). It is based on Excel's built-in BetaDist()
function which has obviousely higher precision than BetaInv()
.'-----------------------------------------------------------------
Public Function NCtInv(p As Double, df As Double) As Double
Dim X As Double, y As Double, Incr As Double
If (p <= 0) Or (p >= 1) Or (df <= 0) Then
MsgBox "Illegal argument", vbCritical, "NCtInv"
Exit Function
End If
If p = 0.5 Then p = 0.5000000001
If p < 0.5 Then
'X = Application.BetaInv(1 - 2 * p, 0.5, df / 2)
X = BetaInverse(1 - 2 * p, 0.5, df / 2)
NCtInv = -Sqr(df * X / (1 - X))
Else
'X = Application.BetaInv(1 - 2 * (1 - p), 0.5, df / 2)
X = BetaInverse(1 - 2 * (1 - p), 0.5, df / 2)
NCtInv = Sqr(df * X / (1 - X))
End If
End Function
'-----------------------------------------------------------------
'-----------------------------------------------------------------
Public Function BetaInverse(p As Variant, Alpha As Variant, Beta As Variant) As Variant
Dim X As Variant
Dim A As Variant
Dim B As Variant
Dim Precision As Variant
X = 0
A = 0
B = 1
Precison = 10 ^ (-15)
Do While B - A > Precision
X = (A + B) / 2
If Application.BetaDist(X, Alpha, Beta) > p Then
B = X
Else
A = X
End If
Loop
BetaInverse = X
End Function
'-----------------------------------------------------------------
How to quickly leave a blank at the head of each row in BEBAC? My code is with blank at several rows, but after I submit the post, the blanks disappeared. *) See note. [HS]
Using the improved BetaInverse(), I get:
Excel's BetaInv(0.001,2,4)= 0.010101318359375
My BetaInverse(0.001,2,4)= 0.0101017878834719
It can be seen that my BetaInverse(0.001,2,4) is much closer to R's qbeta(0.001,2,4) 0.01010178788373775407572 than Excel's BetaInv(0.001,2,4).
Thus, the new results for ElMaestro's df 5.123 example were refreshed as:
df p Excel 2003 OO Calc R
5.123 0.000 #NA #NA -Inf
5.123 0.025 -2.55212890 -2.55212890 -2.55212890
5.123 0.500 0.00000006 0.00022634 0
5.123 0.975 2.55212890 2.55212890 2.55212890
5.123 1.000 #NA #NA Inf
0.000 0.975 #NA #NA NaN
Excel's calculation is performed with improved NctInv() containing improved BetaInverse()
Now, the bias has been greatly decreased. May be acceptable in most cases.
- Multiple blanks and tabs are removed (that’s the standard behaviour of HTML). If you want to keep multiple blanks highlight the text and click on the Code button right of the text area. Suitable BBCodes will be inserted before/after the text and in HTML rendered in a monospaced font. However, tabs are still removed. For details see here. [Helmut]
Complete thread:
- Welch/Satterthwaite in practice ElMaestro 2011-09-21 17:04 [General Statistics]
- Not to round is the question d_labes 2011-09-22 09:11
- my problem ElMaestro 2011-09-22 12:34
- Ways out d_labes 2011-09-27 11:54
- my problem ElMaestro 2011-09-22 12:34
- Welch/Satterthwaite in practice yicaoting 2011-09-29 16:29
- Welch/Satterthwaite in practice ElMaestro 2011-09-29 17:39
- Spreadshit addiction Helmut 2011-09-30 16:52
- what a divine post! ElMaestro 2011-09-30 18:54
- what a divine post! Helmut 2011-09-30 19:41
- what a divine post! ElMaestro 2011-10-01 13:03
- what a divine post! Helmut 2011-10-01 13:20
- what a divine post! ElMaestro 2011-10-01 13:03
- what a divine post! Helmut 2011-09-30 19:41
- Spreadshit addiction yicaoting 2011-10-01 05:36
- Spreadshit addiction ElMaestro 2011-10-01 11:13
- GNU addiction Helmut 2011-10-02 03:20
- Spreadshit addiction Helmut 2011-10-01 13:03
- Spreadshit addictionyicaoting 2011-10-01 16:39
- Spreadshit addiction yicaoting 2011-10-03 07:23
- Spreadshit addiction Helmut 2011-10-03 15:59
- Spreadshit addiction yicaoting 2011-10-03 17:35
- Spreadshit addiction Helmut 2011-10-03 17:41
- Spreadshit addiction yicaoting 2011-10-03 17:35
- Spreadshit addiction Helmut 2011-10-03 15:59
- Spreadshit addiction ElMaestro 2011-10-01 11:13
- what a divine post! ElMaestro 2011-09-30 18:54
- Not to round is the question d_labes 2011-09-22 09:11