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>
- LevenbergMarquardt