[Date Index]
[Thread Index]
[Author Index]
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>
*Reply-to*: drbob at bigfoot.com
*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[1][y][t] == (-a)*y[t],
y[0] == 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[0]==1}/.a->0.05,y,{t,0,10}];
>
> (* minimizing sum of squares *)
> NMinimize[{squares=Plus@@Table[(Table[{#,y[t]/.Evaluate[model][[1]]/.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)**
| |