Oddity Using Root Finder with Integrator (long)

[Ben Fairbank]

This is truly an odd one. I am writing a program using the Beta distribution to solve some statistical problems. As a step in the program, I have to answer such questions as “what is the ordinate value associated with the (for example) 25% point of the cumulative Beta distribution.” The Beta distribution is defined on the range 0 to 1, and is a probability density (or distribution) function and so always sums or integrates to unity when the limits of integration are 0 and 1. The distribution is defined by

y = xx^(pp-1)*(1-xx)^(qq-1)/Beta(pp,qq)

pp and qq in my application can be any integers greater than unity, and within the calculating range of the Prime. They are both almost always under 100, usually far under.

I need to find that value on the x (or, as I named it above, xx) axis where the definite integral from 0 to xx is equal to 0.25 (or .75, or whatever other value I need for the purposes of the program [note that those two values would allow the computation of the interquartile range of the function]).

I used the “FNROOT” root finder (User Guide, p. 315) with the “int” integration function (ibid. p. 326) to come up with the following expression (in this example I am letting pp = 4 and qq = 6):

FNROOT(int(xx^3*(1-xx)^5/Beta(4,6),xx)-.25,xx,0.50)

When using the CAS system on the Prime, that expression correctly evaluates to 0.29 (and does so surprisingly quickly.)

But then things get interesting. When I enter the same expression into the emulator’s CAS system, the returned value is -1.63, which fails on three counts, it is not between 0 and 1, it is negative, and it is wrong. Why? Do I then have to “tell” the emulator system to restrict xx to values from 0 to 1? I did not have to do so on the Prime, nor do I see how to do so. I checked to be sure all of the CAS settings matched on the Prime and the emulator, and they seemed to do.

Worse yet, when I use the variables p and q on the Prime, to which I have assigned values of 4 and 6, instead of the numerical constants 4 and 6 for the Beta function (I have to use variables, not constants, in my program), so the expression looks like this

FNROOT(int(xx^(p-1)*(1-xx)^(q-1)/Beta(p,q),xx)-.25,xx,0.50)

the root finder then fails with the error message “undef” even though the values are all the same, i.e., the variables pp and qq have the values 4 and 6. Why should it work for one but not the other?

But, weirder and weirder, changing ONLY the (q-1) exponent in the above expression to its numerical value, 5, “fixes” the problem and the correct result, .29, is again delivered even though all of the other p and q variables remain in the expression.

A totally different problem comes up when I use the line in my program; the error message “Bad Guess(es)” is displayed. I tested it only using numerical constants, not the variables pp and qq.

This all seems very odd, to say the least, and I hesitate to invoke the b-word to explain it, so I would appreciate any suggestions. As is probably needless to say, progress on the program has pretty much stopped until I clear this up.

With thanks and confidence that someone out there will come up with a perfectly obvious explanation.

Ben Fairbank March 12, 2014

P.S. And yes, I have upgraded to the most recent software revision a few days ago.