Re: Statistisches

[Martin Haas <Martin_Haas@XXXXXXX.net>]

Hallo Udo!

Vielen Dank für Deine schnelle Antwort auf unsere Frage! Uli
ist gerade auf einer Konferenz, aber ich bin Doktorand in der 
gleichen Gruppe. 

Dein erster Vorschlag für unseren Fehler, dass wir einen
Gauss-Term auf das Modell für den Fit addiert hätten, trifft 
tatsächlich nicht zu (schade ..). Und die Hypothese, ob das 
Rauschen Gauss'scher Natur ist, müssen wir die wirklich testen? 
Denn wir addieren doch gerade Rauschen mit der Statistik 
NormalDistribution. Oder liege ich da falsch?

Unser Problem taucht auch erst auf, wenn man Weights->{Gewichte} 
bei Nonlinear Regress einschaltet. Etwa wie folgt: Wir generieren 
einen Lorentz, wie in Deinem Beispiel, nur fitten wir diesmal
mit NonlinearRegress[ .. , Weights->{Table[1/sigma(n)^2]}]. 

Damit wollen wir der Tatsache Rechnung tragen, dass die Punkte 
verschiedenes statistisches Gewicht haben, denn die Unsicherheit
(Fehlerbalken) sigma sei bei Punkten mit höherem Funktionswert n
grösser:

 sigma[n_]:=Sqrt[n+1] + 0.025 n;

 Lorentz[x_,a_,c_,g_]:=a/((x-c)^2 + (g/2)^2);

 With[{a=20, c=0, g=1},
  data=Table[{x,Lorentz[x,a,c,g] + 
                Random[NormalDistribution[0,sigma[Lorentz[x,a,c,g]]]]
             },{x,-2,2,0.05}];
 ];
 myweights=Table[1/sigma[data[[i,2]]]^2, {i,1,Length[data]}]; 
 NonlinearRegress[data, Lorentz[x1, a1, c1, g1], x1, {a1, c1, g1}, 
 Weights->myweights, RegressionReport -> BestFitParameters]

Out=BestFitParameters -> 
    {a1 -> 14.7161, c1 -> 0.0309276, g1 -> 0.860453}
(vgl Werte bei Generierung:
    {a  -> 20,      c ->  0,         g  -> 1})        
    
Bei jedem Run bekommen wir systematisch zu kleine Amplitude,
und zu kleine Breite. 

Dahingehend nun unsere Frage nach dem grundlegenden Denkfehler 
(oder ist es NonlinearRegress?). Jeder Punkt folgt einer Gaussverteilung 
mit Breite sigma[n], dies soll mit der Gewichtung 1/sigma(n)^2 
(pro Punkt individuell) beim Fit berücksichtigt werden. Das klappt 
aber nicht.

Unser Verdacht ist nun, dass das Problem darin besteht: Wir weisen jedem 
Punkt sein Gewicht auf der Basis seines tatsächlichen Wertes zu (wie bei 
gemessenen Daten nunmal nicht anders möglich), dieser ist mal grösser, mal 
kleiner als die Modellfunktion. Dabei bekommen jedoch bei gegebenem x die 
kleineren Realisierungen n systematisch grössere Gewichte. Das könnte die 
Amplitude sowie die Breite verringern. Wenn das so wäre, wie bezieht man die 
verschiedenen Gewichte korrekt ein?

Viele Grüße,

Martin Haas

On Sunday 01 August 2004 16:53, you wrote:
> 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.
>