Re: Nonlinear Fit zu Bivariate Normal Distribution

[Udo und Susanne Krause <su.krause@XXXXXXX.ch>]

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.