Message Boards Message Boards

Get the global minimum fitting?

Posted 7 years ago

Hello, guys.

I have a problem in fitting. Whenever I try to fit, I think it goes to local minimum. Also I saw the error messages which I couldn't fix properly. The message is: " The step size in the search has become less than the tolerance prescribed by the PrecisionGoal option, but the gradient is larger than the tolerance specified by the AccuracyGoal option. There is a possibility that the method has stalled at a point that is not a local minimum" Can you help me to get a global minimum ? Thank you very much.

The code that I used is following. And about the data file, I attached it.

dat = Import["d:\\Data\\P2d5Hs.txt", "Table"];
graph1 = ListPlot[dat, PlotRange -> All]
response[t_, x_, ?_, t0_, k1_, k2_, k3_, ?1_, ?2_, ?3_] := 
  Exp[-(x - t + t0)^2/(?^2/2)] ( 
    UnitStep[
      x] {k1 (1 - Exp[-x/?1]) (Exp[-x/?2] + 
          Exp[-x/?3]) + k2 (1 - Exp[-x/?2]) + 
       k3 (1 - Exp[-x/?3])});
convol[t_, ?_, t0_, k1_, k2_, 
   k3_, ?1_, ?2_, ?3_] = 
  Integrate[
   response[t, x, ?, t0, k1, k2, 
    k3, ?1, ?2, ?3], {x, -2, 13}, 
   Assumptions -> {{t, ?, t0, k1, k2, 
       k3, ?1, ?2, ?3} ? 
      Reals, ? > 0 && ?1 > 0 && ?2 > 0 && ?3 > 
       0 && k2 > 0 && k3 > 0}];
sol = NonlinearModelFit[dat, 
   convol[t, ?, t0, k1, k2, 
    k3, ?1, ?2, ?3], {{t0, 0.05}, {?, 
     0.31}, {k1, 0.1}, {k2, 0.02}, {k3, 0.005}, {?1, 
     0.05}, {?2, 0.4}, {?3, 0.9}}, t, 
   ConfidenceLevel -> .998, MaxIterations -> 10000];
fitgraph = Table[{t, sol[t]}, {t, dat[[1, 1]], dat[[-1, 1]], 0.05}];

Show[graph1, Plot[sol[t], {t, -1.2, 13}, PlotRange -> All]]
sol["ParameterTable"]
sol["RSquared"]
Attachments:
POSTED BY: hwoarang Polar
6 Replies

Moderator -- Please add this to the SystemModeler thread -- people there may have some insight on modeling the physics at play here.

POSTED BY: Neil Singer

The section of the book that I posted was (a very humorous) discussion about your approach. I hope you were not offended -- I thought it was worth reading -- in fact the whole book is a "must read" -- insightful and humorous. The book was written in 1970 so experimentalists for well over 50 years have been falling into your same numerical trap! While it is tempting to take data, then fit some dynamic responses to the data and use that to investigate the physics behind the problem -- this approach is generally intractable. In certain cases it can work (for example if your system response is underdamped the problem is significantly easier -- i.e. "modal analysis").

How to approach this problem:

  1. The physical problem: Why do assume three damped time constants? Is there some physics reason that you believe that there are two or three overdamped poles at play here? If not, and you are just assuming the dynamics and trying to see what fits then you are squarely trying to solve Acton's unsolvable isotope problem.

  2. Model the physics. I would start by going to first principles and modeling the dynamics of the problem. You will generally find that material properties will pin down the dynamics and it will suddenly become easier to fit. In the isotope example this would be equivalent to assuming certain chemical elements in the sample based on knowledge about where the sample came from (your "model"). You can even use Wolfram's SystemModeler to model the dynamics if you understand the underlying physics. (SystemModeler is an excellent program for this!)

  3. Do additional experiments to isolate particular physical phenomena. In the isotope example this would be equivalent to doing a series of chemical tests to determine the elements in the sample. By adding more information you can inform your model and the data fit will become easier/possible.

I hope this helps -- unfortunately I am not familiar enough with laser dynamics to be more specific. Maybe others in the community will have additional suggestions...

POSTED BY: Neil Singer
Posted 7 years ago

Actually, the experimental graph is the transient reflectivity data of the metal induced by the laser pump. The time scale (x-axis of the data) is picosecond. The ideal goal is to compare fitting results between one-free parameter and two-free parameter equations. Then I want to endow physical meanings on parameters for each equation. So, to get parameter values as accurate as possible is quite important. Now I got the fitting value for one-free parameter equation. But I am stuck in getting fitting values for two-free parameter equation. Maybe, as you suggested, I need to think of other ways like modifying the equation a bit or .... If you or other people have some suggestions or comments, always welcomed. Thank you again,

Best regards,

POSTED BY: hwoarang Polar

hwoarang,

Yes, the problem gets worse with noise. However, with three time constants there will be multiple solutions that are all nearly as good. If you have information from other sources that helps you identify the time constants, then you can do your curve fit easily. You can even be missing one time constant and probably solve it.

What are you physically modeling and what is your ideal goal? -- it might help someone here to make a recommendation on an approach to navigate around this numerical problem.

Regards

POSTED BY: Neil Singer
Posted 7 years ago

Thank you very much for valuable information. At the same time, I am also aware of this possibility even though I was not sure 100 %. And one question that I have now is... If we want to fit the data curve which is exactly same to the equation with which we are trying to fit, we should get the perfect fitting without any problem (meaning that free parameters are well determined by fitting). Am I correct ? If this is right, I think that the noise in the experimental data increases the number of local minima around the global minimum. The extent of this aspect becomes severe when the tolerance range of the function for the fit is narrow (like two exponential functions in the book above). Then we call it 'ill-conditioned' (kind of rate ~ noise level of the data/tolerance range of the equation) Do I say correctly ? I am not sure if these messages will be transferred well to you.

Thank you again.

Have a good day!

POSTED BY: hwoarang Polar

hwoarang,

The problem you are seeing is not a Mathematica problem. It appears that you are trying to fit the time constants of exponentials to real data. This is recognized as an ill conditioned problem. You need to add more information in order to "solve" it. For an excellent discussion of this problem look at Forman Acton's Book "Numerical Methods that (usually) work". The section entitled "Interlude: What Not to Compute" (Sorry!) covers your situation: I reproduced the two relevant pages below.

enter image description here

enter image description here

POSTED BY: Neil Singer
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