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.

[Han]

Are the settings of your calculator and the emulator set to be the same?

[Ben Fairbank]

I believe so -- I went into the settings via shift-CAS on the Prime and on the emulator and they all checked as the same. Are there others I should have a look at?
Thank you for the response.
Ben

[Werner]

Wolframalpha confirms the emulator's result:
(integrate (x^3*(1-x)^5) dx from x=0 to y)*beta(4,6)-0.25=0
y = -1.62561, and it's the only root.

UPDATE: and that shows the error you (and I) made...
If I change it to *divide* by beta(4,6):

(integrate (x^3*(1-x)^5) dx from x=0 to y)/beta(4,6)-0.25=0
There are three roots,
y = -0.17836
y = 0.290986
y = 1.3936

[parisse]

Your problem is not well defined, because an integral is defined up to a constant, you should use a integral with boundaries. Should be e.g.
fsolve(int(x^3*(1-x)^5/Beta(4,6),x,0,xx)-.25,xx,.5)
Anyway, I'm unable to reproduce your problem with the latest source build. I get 0.29...
If you want to restrict to an interval, you can enter e.g.
fsolve(int(xx^(p-1)*(1-xx)^(q-1)/Beta(p,q),xx)-.25,xx=-1..1)
then dichotomy is used (and here it will find more than one root), otherwise an iterative Newton-like method is used.

[Ben Fairbank]

Thank you for your comments, but I do not understand your reference to the "error." My post shows division, not multiplication, by the Beta function, no?
Since the only x-axis region of interest is that from 0 to 1, the only root of interest is the one that I found, .29. It is still not clear to me how to restrict the computational area of interest to the 0 to 1 interval.

Ben
Ben

(03-13-2014 07:15 AM)Werner Wrote:  Wolframalpha confirms the emulator's result:
(integrate (x^3*(1-x)^5) dx from x=0 to y)*beta(4,6)-0.25=0
y = -1.62561, and it's the only root.

UPDATE: and that shows the error you (and I) made...
If I change it to *divide* by beta(4,6):

(integrate (x^3*(1-x)^5) dx from x=0 to y)/beta(4,6)-0.25=0
There are three roots,
y = -0.17836
y = 0.290986
y = 1.3936

[DrD]

No expert here, but it looked like Parisse answered your question:

If you want to restrict to an interval, you can enter e.g.
fsolve(int(xx^(p-1)*(1-xx)^(q-1)/Beta(p,q),xx)-.25,xx=-1..1)

Does:

fsolve(int(xx^(p-1)*(1-xx)^(q-1)/Beta(p,q),xx)-.25,xx=-0..1)

accomplish your goal?

[Werner]

(03-13-2014 01:09 PM)Ben Fairbank Wrote:  Thank you for your comments, but I do not understand your reference to the "error."
Ben

If you multiply by beta(4,6) instead of dividing by it, the root becomes -1.63 and it is the only root. If you divide, -1.63 is NOT a root.
So I assumed you mistyped it, as I did.
Apologies if you didn't - but you have to admit that it is quite a coincidence?

Werner

[Ben Fairbank]

Parisse --

Thank you for the reply. I do not know how I overlooked the method of restricting an integral to an interval -- now I know. But that was not the real problem.

Consider these two lines:

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

They are identical except that one uses the number 5 and the other uses q-1 (q elsewhere given the value of 6, and p = 4). The lines should give the same result if p = 4 and q = 6. Yet the first line succeeds (results in a value of .29) but the second line fails on a Prime with software version 5447, CAS version 1.1.0-27. This peculiar anomaly is stopping a program I have written from running...

Ben Fairbank

(03-13-2014 07:47 AM)parisse Wrote:  Your problem is not well defined, because an integral is defined up to a constant, you should use a integral with boundaries. Should be e.g.
fsolve(int(x^3*(1-x)^5/Beta(4,6),x,0,xx)-.25,xx,.5)
Anyway, I'm unable to reproduce your problem with the latest source build. I get 0.29...
If you want to restrict to an interval, you can enter e.g.
fsolve(int(xx^(p-1)*(1-xx)^(q-1)/Beta(p,q),xx)-.25,xx=-1..1)
then dichotomy is used (and here it will find more than one root), otherwise an iterative Newton-like method is used.

[parisse]

Yes, I checked your input with the latest source, and it works. Therefore this problem has been fixed (just wait for the next update...)

[toml_12953]

(03-14-2014 08:22 AM)parisse Wrote:  Yes, I checked your input with the latest source, and it works. Therefore this problem has been fixed (just wait for the next update...)

When? WHEN? WHEN???

[Tim Wessman]

[NOTE - this is not directed specifically at anyone, but your post reminded me I need to do this more often]

To clarify in hopes of making sure people aren't getting wrong idea:

HP has an ongoing working relationship with Bernard. Bernard fixes issues in his main source tree and then HP can choose when/if to integrate those with our sources. Bernard also does lots of excellent work with Geogebra, and I assume with reports he gets from his xcas/giac users.

However, Bernard does not speak for HP. Saying that something is "fixed" means that he believes the problem was resolved in his area or source. That does not mean the fix will ever be released by HP in a new firmware, or that HP is under any sort of obligation to do so.

While I know people are excited about things, please do not take comments from Bernard or anyone else (including Cyrille or myself) as official information or promises from HP.

I'm glad people are excited! However, note that asking about when/where is not anything that anyone who really knows will tell you.

[parisse]

Yes, "fixed" means fixed in giac source. "Wait for the next update" means about 1 week for Xcas and nothing about HP firmware.