[WP34s] Regularized incomplete Beta function

[Dieter]

I just tried the regularized incomplete Beta function \(I\beta\) on the 34s. Looking at appendix I of the PDF manual and the given mathematical definition, the argument x and the two parameters a and b of \(I_x(a, b)\) would have to be entered this way:

Z: b
Y: a
X: x

However, this does not work - it throws errors and/or returns wrong results. After some experimenting it turned out that the actual parameter sequence apparently looks like this:

Z: x
Y: b
X: a

Example:
Evaluating \(I_{0.1}(3, 5) = 0.0256915\) requires the sequence 0.1 [ENTER] 5 [ENTER] 3.

Which leads to the question whether the manual or the implementation is wrong. Or... is it me ?-) I am using v3.2 3405. Did this get changed in the meantime?

BTW, a stack diagram would be helpful with multi-argument functions like this.

Dieter

[Thomas Klemm]

(05-01-2014 07:17 PM)Dieter Wrote:  a stack diagram would be helpful with multi-argument functions like this.

You can find them in the WP 34S pocket reference.

Cheers
Thomas

[Dieter]

(05-01-2014 07:48 PM)Thomas Klemm Wrote:  

(05-01-2014 07:17 PM)Dieter Wrote:  a stack diagram would be helpful with multi-argument functions like this.

You can find them in the WP 34S pocket reference.

Hm - while I can find the function there, the reference only says that it takes three arguments in x, y and z. Their order is the same as in the 34s manual.

Dieter

[Thomas Klemm]

(05-01-2014 08:00 PM)Dieter Wrote:  the reference only says that it takes three arguments in x, y and z

And then it shows how these arguments are used:
\[
I_\beta = \frac{\beta_x(x,y,z)}{\beta(y,z)}
\]
Comparing this to the basic definition of the regularized incomplete beta function I assume there's a typo and it should be:
\[
I_\beta = \frac{\beta_x(y,z)}{\beta(y,z)}
\]

But then why not using the following stack diagram: b a z \(\rightarrow I_z(a, b)\) ?
Or maybe change that to: z a b \(\rightarrow I_z(a, b)\)

Cheers
Thomas

[walter b]

(05-01-2014 07:17 PM)Dieter Wrote:  I just tried the regularized incomplete Beta function \(I\beta\) on the 34s. Looking at appendix I of the PDF manual and the given mathematical definition, the argument x and the two parameters a and b of \(I_x(a, b)\) would have to be entered this way:

Z: b
Y: a
X: x

However, this does not work - it throws errors and/or returns wrong results. After some experimenting it turned out that the actual parameter sequence apparently looks like this:

Z: x
Y: b
X: a

Example:
Evaluating \(I_{0.1}(3, 5) = 0.0256915\) requires the sequence 0.1 [ENTER] 5 [ENTER] 3.

Which leads to the question whether the manual or the implementation is wrong. Or... is it me ?-)

Hmmh, let me look at p. 233 (of 244): it clearly states the parameters for Iβ with the respective stack levels. When I compare this definition to Wolfram's definition (see here) I'd derive that their b is our z, their a is our y, and their z is our x. Translating their formula (6) I'd expect e.g.

Iβ(1, 2, 3) = β(2, 3) = 0.083 333
Check: β(2, 3) = Γ(2) Γ(3) / Γ(2+3) = 1! 2! / 4! = 2/24 = 1/12 q.e.d.

This turns out being not the case since Iβ(1, 2, 3) returns a domain error with build 3626. Now let's look whether there's just a permutation error in the parameter sequence of our Iβ:

• Iβ(3, 1, 2) = domain error
• Iβ(2, 3, 1) = 1
• Iβ(1, 3, 2) = domain error
• Iβ(3, 2, 1) = 1
• Iβ(2, 1, 3) = domain error
  

Hmmh, that seems to be a more severe case We'll look into it. Thanks for reporting!

d:-I

[walter b]

