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