Parameterbecslés a Mathematica segítségével 
Kocsis Tibor GVAKWH
Vegyészmérnöki Kar, vegyészmérnöki szak, III. évfolyam

1. Paraméterekben lineáris, változókban lineáris eset

A feladat ez esetben a Q(p, 
          q) = (∑_(i = 1)^n)(((p (Subscript[x, i]) + q - (Subscript[y, i])))^2){(Subscript[x, i]),(Subscript[y, i]) adott} függvény azon (p,q) paramétereinek megkeresése, ahol a függvényérték minimális.
Ennek szükséges feltétele a deriváltak zérus volta, tehát a megoldandó egyenletrendszer(1.1):

(∂_1) Q(p, q)=2 (∑_(i = 1)^n)((Subscript[px, i]) + q - 
                            (Subscript[y, i])) (Subscript[x, i]) = 0
(∂_2) Q(p, 
                    q) = 2 (∑_(i = 1)^n)((Subscript[px, i]) + q - 
                          (Subscript[y, i])) = 0

Kifejezve q-t:
\(TraditionalForm\`q = \(∑\((y\_i\ x\_i)\) - p\ \(∑x\_i\^2\)\)\/\(n \(∑x\_i\)\)\)
melyet 1.1 második egyenletébe helyettesítve p-re kapjuk:
\(TraditionalForm\`p = \(\(∑\(x\_i\) \(∑y\_i\)\  - n \(∑\(x\_i\) y\_i\)\/n\)\/\(∑\(x\_i\) \(∑x\_i\) - n \(∑x\_i\^2\)\/n\) = \(\(x\&_\) y\&_ - \(x\ y\)\&_\)\/\(x\&_\^2 - \(x\^2\)\&_\)\)\)
ebből és a q-ra von. fenti egyenletből következően pedig:
\(TraditionalForm\`q = \(\(∑\((y\_i\ x\_i)\) - \ \(\(∑\(x\_i\) \(∑y\_i\)\  - n \(∑\(x\_i\) y\_i\)\/n\)\/\(∑\(x\_i\) \(∑x\_i\) - n \(∑x\_i\^2\)\/n\)\) \(∑x\_i\^2\)\)\/\(n \(∑x\_i\)\) = \(1\/n\) \((\(x\ y\)\&_\/x\&_ - 
                p \( x\^2\)\&_\/x\&_)\)\)\)
tehát megkaptuk az általános megoldást a lehetséges minimumhelyekre.

A minimum létezéséhez még szükséges az is, hogy az (1.2)
 (({{(∂_11) Q, (∂_12) Q}, {(∂_21) Q, (∂_22) Q}}))
mátrix pozitív definit legyen, ami a mérési pontok ésszerű megválasztása esetén teljesül, ugyanis fenti mátrix éppen a
(({{\(\(x\^2\)\&_\), \(x\&_\)}, {\(x\&_\), 1}}))(({{p}, {q}}))=(({{\(\(x\ y\)\&_\)}, {\(y\&_\)}})) egyenlet együtthatómátrixa, amelynek determinánsa pontosan akkor zérus, ha az (Subscript[x, i]) pontok egyenlők egymással, ami mérési szempontból értéktelen, minden más esetben pedig pozitív (Cauchy-egyenlőtlenség).

Tekintsünk akkor most egy konkrét példát!
Legyen n = 5, x=(1, 4, 5, 6, 9), y=(2.501, 3.982, 4.530, 4.991, 6.497), a keresett egyenes egyenlete pedig y = 0.5  x + 2!

ClearAll[X, Y, p, q, n];
n := 5
X := {1, 4, 5, 6, 9}
Y := {2.501, 3.982, 4.530, 4.991, 6.497}

Minimize[(∑_(i = 1)^n)(p X[[i]] + q - 
              Y[[i]])^2, {p, q}]

Minimize[(p X + q - Y) . (p X + q - Y), {p, q}]

És lám,megkaptuk a minimumot is,ami persze ideális esetben 0 kellene hogy legyen,de mint látható,a kapott eredmény hibája a mérési hibák nagyságrendjében mozog.
Ábrázoljuk a mért pontokat, az eredményül kapott illesztett egyenest és a reziduumokat!

ListPlot[Transpose[{X, Y}], 
      PlotStyle -> PointSize[0.015]] 

    Plot[0.499794 x + 2.00123, {x, 2, 8}] 

    R := 0.499794 X + 2.00123 - Y

    ListPlot[Transpose[{X, R}], 
      PlotStyle -> PointSize[0.015]]

Amint látjuk, a reziduuomok 0 környékén véletlenszerűen ingadoznak.

2. Paraméterekben lineáris, változókban nem lineáris eset
Itt például a \(TraditionalForm\`Q(p, 
          q) = ∑\+\(i = 1\)\%n\((p\ x\_i + q\ \ x\_i\^2 - y\_i)\)\^2\){(Subscript[x, i]),(Subscript[y, i]) adott} függvénnyel kell az előbbi procedúrát eljátszani. Az általános eredményt ezúttal mellőzzük, tekintve hogy ilyen alakú függvény nem olyan gyakori mint az előző.

Nezzunk erre is peldat!
Legyen n = 5, x=(1, 2, 5, 12, 14), y=(7.501, 19.992, 87.502, 419.997, 560.001), a keresett osszefugges pedig y = 5  x + 2.5  (x^2)!

ClearAll[X, Y, p, q, n];
n := 5
X := {1, 2, 5, 12, 14}
Y := {7.501, 19.992, 87.502, 419.997, 560.001}
Minimize[(∑_(i = 1)^n)(p X[[i]] + q X[[i]]^2 - 
              Y[[i]])^2, {p, q}]

Minimize[(p X + q X X - Y) . (p X + q X X - Y), {p, 
        q}]

Az eredmény ezúttal pontosabb, mivel kisebb hibával, és szélesebb tartományban mértünk.
Ábrázolva a mérési pontokat, a leíró függvényt és a reziduumokat:

ListPlot[Transpose[{X, Y}], 
      PlotStyle -> PointSize[0.015]] 

    Plot[4.99932 x + 2.50005 x^2, {x, 0, 15}] 

    R := 4.99932 X + 2.50005 (X^2) - Y

    ListPlot[Transpose[{X, R}], 
      PlotStyle -> PointSize[0.015]]

Amint látható, jelen esetben is szép, véletlen szórást mutató reziduumokat kaptunk.

3. Paraméterekben nem lineáris, változókban nem lineáris, de még linearizálható  eset
Itt a példafüggvényünk a következő: \(TraditionalForm\`Q(p, 
          q) = ∑\+\(i = 1\)\%n\((p\ e\^\(\(\ \)\(q\ x\_i\)\) - y\_i)\)\^2\){(Subscript[x, i]),(Subscript[y, i]) adott}. Az általános eredmény p-re a következő:

\(TraditionalForm\`p = \(∑\(y\_i\) e\^\(q\ x\_i\)\)\/\(∑e\^\(2  q\ x\_i\)\)\)

Megjegyzendő, hogy fenti esetben (lévén csak egy exponenciális tagunk) az egyenlet még linearizálható:
R(p, 
          q) = (∑_(i = 1)^n)(((ln p + q (Subscript[x, i]) - ln (Subscript[y, i])))^2) lenne az új minimalizálandó függvény.

ClearAll[Q, X, Y, p, q]
n := 3
Q[p, 
        q] := (∑_(i = 1)^n)(p E^(q X[[i]]) - 
              Y[[i]])^2
Solve[{D[Q[p, q], p] == 0, D[Q[p, q], q] == 0}, {p, 
        q}]

Az első hibaüzenetek mindössze arra figyelmeztetnek, hogy X és Y még üres, ami rendben is van, elvégre mi ürítettük ki, viszont a Solve is azt mondja, hogy ne algebrai úton próbáljuk az egyenletet megoldani.

Nézzünk erre is példát!
Legyen n = 5, x=(1, 3, 9, 21, 44, ), y=(5.1522, 5.4711, 6.5498, 9.3880, 18.7170), a keresett összefüggés pedig y = 5 (e^(0.03 x))!

ClearAll[X, Y, p, q, n];
n := 5
X := {1, 3, 9, 21, 44}
Y := {5.1522, 5.4711, 6.5498, 9.3880, 18.7170}
Minimize[(∑_(i = 1)^n)(p E^(q X[[i]]) - Y[[i]])^2, {p, q}]

Minimize[(p Exp[q X] - Y) . (p Exp[q X] - Y), {p, 
        q}]

Map[ff, X]

Amint látjuk a konkrét esetek jól kezelhetők.
Próbáljuk ki a NMinimize, FindMinimum és NonlinearFit eljárásokat is, majd ábrázoljuk a mérési eredményeket, a leíró függvényt és a reziduumokat!

NMinimize[(∑_(i = 1)^n)(p E^(q X[[i]]) - Y[[i]])^2, {p, q}]

FindMinimum[(∑_(i = 1)^n)(p E^(q X[[i]]) - Y[[i]])^2, {{p, 4.5}, {q, 0.05}}]

data = {{1, 5.1522}, {3, 5.4711}, {9, 6.5498}, {21, 9.3880}, {44, 
            18.7170}};
<< Statistics`NonlinearFit`;
NonlinearFit[data, 
      p E^(q x), {x}, {p, q}]

NonlinearRegress

ListPlot[Transpose[{X, Y}], 
      PlotStyle -> PointSize[0.015]] 
Plot[
      5.00003 E^(0.0299997 x), {x, 0, 
        45}] 
R := 5.00003 E^(0.0299997 X) - 
        Y
ListPlot[Transpose[{X, R}], 
      PlotStyle -> PointSize[0.015]]

Amint látható, ez esetben viszonylag egyenletesen szórnak a reziduumok.

4. Implicite adott illesztendő függvény esete
Ez a helyzet a legbonyolultabb, itt már csak azt tudjuk az y függvényünkről, hogy megoldása a következő kezdetiérték-problémának: 
y(ξ)=η
y'(x)=f(y(x)) k{(Subscript[x, i]),(Subscript[y, i]) adott}. 

Ebből következően (a továbbiakban eltekintünk y,  k vektor, ill. f mátrix voltának jelölésétől):
y(x) - 
          y(ξ) = ((∫_ξ^x) f(
                  y(s)) ds) k, melyet a transzponálttal balról szorozva (hogy kvadratikus mátrixot kapjunk), majd invertálással k-ra kapjuk, hogy:
\(TraditionalForm\`k = \((\((∫\(f(y(s))\) ds)\)\^T\ \(∫\(f(y(s))\) ds\))\)\^\(-1\)\ \(\((∫\(f(
                      y(s))\) ds)\)\^T\) \((y(x) - 
              y(ξ))\)\).

Példának ezúttal vegyünk egy konszekutív reakciót, melynek egyenlete:
X -> Y -> Z, és keressük a hozzá tartozó (Subscript[k, 1])és (Subscript[k, 2]) sebességi állandókat.
(({{\(x\& . \)}, {\(y\& . 

                  z\& . \)}}))=(({{-x}, {x

                      0}})({{0}, {-y

                      y}}))(({{(Subscript[k, 1])}, {(Subscript[k, 2])}}))
A kezdeti feltétel legyen az, hogy kiinduláskor csak X anyagunk van (mégpedig az egyszerűség kedvéért 1 mol):
y(0)=(1, 0, 0)

ClearAll[x, y, z];
F := {{-x[t], 0}, {x[t], -y[t]}, {0, y[t]}};
K := {k1, k2};
LinearSolve[F, {x'[t], y'[t], z'[t]}]

Amint látjuk, a LinearSolve nem boldogul vele közvetlenül. Próbáljuk meg a k-ra kapott képletet alkalmazni, de előtte tegyünk egy kísérletet a Reduce és a Solve algoritmusokkal is!

Solve[Thread[F . K == {x'[t], y'[t], z'[t]}], K]

K := {k1, k2}; 
    Thread[F . K == {x'[t], y'[t], z'[t]}]

Reduce[Thread[F . K == {x'[t], y'[t], z'[t]}], 
      K]

Solve[F . K == {x'[t], y'[t], z'[t]}, 
      K]

ClearAll[F, x, y, z, t, w, V, M, K, a, b, c, f]
F[x_, y_, 
        z_] := {{-x[t], 0}, {x[t], -y[t]}, {0, y[t]}}
M[F_] := 
      Inverse[(Transpose[Integrate[F, t]]) . (Integrate[F, 
                t])] . (Transpose[Integrate[F, t]])
V[x_, y_, z_, w_] := {x[w], y[w], z[w]} - {x[0], y[0], z[0]}
K[x, y, z] := M[F[x, y, z]] . V[x, y, z, w]
K[x, y, z]

Megkaptuk tehat a k értékét adott x,y,z függvényhármasra, a baj csak az, hogy nekünk x,y,z függvényértékeink vannak (mérési adatokból), és nem ismerjük a leíró függvényt. Módosítani kellene tehát eredményünket oly módon, hogy az x,y,z függvények helyett csak az összetartozó w -> (x(w), y(w), z(w))párokat kelljen megadni, amiből a program előallítja a keresett függvényhármas integráljának egy közelítését. 
Tegyük fel, hogy mérjük percenként az X, Y, Z anyagok koncentrációját. (Ezt megtehetjük pl. egy galváncellában az elektromotoros erő mérésével, amiből a Nernst-egyenlet alapján számolunk vissza koncentrációt, vagy akár titrálással, ha rendelkezésünkre áll egy olyan inhibitor, amely a mérés időpontjában a koncentrációt változtató reakciót "befagyasztja", tehát leállítja.)
Vizsgáljuk jelen esetben tíz percen keresztül a reakcióelegyet, és nyerjük a következő koncentrációkat:
[X] = {1, 0.80, 0.70, 0.65, 0.62, 0.59, 0.55, 
          0.53, 0.51, 0.50} 
[Y] =  {0, 0.18, 0.22, 0.25, 0.27, 0.26, 0.25, 
          0.23, 0.20, 0.20}
[Z] =  {0, 0.02, 0.08, 0.10, 0.11, 0.15, 0.20, 
          0.24, 0.29, 0.30}
ahol figyelembe vettük az anyagmérleget.

<< NumericalMath`ListIntegrate`

datX := {1, 0.80, 0.70, 0.65, 0.62, 0.59, 0.55, 0.53, 0.51, 
        0.50}
datY :=  {0, 0.18, 0.22, 0.25, 0.27, 0.26, 0.25, 0.23, 0.20, 
        0.20} 
datZ :=  {0, 0.02, 0.08, 0.10, 0.11, 0.15, 0.20, 0.24, 0.29, 
        0.30}

ClearAll[F, J, x, y, z, t, w, V, M, K, a, b, c, f] 
{x[w], y[w], z[w]} := {0.5, 0.2, 0.3} ;
{x[0], y[0], z[0]} := {1, 0, 0} ;
J := {{-ListIntegrate[datX, 1], 
            0}, {ListIntegrate[datX, 1], -ListIntegrate[datY, 1]}, {0, 
            ListIntegrate[datY, 1]}};
F[x_, y_, 
          z_] := {{-x[t], 0}, {x[t], -y[t]}, {0, y[t]}};
M[F_] := 
        Inverse[(Transpose[J]) . (J)] . (Transpose[
              J]);
V[x_, y_, z_, 
          w_] := {x[w], y[w], z[w]} - {x[0], y[0], z[0]};
K[x, y, z] := 
        M[F[x, y, z]] . V[x, y, z, w];
K[x, y, z]

Megkaptuk tehát a keresett sebességi állandókat, ami fizikai-kémiai szempontból rendkívül hasznos.

Felhasznált irodalom:
	Szili László - Tóth János: Matematika és Mathematica, ELTE Eötvös Kiadó, Bp. 1996.