[SciPy-User] expanding optimize.brentq

[josef.pktd at gmail.com]

On Mon, Mar 18, 2013 at 2:58 PM, Charles R Harris
<charlesr.harris at gmail.com> wrote:
>
>
> On Mon, Mar 18, 2013 at 11:13 AM, <josef.pktd at gmail.com> wrote:
>>
>> Just a thought, given the problems with fsolve
>>
>> In scipy.stats.distribution, we use an expanding brentq, which looks
>> like it works very well.
>> https://github.com/scipy/scipy/pull/216/files#L0R1175
>>
>> A question given that I have not much experience with root finders:
>>
>> Is there a general algorithm that has similar properties?
>>
>>
>> In these application we know that the function is strictly monotonic,
>> but we don't have a bound for brentq without checking. For some cases
>> I know one of the bounds (e.g. 0 for non-negative solutions).
>>
>> (I could use that right now as a replacement for fsolve.)
>>
>
> I assume you are talking about the 1-D case where you don't have bounding
> intervals to start with. There are various search methods that can be used
> to find such intervals, but they can't be guaranteed to work. They all come
> down to making trials of various points looking for sign reversals, and the
> searches are either subdividing a larger interval or creating larger and
> larger intervals by, say, doubling each time around. The addition of such
> methods would be useful, but they need a failure mode :(

Yes that's pretty much what I would like, and use in the rewritten
scipy.stats.distribution function.
But I was hoping there might be something more systematic than "home-made".

>
> Your application looks to be finding points on a monotone function on the
> interval [0, inf]. For that you could use the interval [0, 1] and map it to
> [0, inf] with the substitution x/1-x, although the point x == 1 will be a
> bit tricky unless your function handles inf gracefully. The accuracy will
> also tend to be relative rather than absolute, but that seems to be a given
> in any case. I suppose this approach falls under the change of variables
> method. Those tend to be problem specific, so I don't know if they would be
> a good fit for scipy, but certainly they could work well in specialized
> domains. It would also seem to be a good idea to match the asymptotic
> behavior if possible, which will tend to linearize the problem. So other
> options would be functions like log, erf, arctan, etc, for the substitution,
> but in those cases you probably already have the the inverse functions.

That might be a good idea to get brentq to influence the selection of
subdivisions, exponential instead of linear. One problem with brentq
is that I would have to give one huge bound, but the most likely case
is much smaller.
(example find sample size for power identity, we should get an answer
below 1000 but in some cases it could be 100.000)

One related problem for the bounded case in (0,1) is the open
interval, I use 1e-8, some packages use 1e-6 to stay away from the
undefined boundary points.
But again, there could be some extreme cases where the solution is
closer to 0 or 1, than the 1e-8.
(I don't remember if we found a solution to a similar problem in the
stats.distributions, or if I mix up a case there.)

Josef

>
> Chuck
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>