Message Boards Message Boards

Best solver question?

Is there documentation focusing on differences among various solvers in Mathematica? It's difficult to pick out key distinctions by reading documentation for functions separately.

More immediately, which of the solvers are best for finding bounds of integration, where performance matters. Depending on df, equations like the one below can take an inconveniently long time to solve.

Function??[Integrate[PDF[StudentTDistribution[df], x], {x, -t, t}] == 95/100, t, Reals]
POSTED BY: Jay Gourley
4 Replies

Depends on what exactly you want. Is this where you plan to give df as a numeric value and you want to solve numerically for t? For that maybe like so.

ifunc[df_?NumericQ, t_?NumberQ /; t > 0] := 
 NIntegrate[PDF[StudentTDistribution[df], x], {x, -t, t}]

FindRoot[ifunc[3, t] == .95, {t, 1}]

(* Out[216]= {t -> 3.18245} *)

I do not see much hope for getting an exact or symbolic solver for this.

POSTED BY: Daniel Lichtblau

One can go to higher precision and also (in this case) avoid NIntegrate like so.

stdInt[df_, t_] = 
 Integrate[PDF[StudentTDistribution[df], x], {x, -t, t}, Assumptions -> t > 0 && df > 0]

(* Out[218]= (2 t Hypergeometric2F1[1/2, (1 + df)/2, 3/2, -(t^2/df)])/(Sqrt[df] Beta[df/2, 1/2]) *)

FindRoot[stdInt[300, t] == 95/100, {t, 1}, WorkingPrecision -> 20]

(* Out[227]= {t -> 1.9679030112610870301} *)

This seems to work just fine. However, if ever you need a better starting point for the root search, your approximation is a good way to get one.

normInt[t_] = 
 Integrate[PDF[NormalDistribution[0, 300/298], x], {x, -t, t}, Assumptions -> t > 0]

(* Out[228]= Erf[(149 t)/(150 Sqrt[2])] *)

FindRoot[normInt[t] == 95/100, {t, 1}, WorkingPrecision -> 20]

(* Out[230]= {t -> 1.9731181052416653378} *)

Unfortunately the higher precision is still necessary because FindRoot will resort to the precision of its inputs if that's lower than the requested WorkingPrecision. (I myself find this really annoying when it's due to a starting point, But we've yet to change that.)

Last is a method that works for this particular function. We take the log of the exact integral after function-expanding, and only evaluate for t set to an explicit real value. This avoids some conversion of the hypergeometric to numerically unusable expressions (when using machine precision).

logstdInt[df_, t_] = FunctionExpand[Log[stdInt[df, t]], Assumptions -> {t > 0, df > 0}] /. 
  Log[aa_.*Gamma[xx__]^bb_.] :> Log[aa] + bb*LogGamma[xx]

(* Out[245]= 
Log[2] - 1/2 Log[df \[Pi]] + Log[t] + 
 Log[Hypergeometric2F1[1/2, (1 + df)/2, 3/2, -(t^2/df)]] - LogGamma[df/2] + LogGamma[(1 + df)/2] *)

logstdIntNumeric[df_?NumberQ, t_Real] := logstdInt[df, t]

FindRoot[logstdIntNumeric[300, t] == Log[.95], {t, 1}]

(* Out[249]= {t -> 1.9679} *)
POSTED BY: Daniel Lichtblau

Thanks, Daniel Lichtblau. I don't need an exact solution. My original post was unclear on that point. Also, I should have been more precise about performance. I'm not worried about a few seconds.

Your FindRoot idea works great for small degrees of freedom. It's blazing fast and precise enough for my purposes. But it falls apart for larger df, say 300 or above. It hits singularities when I try.

In my clumsy attempts, I've gotten best results with NSolve and Integrate. But in the same cases where FindRoot fails, NSolve becomes very slow, I think because it relies on Integrate to find an exact solution. By slow, I don't mean a few seconds. I mean many minutes. In most cases though, it ultimately returns a credible result. To avoid long waits, I substitute a Normal distribution with a df/(df-2) adjustment. But I can't help thinking there's a better way.

Another thing I don't understand about FindRoot with Nintegrate, is NIntegrate throws off errors "x = t is not a valid limit of integration" even in cases when FindRoot ultimately returns a credible result.

POSTED BY: Jay Gourley

Great contribution, Daniel Lichtblau. Thanks. No telling how long it would have been before I stumbled onto either of those approaches.

POSTED BY: Jay Gourley
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract