Re: Statistisches

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

Hallo Uli,

was genau mit den Gewichten passiert ist in Mathematica bei Dir, kann 
man mangels Beispiel nicht nachvollziehen.
Nur zum Sagen, auf eine Idee wie die folgende wirst Du als Physiker 
sicher nicht gekommen sein:

In[1]:= <<Statistics`NonlinearFit`
In[2]:= <<Statistics`ContinuousDistributions`

In[8]:= With[{? = 0., ? = 0.75, ? = 0.25},
pic = Plot[1/(x^2 + (?/2)^2) + Random[NormalDistribution[?, ?]], {x, -2, 
2},
PlotDivision -> 30, DisplayFunction -> Identity]; data = pic[[1,1,1,1]];
Print["There are ", Length[data], " datapoints."];
NonlinearRegress[data, 1/(x1^2 + (?1/2)^2) + Exp[-((x1 - ?1)^2/(2*?1^2))]/
(Sqrt[2*Pi]*?1), {x1}, {?1, ?1, ?1}, RegressionReport ->
BestFitParameters]]

"There are "860" datapoints."
Out[8]=
{BestFitParameters -> {?1 -> 0.0037703537868760337, ?1 -> 
0.07695177010624978,
?1 -> 0.26111019517895107}}

Dieser Ansatz für die Fitfunktion ist falsch, weil (1) 
NormalDistribution[my, sigma] positive und negative Werte abgibt, aber 
die Normalverteilung im Ansatz bekanntermassen nicht; (2) die 
funktionalen Abhängigkeiten verkehrt sind, das Rauschen ist unabhänhig 
von x1. Mit dem Ansatz erzeugt man nur systematische Fehler.

Man kann testen, inwiefern die Lorentzdaten von einer Normalverteilung 
verrauscht wurden, das ist ein Test von 2 Hypothesen; Fitter testen 
immer Hypothesen in Form der Ansätze, in dem Fall habe ich keinen Weg 
gefunden, beide Hypthesen mit einem Aufruf von NonlinearRegression zu 
testen:

In[161]:= Clear[data, ?1, ?2, gauss];
With[{? = 0., ? = 0.75, ? = 0.25, dy = 0.1},
pic = Plot[1/(x^2 + (?/2)^2) + Random[NormalDistribution[?, ?]], {x, -2, 
2},
PlotDivision -> 30, DisplayFunction -> Identity]; data = pic[[1,1,1,1]];
Print["There are ", Length[data], " data points."];
?2 = NonlinearRegress[data, 1/(x1^2 + (?1/2)^2), {x1}, {?1},
RegressionReport -> BestFitParameters][[1,2,1,2]];
Print["LorentzFit: ?2 = ", ?2];
dist = Sort[Last[Transpose[data]] - (1/(#1^2 + (?2/2)^2) & ) /@
First[Transpose[data]]]; val = First[dist] + dy/2.; oo = 0; gauss = {};
For[o = 1, o <= Length[dist], o++, If[dist[[o]] < val, oo++,
AppendTo[gauss, {val - dy/2., oo}]; val += dy; oo = 0; --o]];
gauss = Transpose[{First[Transpose[gauss]], Last[Transpose[gauss]]/
(dy*Length[dist])}];
Print["Gaussian Fit:"];
NonlinearRegress[gauss, Exp[-((x1 - ?1)^2/(2*?1^2))]/(Sqrt[2*Pi]*?1), 
{x1}, {?1, ?1},
RegressionReport -> BestFitParameters]
]

From In[161]:= "There are "860" data points."
From In[161]:= "LorentzFit: ?2 = "0.25040192109887016
From In[161]:= "Gaussian Fit:"

Out[163]= {BestFitParameters -> {?1 -> -0.0067614107281022936, ?1 -> 
0.7576946340056632}}

Also ganz gut getroffen diesmal. Andere Runs geben auch schon wesentlich 
schlechtere Ergebnisse. Die Ansicht
ListPlot[gauss] ist instruktiv in solchen Fällen.

Schönen Sonntag
Udo.

Uli Jentschura wrote:

> Lieber DEMUG-Mitglieder,
>
> hier eine Frage betreffs
>
>  Statistik,      insbesondere      NonlinearRegress.
>
> Bevor ich das Problem im Detail schildere, moechte ich mich vergewissern,
> nicht von Vornherein einen Denkfehler gemacht zu haben:
>
> Nehmen wir folgendes an:
>
> 1. Wir nehmen eine wohldefinierte Kurve, zum Beispiel eine Lorentzkurve
>   der Form
>
>
>      f(x) = 1/(x^2 + gamma^2/4)
>
>   (bis auf Normierung).
>
> 2. Wir addieren ein statistisches Rauschen, mit einer Gauss-schen
>   Verteilung ("Random" mit "NormalDistribution"), so dass die Linie
>   varrauscht wird, wobei sich die Standardabweichung der addierten
>   Zufallszahl ergibt zu
>
>     sigma(f) = Sqrt[f] + constant * f    (*)
>
>   (man/frau koennte sich auch andere Abhaengigkeiten sigma(f)
>    vorstellen.)
>
> 3. Wir fitten mit "NonlinearRegress" die entstandene verrauschte Linie,
>   wobei wir die verrauschten Daten statistisch gewichten, und zwar genau
>   mit demjenigen Gewicht 1/sigma(f)^2, das wir vorher benutzt haben, um
>   die "Breite" des statistischen Rauschens festzulegen.
>
> Ueberlegung (richtig oder falsch?):
>
> Dann muesste doch im statistischen Mittel eine nichtlineare Anpassung an
> die verrauschte Linie (unter Zuhilfenahme von "NonlinearRegress") mit
> einem Modell f(x) = 1/(x^2 + gamma^2/4) unter Gewichtung der f_i mit
> Gleichung (*) fuer den Parameter gamma den "Startwert" von gamma ergeben.
>
> Beobachtung:
>
> Die Breite gamma der Lorentzkurve wird nicht richtig gefittet. Vielmehr
> ist diese Breite typischerweise kleiner als gamma. Aber vielleicht liegt
> ein Denkfehler vor. (Hoffentlich liegt ein Denkfehler vor!)
>
> Vielen Dank an alle Diskussionsteilnehmer im voraus!
>
> Gruss,
>
> Ulrich Jentschura
>
> ------------------------------------------------------------------------------
> Ulrich Jentschura                          Theoretische Quantendynamik/ Physik
> Albert-Ludwigs-Universitaet Freiburg, Hermann-Herder-Strasse 3, 79104 Freiburg
> Tel. +49-761-203 5956 (fax extension 5883)   jentschura@XXXXXXX.de
> internet homepage: http://tqd1.physik.uni-freiburg.de/~ulj
>
>
>
>