MathGroup Archive 2008

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

Search the Archive

Intermediate Evaluation in FindMinimum

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,



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


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:


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!


  • Prev by Date: Re: Increasing RandomInteger speed
  • Next by Date: RE: Re: LevenbergMarquardt
  • Previous by thread: LevenbergMarquardt
  • Next by thread: RE: Intermediate Evaluation in FindMinimum