Re: Statistisches

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

Hallo Martin,

Das Ganze ist eine Definitionssache. Geschrieben steht (Mma 4.2):
The Weights option allows you to implement weighted least squares by 
specifying a list of weights, one for each data point; the default 
Weights -> Automatic implies a weight of unity for each data point. When 
Weights -> {w_1, w_2, ..., w_N} the parameter estimates are chosen to 
minimize the weighted sum of squared residuals Sum[w_i (e_i)^2, {i, 1, N].

Also nicht die Gewichte werden quadriert, sondern die Abweichungen von 
den Messpunkten. Die Gewichte stehen linear drin (soetwas wie eine 
Massendichte, deswegen heissen die auch Gewichte). Tippt man etwa

In[67]:=
With[{a = 20., c = 0., g = 1., xL = -2., xR = 2., anZ = 100},
 data = ({#, Lorentz[#, a, c, g] + Random[NormalDistribution[0., 
sigma[Lorentz[#, a, c, g]]]]} & ) /@ Sort[Table[Random[Real, {xL, xR}, 
8], {anZ}]];
 ListPlot[data];
 Table[{j, NonlinearRegress[data, Lorentz[x1, a1, c1, g1], {x1}, {a1, 
c1, g1},
   Weights -> (sigma[#]^j &), RegressionReport -> BestFitParameters]},
   {j, -2, 2, 1/2}
 ]
]

dann bekommt man neben dem überzeugend Verrauschung zeigenden (hier 
unterdrückten) Bildchen das Resultat:

Out[67]=
{{-2, {BestFitParameters -> {a1 -> 18.4578, c1 -> 0.00484121, g1 -> 
0.962976}}},
 {-3/2, {BestFitParameters -> {a1 -> 18.882, c1 -> -0.000956114, g1 -> 
0.971776}}},
 {-1, {BestFitParameters -> {a1 -> 19.2635, c1 -> -0.00574937, g1 -> 
0.97938}}},
 {-1/2, {BestFitParameters -> {a1 -> 19.6093, c1 -> -0.00975918, g1 -> 
0.98604}}},
 {0, {BestFitParameters -> {a1 -> 19.9292, c1 -> -0.0131599, g1 -> 
0.992069}}},
 {1/2, {BestFitParameters -> {a1 -> 20.2349, c1 -> -0.0160729, g1 -> 
0.997817}}},
 {1, {BestFitParameters -> {a1 -> 20.5375, c1 -> -0.0185851, g1 -> 
1.00359}}},
 {1/2, {BestFitParameters -> {a1 -> 20.845, c1 -> -0.0207461, g1 -> 
1.00958}}},
 {2, {BestFitParameters -> {a1 -> 21.1619, c1 -> -0.0225848, g1 -> 
1.01589}}}}

Also für Gewichte konstant 1 (Weights -> Automatic, j = 0) ist das 
Ergebnis normalerweise am besten.Wenn die Gewichte 1/sigma[#]^2  (j = 
-2) sind, dann werden gerade bei den grossen Abweichungen die Gewichte 
klein, d.h. die Abweichungen um x = 0 werden unwichtiger, der Fit wird 
am Rand besser. Werden die Gewichte bei den grossen Abweichungen um x = 
0 sehr gross genommen (j = 2), dann bemüht er sich übermässig, bei x = c 
= 0 richtig zu liegen. Man sieht ganz klar die Tendenz. Eure Überlegung 
wäre korrekt gewesen, wenn NonlinearRegress[] den Ausdruck 
Sum[(e_i/w_i)^2, {i, 1, N}] minimieren würde. Das macht es aber lt. 
Beschreibung und lt. numerischem Experiment nicht.

Mit den besten Grüssen
Udo.

Martin Haas wrote:

<snip>

> 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.
> >    
<snip>

> > 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.