In der vorigen Zuschrift war die Normierung voellig verhauen; die
zweidimensionale Normalverteilung ist auf 1 normiert, die BMP Daten
natürlich nicht. Also muss man einen Normvorfaktor n0 in die Funktion fb
einführen und diesen in NonlinearFit[] mit anpassen, mit anderen Worten
In[15]:=
(* das Raster ist \[DoubleStruckCapitalZ] \{CirclePlus]
\[DoubleStruckCapitalZ], die (sichtbare) Norm also die Summe der
normierten Werte *)
unorm = Plus @@ Flatten[udata/Max[udata]] // N
Out[16]= 11976.5
und damit sollte der Vorfaktor etwa zwischen 10000 und 14000 liegen,
In[22]:=Clear[nlf0];
nlf0 = NonlinearFit[latticeData, fb[x, y], {x, y}, {{q, -19/20,
19/20}, {n0, 10000, 14000}, \ {\[Sigma]_1, 1, 10}, {\[Sigma]_2, 1, 10},
{\[Mu]_1, 30, 90}, {\[Mu]_2, 30, 90}}, MaxIterations \[Rule] 71]
Out[23]=
0.9574684344360857753911983526`15.476187252228081/
E^(0.5044318607285512755799524979`17.705943674654403*
(0.0001221615646610996069181002703299863`15.653559774527018*
(-61.063304574417741593099627641`15.954589770191 + x)^2 +
0.00002235686965178924909977430389381796`15.47746851547134*
(-61.063304574417741593099627641`15.954589770191 + x)*
(-56.3084309940831257594022233629`15.954589770191 + y) +
0.00011642434863920711193597886071492704`15.653559774527018*
(-56.3084309940831257594022233629`15.954589770191 + y)^2))
das gibt anscheinend einen passablen Fit, wie das Bildchen zeigt.
Gruss
Udo.
Udo und Susanne Krause wrote:
Frohe Weihanchten, Ulenia!
Ich habe folgendes Problem. Ich möchte zu einem Bild eine Bivariate
Gauss Funktion fitten. Leider geht es nicht und ich bin völlig ratlos.
Bitte, hat jemand eine Idee wie ich das Problem lösen kann?
Das koennen Sie durch Lesen der Dokumentation von NonlinearFit und
durch Betrachten Ihres Modells loesen. Das Modell der
zweidimensionalen Normalverteilung mit dem Korrelationskoffizienten q
(-1 <= q <= 1) hat den Vorfaktor
1/(2 Pi Sqrt[1 - q^2] \[Sigma]_1 \[Sigma]_2).
Nun lesen Sie in der Dokumentation von NonlinearFit:
"The simplest way to specify a parameter is as a symbol, assuming a
starting point of 1.0, or as a {symbol, start} pair. "
Bei Ihrem Input fuer NonlinearFit geht es also mit q=1.0 als Startwert
los, es wird sofort durch Null geteilt und die Iteration ist zuende.
Bad luck.
Leider hat auch nicht jeder das application package ImageProcessing,
versuchen wir also, den Tag ohne das zu retten, indem wir annehmen,
das die (normierten) BMP-Datenpunkte direkt gefittet werden sollen:
In[10]:= Clear[u1];
u1 = Import["C:\\Udo\\Abt_N\\test\\ulenia.bmp"];
Show[u1]
Davon besorgt man sich die Daten (udata)
In[13]:= udata = u1[[1, 1]];
schaut sich das Datenhaeufchen auch einmal an
In[16]:= ListPlot3D[udata]
fertigt den Input fuer NonlinearFit (latticeData)
In[19]:= (* represent udata on the implicit integer lattice *)
Remove[latticeData];
latticeData = Flatten[Table[{o1 - 1, o2 - 1, N[udata[[o1,
o2]]/Max[udata], $MachinePrecision]}, {o1, Dimensions[udata][[1]]},
{o2, Dimensions[udata][[2]]}], 1];
uebernimmt das Modell der zweidim. Normalverteilung
In[17]:= Remove[fb];
fb[x_, y_] := (1/(2*Pi*Sqrt[1 - q^2]*Subscript[\[Sigma],
1]*Subscript[\[Sigma], 2]))*
Exp[(-(1/(2*(1 - q^2))))*((x - Subscript[\[Mu],
1])^2/Subscript[\[Sigma], 1]^2 -
(2*q*(x - Subscript[\[Mu], 1])*(y - Subscript[\[Mu],
2]))/(Subscript[\[Sigma], 1]*
Subscript[\[Sigma], 2]) + (y - Subscript[\[Mu],
2])^2/Subscript[\[Sigma], 2]^2)]
gibt sinnvollen Input fuer NonlinearFit
In[31]:= Clear[nlf0];
nlf0 = NonlinearFit[latticeData, fb[x, y], {x, y}, {{q, -19/20, 19/20},
{Subscript[\[Sigma], 1], 1, 10}, {Subscript[\[Sigma], 2], 1, 10},
{Subscript[\[Mu], 1], 30, 90},
{Subscript[\[Mu], 2], 30, 90}}, MaxIterations -> 71]
Out[32]=
0.0244505167323353093677083653`15.330772284013857/
E^(1.6027492682434579028691143233`15.310053000599336*
(0.1415864513351832032958684763`15.653559774527018*
(-47.2732607821201268020163708695`15.954589770191 + x)^2 -
0.1423494226013671832067579281`15.477468515471339*
(-47.2732607821201268020163708695`15.954589770191 + x)*
(-58.9086438485030340680343468127`15.954589770191 + y) +
0.0520018196274585137088455356`15.653559774527018*
(-58.9086438485030340680343468127`15.954589770191 + y)^2))
und schon hat man "die" Lösung. Bingo. Sie muessen also den Paramter q
einschränken, um nicht durch Null zu teilen.
Die Hilfe sagt weiter: "When elements in the parameter are specified
as {symbol, min, max}, the starting parameter values are taken to be
those minimizing the merit function \[Xi]^2 out of the set forming a
2^p factorial design based on the parameter ranges, where p is the
number of parameters. For example, if the parameter list is specified
as {{a, 0, 1}, {b, 0, 3}}, then the value of {a, b} in {{1/3, 1},
{1/3, 2}, {2/3, 1}, {2/3, 2}} that yields the minimum \[Xi]^2 gives
starting values for parameters a and b."
Deshalb waere es gut, wenn Sie nicht gleich zweidimensional beginnen
wuerden, sondern erst einige eindimensionale Gaussfits für - nicht
unbedingt rechtwinklige (-> q) - Schnitte durch Ihre Datenmenge
anfertigen wuerden, damit die Intervallgrenzen für die Parameter
\[Sigma], \[Mu] halbwegs sinnvoll angegeben werden: Das verbessert die
Startwerte, wie in der Hilfe beschrieben. Bei verbesserten Startwerten
konvergiert NonlinearFit schneller. Das könnte eine Rolle spielen,
wenn Sie den Fit routinemaessig ausfuehren wollen, etwa innerhalb
anderer Programme. Fuer dieses Posting wurde einfach ein sinnvoll
scheinender Satz von Grenzen für die Parameter eingetragen.
Gruss
Udo.
