Services & ResourcesWolfram ForumsMathGroup Archive

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

Intermediate Evaluation in FindMinimum

• To: mathgroup at smc.vnet.net
• Subject: [mg86617] Intermediate Evaluation in FindMinimum
• From: "Ausman, Kevin" <ausman at okstate.edu>
• Date: Sat, 15 Mar 2008 03:10:01 -0500 (EST)
• References: <200803130933.EAA01644@smc.vnet.net>

As I promised in the LevenbergMarquardt thread, I am now posting a full 
working example of the problem I am having. Fortunately, as I was 
writing this up I discovered some errors that I have fixed. 
Unfortunately, there is still something wrong (very different from what 
I was seeing earlier). Fortunately, I was able to distill the problem 
down to a much simpler one. So, here we go.

First, let's define our model:

testmodel[k1_?NumericQ, a0_?NumericQ, totalTime_?NumericQ, t_] :=
  testmodel[k1, a0, totalTime, t] =
   Piecewise[{{0, t < 0},
        {a[t] /. First[NDSolve[{a'[t] == -k1 a[t], a[0] == a0}, {a},
                {t, 0, totalTime}]], t >= 0}}];

Now let's convolve the model with an experimentally-determined instrument function:

testsimulation[xptXAxis_, instfn_, k1_?NumericQ, a0_?NumericQ] :=
  testsimulation[xptXAxis, instfn, k1, a0] = (
    instFnY = Transpose[instfn][[2]];
    normalizedinstfn = instFnY/Total[instFnY];
    xAxisStep = xptXAxis[[2]] - xptXAxis[[1]];
    modelXAxis = Join[Table[t, {t, xptXAxis[[1]] -
         xAxisStep (Length[normalizedinstfn] - 1),
        xptXAxis[[1]] - xAxisStep, xAxisStep}], xptXAxis];
    rawModelYAxis = Table[testmodel[k1, a0, Last[xptXAxis], t]
        /. (t -> modelXAxis[[i]]), {i, 1, Length[modelXAxis]}];
    ListConvolve[Transpose[instfn][[2]], rawModelYAxis]);

I need a function to evaluate the square error:

chiSquare[xpt_, mdl_] := chiSquare[xpt, mdl] = Total[(xpt - mdl)^2];

And I need both a sample instrument function and some sample data to fit to:

tspacing = 10^-6;
testInstFunc = Table[{i, 0}, {i, -3 tspacing, -1 tspacing, tspacing}]~Join~
        Table[{i, Exp[-i/tspacing]}, {i, 0, 20 tspacing, tspacing}];
xptX = Table[i, {i, -50 tspacing, 500 tspacing, tspacing}];
mdlY = testsimulation[xptX, testInstFunc, 7000, 0.003];
xptY = (mdlY + Table[RandomReal[NormalDistribution[0, 0.00002]], {i, 1,
      Length[mdlY]}]);

Show[ListPlot[Transpose[{xptX,xptY}]],ListLinePlot[Transpose[{xptX,startingpoint=testsimulation[xptX,testInstFunc,6500,0.0025]}],PlotStyle=AERed=
]]

chiSquare[xptY,startingpoint]

Performing a fit, however, ends up taking us further away from the minimum:

testSolution=FindMinimum[csq=chiSquare[xptY,testsimulation[xptX,testInstFunc,k,a0]],{{k,6500},{a0,0.0025}},StepMonitor:>Print[ListPlot[csq]]]

And as you can see from the StepMonitor output, chiSquare is returning a list instead of a scalar. I can reproduce this problem in an even simpler format like this:

chiSquare[xptY,testsimulation[xptX,testInstFunc,k,a0]]/.Last[testSolution]

This returns a long list. But move the replacement rule slightly and it returns a scalar:

chiSquare[xptY,testsimulation[xptX,testInstFunc,k,a0]/.Last[testSolution] ]

What is up with that?

Any help at all is greatly appreciated!

Kevin

• References:
  • LevenbergMarquardt
    • From: "Ausman, Kevin" <ausman@okstate.edu>