5.4. Mathematik der linearen LSQ-Kurvenanpassung mit Messwerten, die miteinander korreliert sind¶
5.4.1. Basics¶
Das Auswertemodell \(\mathbf{M}\) wird in folgender Form formuliert (entsprechend Weise & Wöger; Hinweis: Die Bedeutung von X und Y ist im Vergleich zu Weise & Wöger gerade vertauscht):
\(\mathbf{M}\left( \mathbf{X,\ Y} \right)\mathbf{= X - A\ Y = 0}\)
Der nr-Vektor \(\mathbf{Y}\) charakterisiert die zu bestimmenden Parameter, der n-Vektor \(\mathbf{X}\) die eingehenden Messwerte, welche die Fitkurve repräsentieren. Die (n x nr)-Matrix \(\mathbf{A}\) enthält die partiellen Ableitungen der Fitfunktion, nach den zu bestimmenden Parametern, das sind die Funktionen \(X_{i}(t)\). Die Varianzen sowie Kovarianzen zwischen den gemessenen x-Werten seien in der nicht-diagonalen Unsicherheitsmatrix \(\mathbf{U}_{\mathbf{x}}\) zusammengestellt.
Nach Weise & Wöger erhält man für die Lösungen folgende Gleichungen (siehe auch DIN ISO 11929:2010, Anhang C.5.4, Gl. (C.32); gleiche Notation wie hier!):
Hierin ist \(\mathbf{y}\) der Ergebnisvektor der angepassten Fitparameter und \(\mathbf{U}_{\mathbf{y}}^{\mathbf{- 1}}\) die dazugehörige Kovarianzmatrix.
Für den Wert der Minimumfunktion erhält man unter Verwendung von \(\mathbf{\ A}^{\mathbf{T}}\mathbf{\ U}_{\mathbf{x}}^{\mathbf{- 1}}\mathbf{\ x =}\mathbf{U}_{\mathbf{y}}^{\mathbf{- 1}}\mathbf{\ y}\):
\(\mathbf{\min}\mathbf{\chi}^{\mathbf{2}}\mathbf{=}\left( \mathbf{Ay - x} \right)^{\mathbf{T}}\mathbf{\ }\mathbf{U}_{\mathbf{x}}^{\mathbf{- 1}}\mathbf{\ }\left( \mathbf{Ay - x} \right)\mathbf{=}\mathbf{x}^{\mathbf{T}}\mathbf{\ }\mathbf{U}_{\mathbf{x}}^{\mathbf{- 1}}\mathbf{\ x -}\mathbf{y}^{\mathbf{T}}\mathbf{\ }\mathbf{U}_{\mathbf{y}}^{\mathbf{- 1}}\mathbf{\ y}\)
\(\chi_{R}^{2} = \min{\chi^{2}/(n - nr)}\)
Die Kovarianzmatrix \(\mathbf{U}_{\mathbf{y}}^{\mathbf{- 1}}\) kann mit \(\chi_{R}^{2}\) multipliziert werden, falls dessen Wert größer als eins ist, obwohl dies bei strikter Befolgung der Bayes-Statistik nicht erlaubt ist und daher in UncertRadio nicht gemacht wird.
5.4.2. Chi-Quadrat-Optionen¶
Hierzu kann folgende Literatur empfohlen werden:
S. Baker, R.D. Cousins, 1983: Clarifications of the use of Ch-square and likelihood functions. Nucl. Instr. & Meth 221, S 437-442.
T. Hauschild, M. Jentschel, 2001: Comparison of maximum likelihood estimation and chi-square statistics applied to counting experiments. Nucl. Instr. & Meth A 457 (1-2), S 384-401.
Die bis hierher durchgeführte Rechnung arbeitet mit dem sogenannten „Neyman Chi-square”-Ausdruck, das wie folgt definiert ist (Verfahrens-Kürzel: WLS):
\(\chi_{N}^{2} = \sum_{j}^{}\frac{\left( x_{j} - f\left( t_{j};\mathbf{y} \right) \right)^{2}}{u^{2}\left( x_{j} \right)}\) ,
d.h. die Diagonalelemente von \(\mathbf{U}_{\mathbf{x}}\) enthalten Varianzen der gemessenen Werte \(x_{j}\) . Dies kann bei niedrigen Zählraten zu einem Bias im Ergebnis für \(\mathbf{y}\) führen. Zur Reduktion des Bias kann das sogenannte „Pearson Chi-square“-Ausdruck herangezogen werden, das definiert ist als (Verfahrens-Kürzel: PLSQ):
\(\chi_{P}^{2} = \sum_{j}^{}\frac{\left( x_{j} - f\left( t_{j}\mathbf{;y} \right) \right)^{2}}{{\widetilde{u}}^{2}\left( f\left( t_{j};\mathbf{y} \right) \right)}\)
Die im Nenner stehende Varianz wird mit Hilfe der anzupassenden Funktion \(f\) berechnet.
\({\widetilde{\mathbf{U}}}_{\mathbf{x}}(f) = diag\left\lbrack \frac{\left( \mathbf{Ay} + R_{0,i} + R_{bl} \right)}{t_{m}} + var\left( R_{0,i} \right) + var\left( R_{bl} \right) \right\rbrack\) ,
worin \(R_{0,i}\) und \(R_{bl}\) die Nulleffekt- und die Netto-Blindwert-Zählrate darstellen.
Da dieser Funktionswert nur mit Hilfe der dabei anzupassenden Parameter y berechnet werden kann, ist zur Minimierung dieses Chi-Quadrat-Ausdrucks ein iteratives Vorgehen erforderlich. Die Iteration lässt sich in der zu Gl. (5.4.1) analogen Schreibweise, mit (m) und (m‑1) als aufeinander folgender Iterationsschritte, wie folgt ausführen (m=3):
\(\mathbf{y}_{\mathbf{(m)}} = \mathbf{U}_{\mathbf{y,(m)}}\mathbf{\ }\mathbf{A}^{T}\mathbf{U}_{\mathbf{x,(m - 1)}}^{\mathbf{- 1}}\left( \mathbf{x} \right)\mathbf{\ x}\) ; \(\mathbf{\ }\mathbf{U}_{\mathbf{y,(m)}} = \left( \mathbf{A}^{T}\mathbf{U}_{\mathbf{x,(m - 1)}}^{\mathbf{- 1}}\left( \mathbf{x} \right)\mathbf{A} \right)^{- 1}\)
Hierzu werden die „gemessenen” Varianzen vor der Wiederholung der Rechnung durch die Varianzen der zuvor angepassten Funktion ersetzt, indem die diagonalen Elemente der Unsicherheitsmatrix \(\mathbf{U}_{\mathbf{x}}\) der Messwerte durch die mit Hilfe der aus dem ersten Auswertelauf (Neyman-Chi-Quadrat) erhaltenen angepassten Werte \(\mathbf{Ay}\) neu berechnet werden (die nicht-diagonalen Elemente bleiben unverändert):
\(\mathbf{U}_{\mathbf{x}}^{\mathbf{'}} = diag\left\lbrack \frac{\left( \mathbf{Ay} + R_{0,i} + R_{bl} \right)}{t_{m}} + var\left( R_{0,i} \right) + var\left( R_{bl} \right) \right\rbrack\)
Anmerkung: Sind die Kovarianzen zwischen den Eingangswerten gleich null, d.h. die Matrix \(\mathbf{U}_{\mathbf{x}}\) diagonal mit den Varianzen als Elementen, entspricht das Verfahren mit dem Neyman-Chiquadrat-Ausdruck vollständig der bekannteren Methode der gewichteten linearen Least-squares-Anpassung.
Alternativ zum Pearson-Chi-Quadrat-Ausdruck kann zur Verringerung des Bias auch das sogenannte „Poisson Maximum Likelihood Estimation”-Verfahren (oder Poisson MLE) verwendet werden. Der dazugehörige Chi-Quadrat-Ausdruck lautet (Verfahrens-Kürzel: PMLE):
\(\chi_{\lambda}^{2} = 2\left\lbrack \sum_{j}^{}\left( f\left( t_{j};\ \mathbf{y} \right) - x_{j} \right) - \sum_{j;x_{j} \neq 0\ }^{}{x_{j}\ln\left( \frac{f\left( t_{j};\ \mathbf{y} \right)}{x_{j}} \right)} \right\rbrack\) .
Die Minimierung dieses Ausdrucks erfordert ein nicht-lineares Anpassungsverfahren: In UncertRadio wird eine modifizierte Version der Levenberg-Marquardt-Anpassungsroutine LSQmar aus dem Lehrbuch „Datenanalyse” von S. Brandt verwendet. Die für Poisson MLE notwendige Modifikation von LSQmar kann dem Bericht LLNL-JRNL-420247 (2009) entnommen werden.
Werden mit den anderen Fitverfahren Nettozählraten angepasst, so erfordert das PMLE-Verfahren die Anpassung an Brutto-Impulsanzahlen. Dies wird intern in UncertRadio berücksichtigt. Dies hat den Vorteil, dass es zwischen Bruttoimpulsanzahlen keine Kovarianzen gibt, was bei der Verwendung von Nettozählraten jedoch der Fall ist.
Das PMLE-Verfahren kann in dem Beispiel-Projekt vTI-Y90-16748_Blw_V2_DE.txp aktiviert werden.
Nach Baker & Cousins bleibt bei Anwendung der PMLE-Methode die Fläche unter der anzupassenden Kurve erhalten, während PLSQ zu Überschätzung und WLS zu Unterschätzung neigt.
Statistische Tests der Chi-Quadrat-Optionen
Es wird die Summe zweier gleichhoher Gauß-Peaks mit identischer Breite (sigma = 5,1 Kanäle, Nettopeakfläche je 150 Impulse, bei Kanal 30 und Kanal 90) auf einem konstanten Untergrund von 4 Impulsen/Kanal als die anzupassende Funktion verwendet. Die Fitfunktion hat also 6 Parameter: 4.0, 150, 30, 5.1, 150, 90. Diese Parameter werden als wahre Werte vorgegeben, um danach eine große Anzahl Nsp=400 entsprechender Spektren zu erzeugen deren „wahre” Kanalinhalte als Poisson-verteilte Werte statistisch verrauscht werden. Jedes dieser Spektren wird mit den verschiedenen oben angeführten Optionen angepasst. Es wird sowohl ein nicht-lineares als auch ein lineares Fitverfahren verwendet. Die nachfolgende Abbildung zeigt ein Spektrum mit der „wahren” Fitfunktion und den mit nicht-linearem Fitten erhaltenen Kurven für die Verfahren WLS und PMLE.
Die Auswertung von jeweils 400 Spektren, die nicht-linear mit den Verfahren WLS, PLSQ und PMLE angepasst werden, ist im Folgenden tabelliert. Aus jeweils 400 angepassten Werten der 6 Parameter werden ein Mittelwert (pa=), eine Standardabweichung zwischen den 400 Werten (sd=) und der Mittelwert der von der Anpassungsroutine berechneten individuellen Standardunsicherheiten (mean(sds..)) abgeleitet. Zusätzlich werden Summen (über die 6 Parameter) der absoluten Abweichungen (Fitwert minus wahrer Wert) und der Differenzen (sd minus sds) berechnet, die in zwei zusätzlichen Zeilen unter jeder Einzeltabelle angegeben sind. Diese Werte dienen als gute Indikatoren für die Konsistenz der Anpassung.
Nichtlineare Anpassung:
mean values: for method WLS
i=1 pa= 2.75534E+00 sd(pa(i))= 2.34472E-01 mean(sds(i))= 1.73935E-01i=2 pa= 1.57094E+02 sd(pa(i))= 1.91809E+01 mean(sds(i))= 1.59540E+01i=3 pa= 2.99839E+01 sd(pa(i))= 7.59367E-01 mean(sds(i))= 5.81371E-01i=4 pa= 5.18531E+00 sd(pa(i))= 5.09876E-01 mean(sds(i))= 3.88957E-01i=5 pa= 1.58079E+02 sd(pa(i))= 1.84934E+01 mean(sds(i))= 1.59980E+01i=6 pa= 8.99873E+01 sd(pa(i))= 7.84781E-01 mean(sds(i))= 5.79122E-01sum of absolute deviations from true(pa) : 16.5316238sum of absolute deviations of two sd types : 6.28735828mean values: for method PLSQ
i=1 pa= 3.99617E+00 sd(pa(i))= 2.11079E-01 mean(sds(i))= 1.73482E-01i=2 pa= 1.49851E+02 sd(pa(i))= 1.75574E+01 mean(sds(i))= 1.57870E+01i=3 pa= 3.00262E+01 sd(pa(i))= 6.38806E-01 mean(sds(i))= 5.91121E-01i=4 pa= 5.07812E+00 sd(pa(i))= 4.50239E-01 mean(sds(i))= 3.95102E-01i=5 pa= 1.50346E+02 sd(pa(i))= 1.63419E+01 mean(sds(i))= 1.58062E+01i=6 pa= 8.99958E+01 sd(pa(i))= 6.39568E-01 mean(sds(i))= 5.90819E-01sum of absolute deviations from true(pa) : 0.551150203sum of absolute deviations of two sd types : 2.49521804mean values: for method PMLE
i=1 pa= 4.00893E+00 sd(pa(i))= 2.10194E-01 mean(sds(i))= 2.07077E-01i=2 pa= 1.51928E+02 sd(pa(i))= 1.69198E+01 mean(sds(i))= 1.68943E+01i=3 pa= 2.99583E+01 sd(pa(i))= 6.46945E-01 mean(sds(i))= 6.26552E-01i=4 pa= 5.08396E+00 sd(pa(i))= 4.22725E-01 mean(sds(i))= 4.34723E-01i=5 pa= 1.49174E+02 sd(pa(i))= 1.73843E+01 mean(sds(i))= 1.67870E+01i=6 pa= 8.99848E+01 sd(pa(i))= 6.22015E-01 mean(sds(i))= 6.35672E-01sum of absolute deviations from true(pa) : 2.83508420sum of absolute deviations of two sd types : 0.671966136Die beste absolute Abweichung (Fitwert minus wahrer Wert) wird für das Verfahren PLSQ erhalten; die beste Konsistenz der Unsicherheiten (abs. Abweichung (sd minus sds) ) hat sich für das PMLE-Verfahren ergeben.
Dass in diesem Beispiel die Verfahren PLSQ und PMLE deutlich besser als der klassische Fit (WLS) abschneiden, ist damit begründet, dass der konstante Untergrund mit 4 Impulsen/Kanal sehr niedrig ist. Wenn der Untergrund deutlich zunimmt, fallen die Unterschiede zwischen den Fitverfahren geringer aus.
Die Spektren können auch mit linearen Anpassungsmethoden behandelt werden, wenn der Breitenparameter und die beiden Peakpositions-Parameter auf ihren wahren Werten fixiert gehalten werden. In diesem Fall werden nur der Untergrund-Parameter und die beiden Peakflächen angepasst. Die Ergebnisse einer solchen Auswertung sind im Folgenden angegeben.
Lineare Anpassung:
mean values: method = WLS
i=1 pa= 2.77091E+00 sd(pa(i))= 2.19480E-01 mean(sds(i))= 1.61749E-01i=2 pa= 1.53462E+02 sd(pa(i))= 1.78632E+01 mean(sds(i))= 1.50700E+01i=5 pa= 1.54138E+02 sd(pa(i))= 1.74788E+01 mean(sds(i))= 1.51000E+01sum of absolute deviations from true(pa) : 8.82905006sum of absolute deviations of two sd types : 5.22971344mean values: method = PLSQ
i=1 pa= 4.00339E+00 sd(pa(i))= 1.95047E-01 mean(sds(i))= 1.92895E-01i=2 pa= 1.49325E+02 sd(pa(i))= 1.63338E+01 mean(sds(i))= 1.58374E+01i=5 pa= 1.49789E+02 sd(pa(i))= 1.53364E+01 mean(sds(i))= 1.58544E+01sum of absolute deviations from true(pa) : 0.889460266sum of absolute deviations of two sd types : 1.01658702Ebenso zeigen wie beim nichtlinearen Fit die Ergebnisse für die PLSQ-Methode eine bessere Konsistenz als für die WLS-Methode.
5.4.3. Export der Eingangsdaten nach R¶
Seit dieser Version von UncertRadio ist es möglich, die Eingangsdaten, abhängig vom gewählten Anpassungsverfahren (außer WTLS: R unterstützt dies noch nicht!), in einem Format, das mit der entsprechenden R-Routine kompatibel ist, in eine Textdatei zu exportieren (URExport-to-R.txt oder ähnlich). Diese Datei(en) können leicht vom Statistikpaket R importiert werden. Dies erlaubt einen Vergleich der Ergebnisse zwischen UR und R. Diese Option wird über den Menüpunkt „options – LSQ export to R” aufgerufen und kann für die Fälle Berechnung der Ergebnisgröße und der Erkennungsgrenze verwendet werden. Die dafür benötigten Daten beziehen sich nur auf Zählraten, nicht auf das Ergebnis der Ergebnisgröße. Für den Fall der Ergebnisgröße (siehe URExport-to-R.txt) erhält man:
(Hinweis: die Kovarianz-Matrix ist hier rechts abgeschnitten worden)
Fall: Ergebnisgröße
Blank count rate= 4.66670009E-08 background rate= 1.88333332E-03
Input data: variance-covariance matrix: (rank= 18 )
2.87780E-07 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-082.61574E-08 2.46788E-07 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-082.61574E-08 2.61574E-08 1.96152E-07 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-082.61574E-08 2.61574E-08 2.61574E-08 1.64805E-07 2.61574E-08 2.61574E-08 2.61574E-082.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 1.41898E-07 2.61574E-08 2.61574E-082.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 1.35870E-07 2.61574E-082.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 1.43104E-072.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-082.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-082.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-082.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-082.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-082.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-082.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-082.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-082.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-082.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-082.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08Arrays y, X1, x2, X3:
y X1 X3 (Eingangsdaten-Matrix)
1 5.65134E-03 8.80220E-01 2.73093E-012 4.47079E-03 8.07210E-01 1.10866E-013 3.01245E-03 7.40257E-01 4.50075E-024 2.10968E-03 6.78857E-01 1.82714E-025 1.44995E-03 6.22550E-01 7.41751E-036 1.27634E-03 5.70913E-01 3.01124E-037 1.48468E-03 5.23559E-01 1.22245E-038 9.98564E-04 4.80133E-01 4.96272E-049 8.47116E-04 4.40700E-01 2.03089E-0410 8.24953E-04 4.03787E-01 8.17887E-0511 1.24162E-03 3.70295E-01 3.32032E-0512 5.12453E-04 3.39582E-01 1.34793E-0513 9.63842E-04 3.11415E-01 5.47210E-0614 3.38842E-04 2.85585E-01 2.22147E-0615 1.65231E-04 2.61897E-01 9.01838E-0716 -7.78244E-05 2.40175E-01 3.66113E-0717 2.69398E-04 2.20253E-01 1.48629E-0718 1.99953E-04 2.01985E-01 6.03378E-08Parameter values and std uncertatinties obtained by UR:
1 2.83190E-03 3.55440E-043 1.45234E-02 2.01819E-03Chisqr= 1.23143363
Für den Import in R werden zusätzlich jeweils die Kovarianz-Matrix und die Eingabedaten in separate Dateien geschrieben:
Ergebnisgröße: covmat1.txt und data1.txt
Erkennungsgrenze: covmat2.txt und data2.txt
Mit diesen Dateien kann eine statistische Auswertung mit R wie folgt erfolgen:
(load package MASS) (R)
One obtains with R for the output quantity:
> covmat <- read.table("covmat1.txt")> data <- read.table("data1.txt")> res <- lm.gls(formula = y ~ X1 + X3 - 1, data = data, W = covmat, inverse = TRUE)> summary.lm(res)Call:
lm.gls(formula = y ~ X1 + X3 - 1, data = data, W = covmat, inverse =TRUE)Residuals:Min 1Q Median 3Q Max-8.076e-04 -4.422e-04 -3.702e-04 -3.134e-05 5.747e-04Coefficients:Estimate Std. Error t value Pr(>|t\|)X1 2.832e-03 1.643e-07 17231 <2e-16 \**\*X3 1.452e-02 9.332e-07 15564 <2e-16 \**\*---Signif. codes: 0 ‘\**\*’ 0.001 ‘\*\*’ 0.01 ‘\*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.0004624 on 16 degrees of freedomMultiple R-squared: 0.9625, Adjusted R-squared: 0.9578F-statistic: 205.1 on 2 and 16 DF, p-value: 3.95e-12Warning message:In summary.lm(res) : calling summary.lm(<fake-lm-object>) ...Um die Unsicherheit mit derjenigen aus UR vergleichen zu können, teilt man die Unsicherheit aus R (Std. error) durch den Wert Residual standard error.
Damit erhält man aus R eine Unsicherheit des angepassten Parameters X1:
1.643E-07 / 0.0004624 = 3.5532E-04
Fall: Erkennungsgrenze
Blank count rate= 4.66670009E-08 background rate= 1.88333332E-03
Input data: variance-covariance matrix: (rank= 18 )
2.29269E-07 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-082.61574E-08 1.47460E-07 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-082.61574E-08 2.61574E-08 1.14249E-07 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-082.61574E-08 2.61574E-08 2.61574E-08 1.00767E-07 2.61574E-08 2.61574E-08 2.61574E-082.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 9.52931E-08 2.61574E-08 2.61574E-082.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 9.30711E-08 2.61574E-082.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 9.21690E-082.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-082.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-082.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-082.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-082.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-082.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-082.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-082.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-082.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-082.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-082.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08 2.61574E-08Arrays y, X1, x2, X3:
y X1 X3
1 3.96624E-03 8.80220E-01 2.73093E-012 1.61015E-03 8.07210E-01 1.10866E-013 6.53662E-04 7.40257E-01 4.50075E-024 2.65363E-04 6.78857E-01 1.82714E-025 1.07728E-04 6.22550E-01 7.41751E-036 4.37335E-05 5.70913E-01 3.01124E-037 1.77542E-05 5.23559E-01 1.22245E-038 7.20756E-06 4.80133E-01 4.96272E-049 2.94955E-06 4.40700E-01 2.03089E-0410 1.18786E-06 4.03787E-01 8.17887E-0511 4.82232E-07 3.70295E-01 3.32032E-0512 1.95772E-07 3.39582E-01 1.34793E-0513 7.94798E-08 3.11415E-01 5.47210E-0614 3.22691E-08 2.85585E-01 2.22147E-0615 1.31030E-08 2.61897E-01 9.01838E-0716 5.32199E-09 2.40175E-01 3.66113E-0717 2.16298E-09 2.20253E-01 1.48629E-0718 8.80329E-10 2.01985E-01 6.03378E-08Parameter values and std uncertatinties obtained by UR:
1 1.98944E-11 3.10543E-043 1.45234E-02 1.73864E-03Chisqr= 8.88178420E-16
For the case of the decision threshold one obtains with R
(in diesem Fall ist nur die Unsicherheit zu X1 von Interesse; der Wert von X1 sollte nahe null liegen):
> covmat <- read.table("covmat2.txt")> data <- read.table("data2.txt")> res <- lm.gls(formula = y ~ X1 + X3 - 1, data = data, W = covmat, inverse = TRUE)> summary.lm(res)Call:
lm.gls(formula = y ~ X1 + X3 - 1, data = data, W = covmat, inverse =TRUE)Residuals:Min 1Q Median 3Q Max-1.370e-09 -5.200e-13 -4.300e-13 3.721e-11 8.881e-10Coefficients:Estimate Std. Error t value Pr(>|t\|)X1 2.193e-11 1.334e-13 1.644e+02 <2e-16 \**\*X3 1.452e-02 7.467e-13 1.945e+10 <2e-16 \**\*---Signif. codes: 0 ‘\**\*’ 0.001 ‘\*\*’ 0.01 ‘\*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 4.295e-10 on 16 degrees of freedomMultiple R-squared: 1, Adjusted R-squared: 1F-statistic: 5.105e+13 on 2 and 16 DF, p-value: < 2.2e-16Warning message:In summary.lm(res) : calling summary.lm(<fake-lm-object>) ...Damit erhält man aus R eine Unsicherheit des angepassten Parameters X1:
1.334E-13 / 4.295E-10 = 3.10594E-04
5.4.4. Hinweis zur Berechnung von Erkennungs- und Nachweisgrenze:¶
Diese speziellen Ergebnisgrößen beziehen sich auf die Komponente \(y_{1}\) des Ergebnisvektors \(\mathbf{y}\) und deren Unsicherheit \(u\left( y_{1} \right)\), die iteriert wird, bis \(y_{1}\) und \(u\left( y_{1} \right)\) die Abbruchbedingung der Iteration erfüllen. \(y_{1}\) wird durch einen neuen, von der Iteration bestimmten Wert \(y_{1}^{'}\) ersetzt. Daraus erhält man gemäß dem Modell \(\mathbf{x}^{\mathbf{'}}\mathbf{= A\ }\mathbf{y}^{\mathbf{'}}\), woraus eine geänderte Kovarianzmatrix \(\mathbf{U}_{\mathbf{x}}^{\mathbf{'}}\) resultiert. Dazu werden zunächst den Bruttozählraten der Zerfallskurve neue Werte zugewiesen:
\(R_{b,i}^{'} = \left( R_{0,i} + R_{bl} \right) + y_{1}^{'} \bullet X_{1}\left( t_{i} \right) + y_{2} \bullet X_{2}\left( t_{i} \right) + y_{3} \bullet X_{3}\left( t_{i} \right)\)
Daraus folgt die Unsicherheits-Funktion (Standardunsicherheit) der Bruttozählrate:
\(u\left( R_{b,i}^{'} \right) = \sqrt{R_{b,i}^{'}/t_{m,i}}\) .
Die Nettozählraten der modifizierten Abklingkurve lauten:
\(R_{n,i}^{'} = R_{b,i}^{'} - R_{0,i} - R_{bl}\) ,
woraus sich die Diagonalelemente der geänderten Kovarianzmatrix \(\mathbf{U}_{\mathbf{x}}^{\mathbf{'}}\) (Varianzen der Nettozählraten) ergeben:
\(var\left( R_{n,i}^{'} \right) = \frac{R_{b,i}^{'}}{t_{m,i}} + var\left( R_{0,i} \right) + var\left( R_{bl} \right)\) ,
während die nicht-diagonalen Elemente unverändert bleiben.
Aus der rechten Formel von Gl. (5.4.1) wird die Unsicherheit \(u\left( y_{1}^{'} \right)\) bestimmt. Mit dem Wertepaar \(y_{1}^{'}\) und \(u\left( y_{1}^{'} \right)\) kann der nächste Iterationsschritt gestartet werden; diese Iteration kann wiederholt werden, falls das Konvergenzkriterium noch nicht erfüllt ist.
5.4.5. PMLE-Literatur¶
ISO 11929:2010, Appendix C.5.4Klaus Weise u. Wolfgang Wöger, 1999: Meßunsicherheit und Meßdatenauswertung.Verlag Wiley-VCH Weinheim,S. 200 oben (Section 5.4.2 Lineare Kurvenanpassung)Roger J. Barlow, 1999: Statistics. A Guide to the Use of Statistical Methods in the PhysicalSciences. The Manchester Physics Series. John Wiley & Sons Ltd., Chichester, New York.Section 6.6, pp. 111-113.
Zur Durchführung der Berechnung werden Matrix-Routinen aus der Datan-Library verwendet (in Fortran 90 konvertiert):
Datan-Library aus:Siegmund Brandt: Datenanalyse. Mit statistischen Methoden undComputerprogrammen; 4. Auflage. Spektrum, Akademischer Verlag,Heidelberg-Berlin, 1999. In German.Dieses Lehrbuch ist auch in englischer Sprache verfügbar.
Weitere:
T. Hauschild, M. Jentschel, 2001: Comparison of maximum likelihood estimation and chi-square statistics applied to counting experiments. Nucl. Instr. & Meth A 457 (1-2), S 384-401.
Pommé, S., Keightley, J., 2007. Countrate estimation of a Poisson process: unbiased fit versus central moment analysis of time interval spectra. Applied Modeling and Computations in Nuclear Science. In: Semkow, T.M., Pommé, S., Jerome, S.M., Strom, D.J. (Eds.), ACS Symposium Series 945. American Chemical Society, Washington, DC, pp.316–334.2007.ISBN0-8412-3982-7.
Laurence, T. A., Chromy, B., 2009. Efficient Levenberg-Marquardt Minimization of the Maximum Likelihood Estimator for Poisson Deviates. Report LLNL-JRNL-420247, November 13, 2009.
