       Re: nonlinear programming with differential-algebraic constraints (2)

• To: mathgroup at smc.vnet.net
• Subject: [mg52129] Re: [mg52084] nonlinear programming with differential-algebraic constraints (2)
• From: DrBob <drbob at bigfoot.com>
• Date: Thu, 11 Nov 2004 04:53:09 -0500 (EST)
• References: <200411100945.EAA11228@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

```>> The main problem seems to be that the differential
>> equation does not get evaluated in each step but only once.

model (as written) has the same value every time, so why do you want it evaluated more than once? I think you wanted it to depend on "a", but you've written it with a->0.05 hardwired in.

I think this does what you want:

f[t_] := Exp[-0.1*t] + Random[Real, {0, 0.1}];
data = {#, f@#} & /@ Table[t, {t, 0, 10, 1}];
{xData, yData} = Transpose@data;

Clear[model, square]
model[(a_)?NumericQ] :=
y /. First[NDSolve[
{Derivative[y][t] == (-a)*y[t],
y == 1}, y, {t, 0, 10}]]
square[a_] := (#1 . #1 & )[
yData - model[a] /@ xData]
NMinimize[{square[a], a >= 0}, a]

{0.00638361, {a -> 0.0818508}}

Here's a more direct (but fully equivalent) approach:

Clear[f, model, square]
model[a_][t_] := Exp[(-a)*t]
f[t_] := model[0.1][t] + Random[Real, {0, 0.1}]
data = Table[{t, f[t]}, {t, 0, 10, 1}];
{xData, yData} = Transpose[data];
square[a_] := (#1 . #1 & )[
yData - model[a] /@ xData]
NMinimize[{square[a], a >= 0}, a]

This always underestimates a, however, because the noise you've added is always positive, and the error function doesn't take that into account. You can do better this way:

Clear[f, model, square]
model[a_][t_] := Exp[-a*t]
f[t_] := model[0.1]@t + 0.1Random[]
{xData, yData} = Transpose@Table[{t, f@t}, {t, 0, 10, 1}];
square[a_] := #.# &[yData - a/2 - model[a] /@ xData]
NMinimize[{square@a, a >= 0}, a]

{0.0139239, {a -> 0.0999625}}

and this does pretty well with no minimization at all:

Clear[f, model, inverse]
model[a_, t_] = Exp[-a*t];
inverse[a_]@{x_, y_} := -Log[y - .5a]/x
f[t_] := model[0.1, t] + 0.1Random[]
data = Table[{t, f@t}, {t, 1, 10}];
Mean[inverse[0.1] /@ data]

0.100146

I didn't put a lot of "science" into that. I'm just subtracting the average bias and solving for a in y - .5 a == Exp[-a x]. The LHS of that equation is an unbiased estimate of y (given a and x), but the result isn't necessarily an unbiased estimate of a (given y and x).

Bobby

On Wed, 10 Nov 2004 04:45:25 -0500 (EST), Joerg Schaber <schaber at molgen.mpg.de> wrote:

> I made a little example to illustrate the problem:
>
> I want to fit a model which is stated as a system of differenttial
> equations to some data.
>
> (* curve to be fitted *)
> f[t_]:=Exp[-0.1*t]+Random[Real,{0,0.1}];
> data={#,Evaluate[f[#]]}&/@Table[t,{t,0,10,1}];
>
> ListPlot[data,PlotStyle\[Rule]{RGBColor[0,0.2,0.8],PointSize[0.02]}];
>
> (* the parameter 0.1 is to be found *)
> (* the model as a differential equations *)
> model:=NDSolve[{y'[t]==-a*y[t],y==1}/.a->0.05,y,{t,0,10}];
>
> (* minimizing sum of squares *)
> NMinimize[{squares=Plus@@Table[(Table[{#,y[t]/.Evaluate[model][]/.t->#}&/@
> 		Table[t,{t, 0, 10, 1}]][[i, 2]]-data[[i,2]])^2,
> 		{i, Length[data]}], a > 0}, a,
> 		Method -> {"DifferentialEvolution"},
> 		StepMonitor :> Print[squares, a],
> 		WorkingPrecision -> 8]
>
>
> The main problem seems to be that the differential equation does not get
> evaluated in each step but only once.
>
> So if anybody could make this example work that would be great.
>
> best,
>
> joerg
>
>
>
>

--
DrBob at bigfoot.com
www.eclecticdreams.net

```

• Prev by Date: Re: Re: newbie question DSolve (revisited)
• Next by Date: Re: A NIntegrate question
• Previous by thread: nonlinear programming with differential-algebraic constraints (2)
• Next by thread: Re: nonlinear programming with differential-algebraic constraints (2)