RE: Fitting experimental data
- To: mathgroup at smc.vnet.net
- Subject: [mg85501] RE: Fitting experimental data
- From: "Ausman, Kevin" <ausman at okstate.edu>
- Date: Wed, 13 Feb 2008 04:06:29 -0500 (EST)
- References: <fopb53$bhu$1@smc.vnet.net> <47B056F2.9000401@umbc.edu>
Thanks for the suggestion. Unfortunately, I don't see how it helps my
particular problem.
I had hoped to avoid going into the full messy problem, but I probably
need to, so here we go.
GENERAL GOAL: an integrated interactive kinetic fitting routine.
GENERAL APPROACH: A module file that has a function which produces a
panel with a number of dynamic tabs:
1) Defining the kinetic model: specify the number of
time-dependent concentrations to follow (which determines both the
number of differential rate equations and the number of initial
concentrations to specify), the number of rate constants (which are
simply constant parameters in the differential rate equations), and the
rate equations themselves.
2) The data to fit: specify the experimental data to fit (which
is in the form of a list of {x,y} points) and the instrument function
(which is also a list of {x,y} points, and serves as a kernel for a
ListConvolve step. General idea is to create a numerically-integrated
simulation based on a set of parameters using the user-specified rate
equations, convolve it with the instrument function, and then determine
the goodness of fit by comparison with the experimental data).
3) The simulation parameters: specifies interactively the
initial values for all of the variable parameters (showing a fit to the
data at the same time), and determines which of them are variable and
which are to be held constant.
4) A close-up of the fit, the residuals, and a button that says,
"Start automatic fit".
WHAT CURRENTLY WORKS: The user interface works perfectly. All the values
and the graphs update interactively, as does the numerical integration
of the user-defined kinetic model followed by convolution with the
instrument function.
WHAT DOESN'T WORK: The automatic fit.
I have extracted the problem from the user-interface portion of the code
to try to cut down on how much code must be waded through, though what
remains is still pretty long.
Let's start with the model itself:
model[numRateConst_?NumericQ, rateConst_?(VectorQ[#, NumericQ] &),
numInitConc_?NumericQ, ic_?(VectorQ[#, NumericQ] &), rateEq_,
finalTime_?NumericQ, t_] := model[numRateConst, rateConst,
numInitConc,
ic, rateEq, finalTime,
t] =
Module[{systemOfEquations, speciesToSolveFor,
rateReplacementRules = {ToExpression["k1"] ->
rateConst[[1]],
ToExpression["k2"] -> rateConst[[2]],
ToExpression["k3"] -> rateConst[[3]],
ToExpression["k4"] -> rateConst[[4]],
ToExpression["k5"] -> rateConst[[5]],
ToExpression["k6"] -> rateConst[[6]],
ToExpression["k7"] -> rateConst[[7]],
ToExpression["k8"] -> rateConst[[8]],
ToExpression["k9"] -> rateConst[[9]],
ToExpression["k10"] -> rateConst[[10]]}},
systemOfEquations = Table[ReleaseHold[rateEq][[i]], {i, 1,
numInitConc}]~Join~
Table[ToExpression["a" <> ToString[i]
<> "[0]\[Equal]"
<> ToString[ic[[i]]]], {i, 1,
numInitConc}];
speciesToSolveFor = Table[ToExpression["a" <> ToString[i]], {i,
1, numInitConc}];
First[NDSolve[systemOfEquations /. rateReplacementRules,
speciesToSolveFor,
{ToExpression["t"], 0, finalTime}]]
]
The approach may seem strange, but it's the best way that I've found to
convert the use input of k1, k2, etc. for the rate constants into the
internal array representation of k[[1]], etc. The following test code
works for this function:
numRateConst = 1;
rateConst = {7000, 0, 0, 0, 0, 0, 0, 0, 0, 0};
numInitConc = 1;
ic = {0.0025, 0, 0, 0, 0, 0, 0, 0, 0, 0};
rateEq = {a1'[t] == -k1 a1[t], 0, 0, 0, 0, 0, 0, 0, 0, 0};
finalTime = 500 10^-6;
rawsimulation = model[numRateConst, rateConst, numInitConc, ic, =
rateEq,
finalTime, t]
Plot[(a1 /. rawsimulation)[t], {t, 0, finalTime}]
The resulting list of replacement rules has a single time-dependent
function, a1. The general case I am trying to solve, however, has the
experimental signal as a linear combination of multiple time-dependent
functions, a1, a2, etc. So, I use the following intermediate function:
signalFromModel[numRateConst_?NumericQ, rateConst_?(VectorQ[#, NumericQ]
&),
numInitConc_?NumericQ, ic_?(VectorQ[#, NumericQ]
&),
ec_?(VectorQ[#, NumericQ] &), rateEq_,
finalTime_?NumericQ,
t_] := signalFromModel[numRateConst, rateConst,
numInitConc,
ic, ec, rateEq, finalTime, t] =
(Block[{integratedModel},
integratedModel = model[numRateConst, rateConst,
numInitConc, ic,
rateEq, finalTime, t];
Total[Table[(ec[[i]]*ToExpression["a" <>
ToString[i]][t]),
{i, 1, numInitConc}]] /.
integratedModel
])
The apparent weirdness of the ToExpression bit was the way I found to
get around problems with enclosing this function in a module file
(otherwise the a1[t] in my function would refer to the module domain,
but the one entered in my user interface would refer to a global
version). If someone has a better way around this problem, I'd love to
hear it. Regardless, however, the code above works, as the following
test code demonstrates (reusing definitions from the first test code):
ec = {1, 0, 0, 0, 0, 0, 0, 0, 0, 0};
signalsimulation = signalFromModel[numRateConst, rateConst, =
numInitConc,
ic, ec,
rateEq, finalTime, t]
Plot[signalsimulation, {t, 0, finalTime}]
The final step in generating true simulated data is to generate the
simulation as a set of {x,y} pairs corresponding to the experimental x
axis and run ListConvolve with the experimental instrument function. My
following routine accomplishes this:
simulateData[xptXAxis_, instfn_, xoffset_, mdl_, t_, sln_] :=
simulateData[xptXAxis, instfn, xoffset, mdl, t, sln] =
(Block[{instFnX, instFnY, normalizedinstfn, xAxisStep,
modelXAxis, offsetModelXAxis, signalFunction,
rawModelYAxis},
instFnX = Transpose[instfn][[1]];
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];
offsetModelXAxis = modelXAxis - xoffset;
signalFunction = mdl /. sln;
rawModelYAxis = Table[Piecewise[{{0, offsetModelXAxis[[i]] < 0},
{signalFunction /. t -> offsetModelXAxis[[i]],
offsetModelXAxis[[i]] >= 0}}],
{i, 1, Length[offsetModelXAxis]}];
ListConvolve[Transpose[instfn][[2]], rawModelYAxis]
])
There are several steps here that may be confusing, but which I don't
think you need to worry about. (a) I normalize the instrument function.
(b) I prepend the x-axis with enough data points to ensure that the
resulting Convolved model matches the experimental data x-axis. (c) I
allow an offset to the x axis to account for signal delays. (d) My
experimental data starts before time-zero, and the signal has a sudden
onset at that time, so I define my raw Y axis simulation as a piecewise
function. Again, regardless of these details, the function works, as
this test code demonstrates:
tspacing = 10^-6;
testInstFunc = Table[{i, 0}, {i, -3 tspacing, -1 tspacing,
tspacing}]~Join~
Table[{i, Exp[-i/tspacing]}, {i, 0, 20 tspacing,
tspacing}];
ListPlot[testInstFunc, PlotRange -> Full]
xptX = Table[i, {i, -50 tspacing, 500 tspacing, tspacing}];
fullSimulation = Transpose[{xptX, simulationY = simulateData[xptX,
testInstFunc,
0, signalFromModel[numRateConst, rateConst, numInitConc, ic, ec,
rateEq,
finalTime, t], t, {}]}];
ListLinePlot[fullSimulation]
The last parameter in this function, sln, is there to allow the
following alternate notation, which also works:
rateConst2 = {kauto, 0, 0, 0, 0, 0, 0, 0, 0, 0};
ic2 = {icauto, 0, 0, 0, 0, 0, 0, 0, 0, 0};
fullSimulation2 = Transpose[{xptX, simulateData[xptX, testInstFunc, 0, =
signalFromModel[numRateConst, rateConst2, numInitConc, ic2, ec,
rateEq, finalTime, t], t, {kauto -> 7000, icauto -> 0.0025}]}];
ListLinePlot[fullSimulation2]
Now we have built-up enough of the problem to try (and fail) at the
automatic fit. First, let's generate some test experimental data (and
show the sum of the squared residuals as well just for kicks):
xptData = Transpose[{xptX, xptY = (simulationY +
Table[RandomReal[NormalDistribution[0, 0.00002]], {i, 1,
Length[simulationY]}])}];
Show[ListPlot[xptData], ListLinePlot[fullSimulation2, PlotStyle -> Red]]
Sum[(simulationY[[i]] - xptY[[i]])^2, {i, 1, Length[xptY]}]
An automatic fit, however, doesn't seem to work here:
result = FindMinimum[
(sim = simulateData[xptX, testInstFunc, 0,
signalFromModel[numRateConst, rateConst2, numInitConc,
ic2,
ec, rateEq, finalTime, t],
t, {kauto -> k1float, icauto -> ic1float}];
Sum[(xptY[[i]] - sim[[i]])^2, {i, 1, Length[xptY]}]),
{{k1float, 6000}, {ic1float, 0.003}},
Method -> "LevenbergMarquardt"]
I strongly suspect that the problem has to do with signalFromModel
requiring numerical values for the parameters to evaluate. For example,
this also does not work:
fullSimulation3 = Transpose[{xptX, simulateData[xptX, testInstFunc, 0, =
signalFromModel[numRateConst, rateConst2, numInitConc, ic2, ec,
rateEq, finalTime, t], t, {kauto -> ktmp, icauto -> ictmp}]}] /.
{ktmp -> 7000, ictmp -> 0.0025};
Unfortunately, I don't see how to get around this problem, and none of
my dozens of attempts have improved matters. Help?
Thanks!
Kevin Ausman
-------------------------
Prof. Kevin Ausman
Oklahoma State University
Department of Chemistry
342 Physical Science
Stillwater, OK 74078
Email: ausman at okstate.edu
Office: 405-744-4330
Fax: 405-744-6007
--------------------------