Hallo Martin,
das beiliegende Bildchen ist für die Fitparameter
Weight: 0 BestFitParameters: a = 19.0559 c = 0.0201002 g = 0.92473
das sieht zahlenmässig schlecht aus und spricht anscheinend gegen die
Weights -> Automatic. Deshalb habe ich die Bildchen nun doch erzeugt.
Der Fit ist aber gut. Das Problem ist eher, dass im Bereich des Maximums
der Ausgangslorentzkurve nur sehr wenige Punkte zu liegen kommen. Wenn
diese alle auf dieselbe Seite streuen, verzerrt es den Fit ganz klar, i.
e. die Startwerte können nicht getroffen werden.
Mit den besten Grüssen
Udo.
In[43]:=
With[{a = 20., c = 0., g = 1., xL = -2., xR = 2., anZ = 100, j = 0},
data = ({#, Lorentz[#, a, c, g] + Random[NormalDistribution[0.,
sigma[Lorentz[#, a, c, g]]]]} & ) /@ Sort[Table[Random[Real, {xL,
xR}, 8], {anZ}]];
If[Length[Select[Last[Transpose[data]], Negative]] > 0,
Print["Negative data point(s) generated."];
Return[$Failed]
];
pic = ListPlot[data, DisplayFunction -> Identity];
{a2, c2, g2} = {a1, c1, g1} /. BestFitParameters /.
NonlinearRegress[data, Lorentz[x1, a1, c1, g1], {x1}, {a1, c1, g1},
Weights -> (sigma[#]^j &), RegressionReport -> BestFitParameters];
Print["Weight: ", j, " BestFitParameters: a = ", a2, " c = ", c2, "
g = ", g2];
Show[{pic,
Plot[Lorentz[x2, a2, c2, g2], {x2, xL, xR}, PlotStyle -> RGBColor[1,
0, 0],
DisplayFunction -> Identity],
Plot[Lorentz[x2, a, c, g], {x2, xL, xR}, PlotStyle -> RGBColor[0, 0,
1],
DisplayFunction -> Identity]}, DisplayFunction -> $DisplayFunction]
]
Udo und Susanne Krause wrote:
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.
