MathGroup Archive 2005

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: nonLinearFit and nonLinearRegress

John Hawkin wrote:
> I'm trying to do a nonLinear fit on some data, fitting a parameter
> which appears in the limits of the integral that I'm trying to fit. 
> When I do a nonLinearFit, I get about a page of errors, but it does
> work.  It takes about 10 minutes to do it for 20 points with a
> reasonable initial estimate.  However when I use the nonLinearRegress
> command, even with the best fit parameter from nonLinearFit given as
> the initial guess, it takes hours to work.  Is it possible there is a
> way to speed this up, and is this normal?
> The commands that I'm using to integrate this function are:
> h[B]:=2/(B pi) NIntegrate[x Sin[x] Exp[-(x/B)^3/2], {x, 0, inf},
> Method-Oscillatory];
> holts[c_, d_, Q_]:= NIntegrate[h[B], {B, c/Q, d/Q}];
> where holts is the function I am fitting (the h[B] integral is known
> as the Holtsmark integral).  The fitting command I'm using is:
> NonlinearFit[fitData, holts[c,d,Q], {c,d}, {Q, estimateQ}];
> where fitData contains sets of size 3, with the values of c and d as
> the first two values, and the data point that I'm fitting the integral
> to as the 3rd value.  If anyone has any ideas on how I can perform the
> nonlinear regression fast enough so that I can do it many times (maybe
> 50), I would greatly appreciate it.  Thanks very much,
> John Hawkin

One possibility for speeding it would be to use the analytic form for h.

n[9]:= InputForm[i1 = Integrate[x*Sin[x]* Exp[-(x/B)^3/2], {x, 0, 
Infinity}, Assumptions->B>0]]

     ((-1)^(1/3)*BesselI[-2/3, ((2/3 + (2*I)/3)*B^(3/2))/Sqrt[3]] -
      BesselI[2/3, ((2/3 + (2*I)/3)*B^(3/2))/Sqrt[3]] -
      (-1)^(1/3)*BesselJ[-2/3, ((2/3 + (2*I)/3)*B^(3/2))/Sqrt[3]] +
      BesselJ[2/3, ((2/3 + (2*I)/3)*B^(3/2))/Sqrt[3]]) +
    9*HypergeometricPFQ[{1}, {1/3, 2/3, 5/6, 7/6}, -B^6/2916]))/27

This can be evaluated at machine precision for B not very large. If your 
  upper bound d/Q is large for plausible values of Q you might need to 
increase WorkingPrecision in the outer integration holts[...].

I don't know if this approach will help but it might be worth trying. 
Obviously it relies on the fact that we can find an analytic form for 
h[]. Taking it a step further, here is what I get for

Integrate[i1, {B, c/Q, d/Q}, Assumptions->{Q>0,0<c<d}]

(((-2 + 2*I)*(-1)^(2/3)*Pi*
     (c^(5/2)*BesselI[-5/3, ((2/3 + (2*I)/3)*(c/Q)^(3/2))/Sqrt[3]] -
      d^(5/2)*BesselI[-5/3, ((2/3 + 
(2*I)/3)*(d/Q)^(3/2))/Sqrt[3]]))/Q^(5/2) +
        ((-2/3 + (2*I)/3)*(c/Q)^(3/2))/Sqrt[3]] -
      d^(5/2)*BesselI[5/3, ((2/3 - (2*I)/3)*(d/Q)^(3/2))/Sqrt[3]]))/
    ((1 + I)^(1/3)*Q^(5/2)) - (2*(-1)^(11/12)*(1 + I)^(2/3)*2^(1/6)*Pi*
     (c^(5/2)*BesselI[5/3, ((2/3 + (2*I)/3)*(c/Q)^(3/2))/Sqrt[3]] -
      d^(5/2)*BesselI[5/3, ((2/3 + 
(2*I)/3)*(d/Q)^(3/2))/Sqrt[3]]))/Q^(5/2) +
       (((-2*I)/27)*c^3)/Q^3] + Hypergeometric0F1Regularized[-2/3,
      (((-2*I)/27)*d^3)/Q^3]) +
   (3*(-(c^4*HypergeometricPFQ[{1}, {1/3, 5/6, 7/6, 5/3}, 
-c^6/(2916*Q^6)]) +
      d^4*HypergeometricPFQ[{1}, {1/3, 5/6, 7/6, 5/3}, 

As for verification of correctness, you're on your own. Assuming it is 
basically sound, I have no idea offhand whether this will suffer from 
branch cut problems when evaluated at specific {d,c,Q}.

Another possibility for approaching a problem such as this might be to 
use a series or similar such approximation for the first integrand, 
integrated termwise. This will not rely on vagaries of NIntegrate error 
control, but you'd need to control it yourself in choosing how many 
terms in the approximation to use.

Daniel Lichtblau
Wolfram Research

  • Prev by Date: Re: Add new option to Notebook[]
  • Next by Date: findroot
  • Previous by thread: Re: nonLinearFit and nonLinearRegress
  • Next by thread: Re: nonLinearFit and nonLinearRegress