Maybe we also came across 'regularized' and 'incomplete'. Iβ shall be the regularized Beta function (drop incomplete!).

At least I confused those two in my previous post. Wolfram's formula (6) shall read (translated in our function names)

Iβ(1, y, z) β(y, z) = β(y, z) and thus Iβ(1, y, z) = 1.

This points to a permutation error as explained above.

d:-I

P.S. Just found Dieter was writing faster than me. Congrats!

[Dieter]

(05-01-2014 09:08 PM)walter b Wrote:  Hmmh, let me look at p. 233 (of 244): it clearly states the parameters for Iβ with the respective stack levels. When I compare this definition to Wolfram's definition (see here) I'd derive that their b is our z, their a is our y, and their z is our x.

That's what I thought as well. Until I tried. ;-)

Quote:Iβ(1, 2, 3) = β(2, 3) = 0.083 333

Hm - in my book, \(I_1(2, 3) = 1 = I_1(3, 2)\). If (!) the three parameters are entered in the right order (i.e. different from what the manual says), this is what you get:
1 [ENTER] 3 [ENTER] 2 [Iβ] returns 1.
1 [ENTER] 2 [ENTER] 3 [Iβ] also returns 1.

Quote:Hmmh, that seems to be a more severe case We'll look into it. Thanks for reporting!

Looks like a simple parameter mixup.

Dieter

[Paul Dale]

Changing the documentation will be the quickest fix here. Reordering the arguments to this function in the code is straightforward but several of the statistical distributions will also have to be updated to match.

I guess the big question is what order should the arguments be in? I don't see any benefit to reordering them any differently from what they are currently but am willing to be convinced otherwise.


- Pauli

[walter b]

(05-01-2014 09:55 PM)Paul Dale Wrote:  Changing the documentation will be the quickest fix here.

I concur. But which order do you want to see? As mentioned above, Wolfram calls it I(z; a, b) - this would correspond to Iβ(z, y, x) on the WP 34S (assuming I(z; a, b) = I(z; b, a) always).

d:-?

[Thomas Klemm]

(05-01-2014 09:55 PM)Paul Dale Wrote:  I don't see any benefit to reordering them any differently from what they are currently but am willing to be convinced otherwise.

It violates the principle of least astonishment.

With these functions the parameters are entered in the same order as they appear in the function:

COMB:
n k \(\rightarrow C(n, k)\)

PERM:
n k \(\rightarrow P(n, k)\)

/:
y x \(\rightarrow \frac{y}{x}\)

y^x:
y x \(\rightarrow y^x\)

\(\int\):
y x \(\rightarrow \int_y^x\)

P_n:
n x \(\rightarrow P_n(x)\)


With these function the parameters are entered in the reverse order:

LOG_x:
y x \(\rightarrow Log(x, y)\)

\(\gamma_{XY}\):
y x \(\rightarrow \gamma(x, y)\)

\(I\beta\):
b a z \(\rightarrow I_\beta(z, a, b)\)


IMHO that's not consistent. Of course you are free to define it the way you want. But then there are conventions. Thus we have to remember which of the functions follow it. Or then we always have to use a manual or reference guide.

Cheers
Thomas


PS: There are two exceptions: \(\rightarrow POL\) and \(\rightarrow REC\).
But I look at them as coordinate transformations that use x and y from the stack.
For the same reason the complex operations use this order.

[Thomas Klemm]

(05-01-2014 10:14 PM)walter b Wrote:  assuming I(z; a, b) = I(z; b, a) always

I_0.1(3, 5) = 0.0256915
I_0.1(5, 3) = 0.0001765

This function is not symmetric in a and b.

[Dieter]

(05-01-2014 11:01 PM)Thomas Klemm Wrote:  It violates the principle of least astonishment.
...
IMHO that's not consistent. Of course you are free to define it the way you want. But then there are conventions. Thus we have to remember which of the functions follow it. Or then we always have to use a manual or reference guide.

First of all: yes, I am a big fan of the POLA principle as well. I love things that work naturally in the way I would expect if there was no manual. On the 34s, most functions follow this paradigm. Others do not as there are good reasons for doing things differently.

The \(log_xy\) command is such an example. Long time ago this has been discussed in the old forum, and the majority (including me) voted for the current solution, essentially because of the same reason why the original HP 35 had a \(x^ y\) key and the general root function is \(\sqrt[x]y\).

Then consider the Σ+ and Σ- keys. I always found the way data pairs have to be entered somewhat awkward – why can't I simply type x [ENTER] y? Why do I have to key in y first and then x? It's because the X-register is used for x-related data and the Y-register for those that relate to y: think of means, standard deviations, sums, etc.

That's why I think that sometimes the POLA principle might not be the best possible solution. Looking at the regularized incomplete Beta function, the usual notation \(I_x (a, b)\) would suggest using X for x and the higher stack levels for the two parameters a and b. Which means that x is entered as the last value, not the first. This also makes sense if you think of the LastX command: Which of the three values would you like to recover via RCL L? I'd vote for x here. Would you agree?

That's why my suggestion for the argument sequence of this function is a [ENTER] b [ENTER] x or, maybe even better, b [ENTER] a [ENTER] x. Remember, it's RPN after all. ;-)

Code:

Z:  a          b
Y:  b    or    a
Z:  x          x

In other words, the 34s function would evaluate \(I_x (z, y)\) or \(I_x (y, z)\). Either way is fine by me, with some preference for the latter. Which happens to be the way that is stated in the manual, so it seems to be the obvious choice here. ;-)

Dieter

[Dieter]

(05-01-2014 09:55 PM)Paul Dale Wrote:  Changing the documentation will be the quickest fix here. Reordering the arguments to this function in the code is straightforward but several of the statistical distributions will also have to be updated to match.

While we're at it: I did some tests with the Student and Chi² quantile functions and found some cases that IMHO require too much time before they return a result. 15 seconds (or more than a minute in single cases) can be vastly improved. For instance with a Halley-solver (at least for these two distributions, the second derivative can be easily obtained) and a shorter (yet accurate) way of estimating a first guess for the Student quantile (even without slow logs ;-)). All in all, usually not more than three iterations or 2 to 5 seconds are required for 30+ digit accuracy. In user code, that is. Now imagine how fast this would run in compiled C code.

(05-01-2014 09:55 PM)Paul Dale Wrote:  I guess the big question is what order should the arguments be in? I don't see any benefit to reordering them any differently from what they are currently but am willing to be convinced otherwise.

Please read my other post in reply to Thomas Klemm in this thread. The current order as documented in the manual is fine. It just needs to get implemented this way. ;-)

Dieter

[walter b]

(05-02-2014 08:48 PM)Dieter Wrote:  In other words, the 34s function would evaluate \(I_x (z, y)\) or \(I_x (y, z)\). Either way is fine by me, with some preference for the latter. Which happens to be the way that is stated in the manual, so it seems to be the obvious choice here. ;-)

For some strange reason, I can perfectly live with that. Now it's up to Pauli ...

d:-)

[Thomas Klemm]

(05-02-2014 08:48 PM)Dieter Wrote:  the usual notation \(I_x (a, b)\) would suggest using X for x and the higher stack levels for the two parameters a and b

Mathworld uses:



Mathematica uses:



What exactly do you call "usual notation"?

Cheers
Thomas

[Dieter]

(05-02-2014 10:22 PM)Thomas Klemm Wrote:  Mathematica uses:

...something that has to be ASCII-compliant. ;-)

And now take another look at the Wolfram website you quoted and see what's written one line below at "Traditional notation". ;-)

(05-02-2014 10:22 PM)Thomas Klemm Wrote:  What exactly do you call "usual notation"?

Abramowitz & Stegun, p. 944 ff. Note that his was published by the National Bureau of Standards. How much more "standard" can it get ?-)

Just for comparison: does anyone have a Bronstein–Semendjajew at hand?

Dieter

[walter b]

(05-02-2014 10:34 PM)Dieter Wrote:  does anyone have a Bronstein–Semendjajew at hand?

For sure (13th edition of 1973 ). Alas, no Beta in the index.

d:-I

[Thomas Klemm]

(05-02-2014 10:34 PM)Dieter Wrote:  Abramowitz & Stegun, p. 944 ff. Note that his was published by the National Bureau of Standards. How much more "standard" can it get ?-)

Just for comparison: does anyone have a Bronstein–Semendjajew at hand?

So we have at least two, maybe three "standards". Maybe not a good idea to base an argument on the choice of a free variable.

[Manolo Sobrino]

This is actually funny. You usually call the first complex variable z, and the first real variable x. As the definitions of the gamma/beta functions start with real functions and then get analytically continued to the complex plane, it is implied that a definition with x is for real arguments and if we considered complex arguments, it would involve z instead.

The Abramowitz notation is for real values of x, because only real values are considered there. The current DLMF from NIST uses x because it is assumed that a>0, b>0, 0<=x<=1, and then "it is straightforward to continue most results analytically to other real values of a, b, and x, and also to complex values."

The Wolfram notation is just this, with z, as the Mathematica implementation uses complex values.

But alas, the stack has levels traditionally labelled x, y, z, t. Well, they're just levels 1, 2, 3, 4 with input sequence {4, 3, 2, 1}. Label names don't have to mean anything. Whatever complex or real variable is considered: x or z, the variable name by itself shouldn't determine the position in the stack.


If you enter things as they appear in standard notation you have the y^x ordering, that is var1^var2 has the input sequence {var1, var2} and a stack [x = var2, y = var1, z = 0, t = 0].

This way I("x or z"; a, b) would be I(var1; var2, var3) with input sequence {var1, var2, var3} and a stack [x = var3, y = var2, z = var1, t = 0].

Now, if we want to give the first position in the stack to the "main" variable, that is, if we consider the other variables as parameters, then I("x or z"; a, b) would have two choices for the stack [x = "x or z", y = a, z = b, t = 0] or [x = "x or z", y = b, z = a, t = 0]. There's no way to choose one of them as more logical than the other unless we adopt another convention. And that would likely be that the values in the stack must have the same order in which they appear in the standard notation left to right written expression, i.e. the first one. The input sequence would be: {b, a, "x or z"}.

According to this, if we consider the function with value var1^var2 and we think that the "main" variable is var1, our printed y^x is wrong and should be x^y. If the "main" variable is var2 it's all fine.

As this is a source of potentially unending confusion, at least for variable lists I'd just enter things as they appear in written expressions, independently of the names of the variables.

[Paul Dale]

(05-02-2014 09:03 PM)Dieter Wrote:  While we're at it: I did some tests with the Student and Chi² quantile functions and found some cases that IMHO require too much time before they return a result. 15 seconds (or more than a minute in single cases) can be vastly improved. For instance with a Halley-solver (at least for these two distributions, the second derivative can be easily obtained) and a shorter (yet accurate) way of estimating a first guess for the Student quantile (even without slow logs ;-)). All in all, usually not more than three iterations or 2 to 5 seconds are required for 30+ digit accuracy. In user code, that is. Now imagine how fast this would run in compiled C code.

I can only agree with you here. I'm using a fairly generic Newton based solver for all the quantile functions except the normal QF. I know it can be very slow, that is the price for having a single general solver.

The 34S certainly doesn't have space for a C implementation of these. I originally implemented all the distribution functions in C but had to rewrite them as keystroke programs for space reasons. Adding additional custom solvers is likewise going to consume precious bytes. Thus, the 34S is unlikely to see any distribution speed ups. However, the 31S quite possibly will. I'd like to rewrite the distribution code there in native C -- I can't just drop in the original distribution code, we moved on algorithmically.

So, yes I'd love to hear your thoughts on improving these functions.


- Pauli