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

```

• Prev by Date: Re: Increasing RandomInteger speed
• Next by Date: RE: Re: LevenbergMarquardt