jackknife und bootstrap - users.ph.tum.deusers.ph.tum.de/ga32pex/knap_bachelor.pdf · algorithmus...
Post on 07-Feb-2018
218 Views
Preview:
TRANSCRIPT
Jackknife und Bootstrap
Michael Knap
B A K K A L A U R E ATS AR BE I T
eingereicht am
Bachelorstudiengang
Technische Physik
am
Institut fur Theoretische Physik –Computational Physics
an der
im Juni 2007
ii
Diese Arbeit ist entstanden unter Betreuung von:
Ao.Univ.-Prof. Dr. Hans Gerd Evertz
iii
iv
Inhaltsverzeichnis
Kurzfassung vii
Abstract ix
1 Einleitung 1
2 Traditioneller Zugang 3
3 Data-Resampling 93.1 Jackknife . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103.2 Bootstrap . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.2.1 Der Bootstrap-Algorithmus . . . . . . . . . . . . . . . 123.2.2 Konfidenzintervalle . . . . . . . . . . . . . . . . . . . . 15
3.3 Theoretische Hintergrunde . . . . . . . . . . . . . . . . . . . . 163.4 Fehlerfortpflanzung . . . . . . . . . . . . . . . . . . . . . . . . 163.5 Jackknife vs. Bootstrap . . . . . . . . . . . . . . . . . . . . . 17
4 Anwendungen 194.1 Lineare Regression . . . . . . . . . . . . . . . . . . . . . . . . 19
4.1.1 Mikroskop . . . . . . . . . . . . . . . . . . . . . . . . . 194.2 Regression mit nichtlinearen Parametern . . . . . . . . . . . . 22
4.2.1 Modellfunktion y = axb e−c x . . . . . . . . . . . . . . 224.2.2 Fernrohr . . . . . . . . . . . . . . . . . . . . . . . . . . 27
4.3 Integrierte Autokorrelationszeit . . . . . . . . . . . . . . . . . 31
5 Zusammenfassung 39
A Matlab-Quellcode 41A.1 Stichprobenmittelwert – Bootstrap . . . . . . . . . . . . . . . 41A.2 Lineare Regression . . . . . . . . . . . . . . . . . . . . . . . . 43A.3 Regression mit nichtlinearen Parametern . . . . . . . . . . . . 45A.4 Integrierte Autokorrelationszeit . . . . . . . . . . . . . . . . . 54
Literaturverzeichnis 59
v
vi
Kurzfassung
Die statistische Analyse von Daten ist zweifellos ein wichtiger Bereich in derPhysik und in anderen Wissenschaften. Der Grund dafur liegt im umfangrei-chen Anwendungsgebiet. Statistische Analyse beinhaltet immer den Aspektder Unsicherheit von Schatzwerten. Daher ist es ein Ziel diese Unsicherhei-ten zu quantifizieren. Das kann durch analytische Berechnung, basierendauf einem Wahrscheinlichkeitsmodell der Daten, erfolgen. Der Nachteil die-ser Methode ist, dass sie sowohl langwierig als auch schwierig durchzufuhrenist. Ein ganz anderer Zugang zur Berechnung der Qualitat von Schatzwertenerfolgt uber die Data-Resampling Methoden Jackknife und Bootstrap.
Das Ziel dieser Arbeit ist es, eine Einfuhrung in die Theorie von Jackknifeund Bootstrap zu liefern. Um die Methoden besser verstehen zu konnen,wurde die theoretische Beschreibung durch Beispiele in den Bereichen derlinearen und nichtlinearen Regression erganzt.
Als Ergebnis kann festgehalten werden, dass Jackknife und Bootstrapdie statistische Analyse in vielen Fallen vereinfachen. Nichtsdestotrotz istdie Anwendung von Data-Resampling Methoden nur mit einer theoretischenBestatigung gerechtfertigt.
Wegen der einfachen Anwendung von Data-Resampling und der hervor-ragenden Leistung bei anspruchsvollen Problemen, sind Jackknife und Boot-strap zur Bestimmung der Qualitat von Schatzwerten sehr zu empfehlen.
vii
viii
Abstract
It goes without saying that statistical analysis of data is very importantin physics and other sciences due to its wide range of application. Statis-tical procedures always imply the aspect of uncertainty of estimators andtherefore the quality of estimators must be evaluated. This can be done byanalytical calculation based on an assumed probability model for the data.The drawback of this method is that it is both time consuming and difficult.A completely different approach to calculating the quality of estimators isbased on data-resampling methods like Jackknife and Bootstrap.
The aim of this thesis is to provide an introduction to the theory of Jack-knife and Bootstrap. For clarification, these methods have been illustratedwith some examples of linear and nonlinear regression.
As a result, it should be mentioned that Jackknife and Bootstrap simplifythe statistical analysis in many cases. Yet the application of data-resamplingmethods is only justified through theoretical confirmation.
Given their straight-forward usage and their outstanding performancein sophisticated problems, Jackknife and Bootstrap can be highly recom-mended for calculating the quality of estimators in statistical theory.
ix
x
Kapitel 1
Einleitung
Ziel der Statistik ist es, moglichst viele Informationen aus gegebenen Datenherauszufiltern. Dazu muss eine Verteilung der Daten bestimmt werden, ausder schließlich weitere Eigenschaften abgeleitet werden konnen. Eine solcheEigenschaft kann beispielsweise der Stichprobenmittelwert, die Stichproben-varianz oder ein Konfidenzintervall sein. Aus Daten mit endlichem Umfangkann nur ein Schatzwert einer statistischen Große berechnet werden. Da-her ist es fur weitere Aussagen notwendig die Qualitat des Schatzwerts zubestimmen.
Die Qualitat eines Schatzwerts ist durch die Varianz und den Bias ge-geben. Eine analytische Berechnung der Varianz und des Bias eines Schatz-werts aus einer bestimmten Verteilung der Daten ist aber meistens mit sehrgroßem Aufwand verbunden.
Diese Arbeit liefert eine Einfuhrung in die Data-Resampling MethodenJackknife und Bootstrap, die aus den ursprunglichen Daten wiederholt Da-tensatze generieren. Ohne komplizierte und langwierige analytische Rech-nung kann aus diesen Datensatzen direkt auf die Varianz und den Bias ei-nes Schatzwerts geschlossen werden. Da dieser Ansatz eine große Menge anDaten produziert und eine wiederholte Auswertung des ursprunglichen Pro-blems erfordert, wird die Berechnung nach dem Jackknife- und Bootstrap-Algorithmus in den meisten Fallen numerisch durchgefuhrt.
Das erste Kapitel liefert einen Uberblick uber den traditionellen Zugangzur Berechnung von Schatzwerten und deren Varianz. Im zweiten Kapitelwerden die Algorithmen der Data-Resampling Methoden definiert und miteinem kurzen Beispiel erklart. Daruber hinaus wird noch speziell auf dieVorteile von Bootstrap bei der Berechnung von Konfidenzintervallen ein-gegangen. Abschließend enthalt dieses Kapitel eine Gegenuberstellung derbeiden Methoden. Im dritten Kapitel werden Anwendungen von Jackknifeund Bootstrap behandelt, welche den Bereich der linearen und nichtlinearenRegression umfassen.
1
2
Kapitel 2
Traditioneller Zugang
Es wird der Unterschied zwischen parametrischen und nicht-parametrischenSchatzwerten definiert und die empirische Verteilungsfunktion (EDF – Em-pirical Distribution Function) eingefuhrt. Auf dieser Basis erfolgt eine ana-lytische Berechnung der Varianz des Stichprobenmittelwerts. Danach wirdder traditionelle Zugang zur Fehlerfortpflanzung erlautert. Am Ende diesesKapitels werden die Nachteile des traditionellen Zugangs zur Berechnungvon Schatzwerten aufgezeigt.
Ein Ziel der statistischen Analyse ist es, mit Hilfe vorhandener Daten eineunbekannte Große θ zu approximieren1. Die Daten (x1, x2, . . . xn) werdenaus Zufallsvariablen (X1, X2, . . . Xn) erzeugt, die der Wahrscheinlichkeits-dichte f(x) (PDF – Probability Density Function) und der Verteilungsfunk-tion F (x) (CDF – Cumulative Distribution Function) genugen. Im allge-meinen Fall besitzt jede Zufallsvariable Xi eine eigene VerteilungsfunktionFi(x).
θ hangt von den Verteilungsfunktionen der Zufallsvariablen ab und kannbeispielsweise der Mittelwert, der Median, die Varianz oder der Achsenab-schnitt bei der linearen Regression sein. Die unbekannte Große θ kann nurexakt berechnet werden, falls die Verteilungsfunktionen Fi(x) bekannt sind.In der Praxis sind jedoch die Daten (x1, x2, . . . xn) vorgegeben. Die Vertei-lungsfunktionen – und daher auch θ – sind wegen dem endlichen Umfangder Daten nicht exakt festgelegt. Es ist jedoch moglich, einen Schatzwert θfur die unbekannte Große θ zu bestimmen.
Schatzwerte konnen berechnet werden, indem Annahmen uber die Ver-teilung der Daten getroffen werden. Das bedeutet, dass ein von Parame-tern abhangiges Modell eingefuhrt wird, das bestmoglich zu den Daten(x1, x2, . . . xn) passt. Die Schatzwerte sind in diesem Zugang unbekannteVariablen des Modells und werden als parametrische Schatzwerte bezeichnet.Ein Beispiel hierzu ist der Maximum-Likelihood-Schatzwert.
1Die Grundlagen dieser Arbeit basieren auf [2,7,10].
3
Kapitel 2. Traditioneller Zugang
Eine weitere Moglichkeit besteht darin, Schatzwerte direkt aus den Daten(x1, x2, . . . xn) zu berechnen. Diese Schatzwerte heißen nicht-parametrischeSchatzwerte, da fur sie kein von Parametern abhangiges Wahrscheinlich-keitsmodell hergeleitet werden muss. Stattdessen wird die empirische Ver-teilungsfunktion (EDF – Empirical Distribution Function) verwendet, diedurch
F (x) =1
n
n∑
i=1
H(x− xi) (2.1)
definiert ist. Die Funktion H(x) in (2.1) ist die Heavisidesche Stufenfunk-tion. Der Schatzwert θ wird im nicht-parametrischen Fall uber die empirischeVerteilungsfunktion F (x) bestimmt.
Beispiel 2.1 (Stichprobenmittelwert) Der Erwartungswert des Mittel-werts der Verteilung ist durch
〈x〉 =
∫
x dF (x)
definiert. Wird die Verteilungsfunktion F (x) durch die empirische Vertei-lungsfunktion F (x) ersetzt, entsteht ein Schatzwert θ fur den Mittelwert.Das Ergebnis fur den Schatzwert θ ist der Stichprobenmittelwert x.
θ =
∫
x dF (x)(2.1)=
∫
x d( 1
n
n∑
i=1
H(x− xi))
=1
n
n∑
i=1
∫
x dH(x− xi)
=1
n
n∑
i=1
∫
x δ(x− xi) dx =1
n
n∑
i=1
xi = x (2.2)
Dieses Beispiel wurde aus [2, S. 13] entnommen.
Zur Abschatzung der Qualitat des Schatzwerts ist zumindest eine An-gabe der Varianz von θ notwendig.
Beispiel 2.2 (Varianz des Stichprobenmittelwerts) Aus dem zentra-len Grenzwertsatz folgt fur die Varianz des Stichprobenmittelwerts θ die Be-ziehung
var(θ) =var(x)
n.
var(x) ist die Varianz der Stichprobe (Details zum zentralen Grenzwertsatzin [6, Abschn. 8.6]).
Die Stichprobenvarianz kann durch
var(x) =
∫
(
x−
∫
y dF (y))2dF (x)
4
Kapitel 2. Traditioneller Zugang
berechnet werden. Wird in dieser Gleichung erneut die VerteilungsfunktionF (x) durch die empirische Verteilungsfunktion F (x) ersetzt, ergibt sich derSchatzwert der Stichprobenvarianz.
ˆvar(x) =
∫
(
x−
∫
y dF (y))2dF (x)
(2.2)=
∫
(x− θ)2 dF (x)
(2.1)=
1
n
n∑
i=1
∫
(x− θ)2 dH(x− xi)
=1
n
n∑
i=1
∫
(x− θ)2 δ(x − xi) dx =1
n
n∑
i=1
(xi − θ)2
Mit diesem Ergebnis ist die Varianz des Stichprobenmittelwerts durch denSchatzwert
ˆvar(θ) =1
n2
n∑
i=1
(xi − θ)2
festgelegt.
In der Praxis muss haufig die Unsicherheit einer Große x bestimmt wer-den, wobei die Große x eine Funktion von mindestens zwei Variablen u, v, . . .ist.
x = f(u, v, . . .) (2.3)
u, v, . . . sind nicht eindeutig bestimmt, sondern liegen als fehlerbehafteteSchatzwerte u, v, . . . vor. Die Unsicherheiten von u, v, . . . ziehen eine Unsi-cherheit von x nach sich. Deshalb wird von Fehlerfortpflanzung gesprochen.
Da u, v, . . . Schatzwerte sind, kann auch von x nur ein Schatzwert
x = f(u, v, . . .) (2.4)
berechnet werden. Dieser Schatzwert unterliegt der Varianz var(x) = σ2x,
die nun bestimmt werden soll. Durch individuelle Datenpunkte ui, vi, . . .,die von den zugehorigen Verteilungen generiert werden, konnen individuelleErgebnisse fur xi berechnet werden.
xi = f(ui, vi, . . .) (2.5)
Damit ergibt sich fur die Varianz des Schatzwerts x
σ2x = lim
n→∞
1
n
n∑
i=1
(xi − x)2 . (2.6)
Die erste Ordnung der Taylorreihenentwicklung von xi um x liefert
xi ≈ x+ ∆u
(
∂f
∂u
)
+ ∆v
(
∂f
∂v
)
+ . . .
xi − x ≈ (ui − u)
(
∂f
∂u
)
+ (vi − v)
(
∂f
∂v
)
+ . . . . (2.7)
5
Kapitel 2. Traditioneller Zugang
Der Abbruch der Taylorreihenentwicklung nach dem ersten Glied ist nurgerechtfertigt, falls die Unsicherheiten ∆u, ∆v, . . . klein sind. Ansonstenmussen zusatzlich hohere Ordnungen der Reihenentwicklung berucksichtigtwerden.
Wird (2.7) in (2.6) eingesetzt, folgt
σ2x ≈ lim
n→∞
1
n
n∑
i=1
[
(ui − u)
(
∂f
∂u
)
+ (vi − v)
(
∂f
∂v
)
+ . . .
]2
≈ limn→∞
1
n
n∑
i=1
[
(ui − u)2(
∂f
∂u
)2
+ (vi − v)2(
∂f
∂v
)2
+ . . .
+ 2(ui − u)(vi − v)
(
∂f
∂u
)(
∂f
∂v
)
+ . . .
]2
. (2.8)
Mit der Definition der Varianz in (2.6) ist es moglich, die ersten bei-den Terme von (2.8) zu vereinfachen. Der dritte Ausdruck entspricht derKovarianz zwischen den Parametern u und v, die durch
σ2uv ≡ lim
n→∞
1
n
n∑
i=1
(ui − u)(vi − v) (2.9)
definiert ist. Damit ergibt sich fur die Varianz von x
σ2x ≈ σ2
u
(
∂f
∂u
)2
+ σ2v
(
∂f
∂v
)2
+ . . .+ 2σ2uv
(
∂f
∂u
)(
∂f
∂v
)
+ . . . . (2.10)
Sind die Parameter u, v, . . . unkorreliert, entfallt der Term der Kovari-anzen und es entsteht die vereinfachte Fehlerfortplanzungsformel
σ2x ≈ σ2
u
(
∂f
∂u
)2
+ σ2v
(
∂f
∂v
)2
+ . . . . (2.11)
Gleichung (2.10) kann auch in kompakter Matrixnotation angegeben wer-den.
σ2x ≈ JTCJ (2.12)
J =
∂f∂u∂f∂v...
C =
σ2u σ2
uv · · ·σ2
uv σ2v
.... . .
J entspricht der Jacobi-Matrix und C der symmetrischen Kovarianzmatrix.Die Herleitung der Fehlerfortflanzungsformel wurde aus [1, Kap. 3] entnom-men.
Die Nachteile des traditionellen Zugangs zur Berechnung von Schatzwer-ten und der Fehlerfortpflanzung sind:
6
Kapitel 2. Traditioneller Zugang
1. Die Berechnung von komplizierten Schatzwerten θ ist schwierig durch-zufuhren.
2. Nicht alle Schatzwerte sind exakt, wenn sie auf die oben angefuhrteWeise berechnet werden. Beispielsweise muss die Varianz einer Stich-probe auf n
n−1 ˆvar(x) korrigiert werden, um einen systematischen Feh-ler (Bias) zu unterbinden.
3. Ist beim parametrischen Zugang die Wahl des Modells nicht korrekt, sokonnen bei der Berechnung von weiterfuhrenden Schatzwerten großeFehler auftreten.
4. Die analytische Berechnung des Schatzwerts muss beim parametri-schen Fall fur jedes Problem erneut erfolgen.
5. Verlauft die Fehlerfortpflanzung eines Schatzwerts uber viele Schritte,so erfordert das die wiederholte Berechnung der Kovarianzmatrix undder Jacobi-Matrix des Schatzwerts. Bei einer großen Anzahl an Schrit-ten fuhrt das zu einem enormen Aufwand. Ein Beispiel fur wiederholteFehlerfortpflanzung ist
α = f(θ, a)
β = g(α, b). . .
ω = z(ψ, x) .
Es wird der ursprungliche Schatzwert θ uber viele Zwischenschritte aufω fortgepflanzt. Dabei sind a, b, . . . x zusatzliche Parameter, die indie Funktionen f , g, . . . z einfließen. Der traditionelle Zugang erfor-dert eine wiederholte Berechnung der Kovarianzmatrix und der Jacobi-Matrix der Schatzwerte θ bis ψ, um den Fehler von ω ermitteln zukonnen.
Diese Nachteile treten bei Verwendung von Data-Resampling-Methodennicht auf.
7
8
Kapitel 3
Data-Resampling
Angenommen, es soll eine Große θ durch einen Schatzwert θ approximiertwerden, so ist von vornherein nicht klar, wie exakt dieser Schatzwert ist. EineMoglichkeit zur Beurteilung der Qualitat ist die Bestimmung des mittlerenquadratischen Fehlers (MSE - Mean Square Error)
MSE = var(θ) +B2(θ) . (3.1)
Dabei ist var(θ) in (3.1) die Varianz von θ, die durch
var(θ) = 〈(θ − 〈θ〉)2〉 = 〈θ2〉 − 〈θ〉2
definiert ist und B(θ) der Bias (Systematischer Fehler), der durch
B(θ) = 〈θ〉 − θ
gegeben ist. 〈. . .〉 entspricht dem Erwartungswert einer Zufallsvariable. DerErwartungswert einer Zufallsvariable X ist durch
〈X〉 =
∫
x∈S
x f(x) dx (3.2)
definiert, wobei uber alle x aus dem Wahrscheinlichkeitsraum S integriertwird und f(x) die Wahrscheinlichkeitsdichte ist.
Eine Moglichkeit zur Berechnung der Qualitat des Schatzwerts ist durchVerwendung von Data-Resampling gegeben. Data-Resampling besagt, dassdie vorhandenen Daten mit Hilfe bestimmter Algorithmen neu angeordnetund strukturiert werden. Durch die neue Datenstruktur konnen weitere In-formationen uber den Schatzwert, wie z. B. dessen Varianz oder MSE, gewon-nen werden. Der Vorteil dieser Methode liegt darin, dass diese Informationenaus einer einzigen Stichprobe stammen. Ansonsten mussten viele Stichpro-ben erzeugt werden, was aber nicht sinnvoll und oft auch nicht moglich ist.
9
Kapitel 3. Data-Resampling
3.1 Jackknife
Es wird der Jackknife-Algorithmus, der wiederholt Datensatze aus den ur-sprunglichen Daten generiert, definiert. Die Formeln fur Varianz und Biasdes Schatzwerts werden auf Basis der neu generierten Datensatze festgelegt.Daruber hinaus wird in einem Beispiel die Varianz des Stichprobenmittel-werts mit Hilfe des Jackknife-Algorithmus berechnet.
Eine Methode zur Berechnung der Qualitat von Schatzwerten auf Basisvon Data-Resampling ist der Jackknife-Algorithmus. Dabei werden die Zu-fallsvariablen (X1, X2, . . . Xn) zur Bestimmung von θ herangezogen. Die nZufallsvariablen stammen aus der zugrunde liegenden Wahrscheinlichkeits-dichte.
Der delete-1-Jackknife-Algorithmus (vgl. [10, S. 152–153]):
1. Es werden n Großen θ(i) berechnet. Die Berechnung von θ(1) erfolgt
analog zur Berechnung von θ, jedoch wird der 1. Datenpunkt weggelas-sen, d. h. es tragen die n−1 Werte (X2, X3, . . . Xn) zu θ(1) bei. Bei der
Berechnung von θ(2) wird der 2. Datenpunkt ausgelassen und es werdendie n−1 Werte (X1, X3, . . . Xn) verwendet. Allgemein werden zur Be-rechnung von θ(i) die Daten (X1, . . . Xi−1, Xi+1, . . . Xn) verwendet.
Durch diese Vorschrift entstehen n neue Werte θ(1), θ(2), . . . θ(n). DieseWerte, die aus n − 1 Datenpunkten entstanden sind, ahneln sich undunterscheiden sich auch nicht stark von θ.
2. Es sei
θ(·) =1
n
n∑
i=1
θ(i) (3.3)
der Mittelwert der neu erzeugten Werte.
3. Wird nun θ(·) mit θ verglichen, ergibt sich ein Schatzwert des Biasdurch
B = (n− 1)(θ(·) − θ) . (3.4)
4. Der Bias gibt die Differenz zwischen dem Schatzwert und der korrektenGroße an. Damit ist es moglich, einen korrigierten Wert θ zu berechnen,der naher beim tatsachlichen Wert θ liegt als θ.
θ = θ − B (3.5)
5. Ein Schatzwert der Varianz ist durch
ˆvar =n− 1
n
n∑
i=1
(θ(i) − θ(·))2 (3.6)
10
Kapitel 3. Data-Resampling
gegeben. Die Beziehung (3.6) entspricht bis auf den Faktor (n−1)2 derVarianz des Stichprobenmittelwerts. Dieser zusatzliche Faktor entstehtdurch die Abhangigkeit der θ(i). Die θ(i) sind voneinander abhangig,da sie aus einer einzigen Stichprobe durch Data-Resampling erzeugtwurden.
Diese Behandlung wird als delete-1-Jackknife bezeichnet, da die Berechnungvon θ(i) durch Weglassen eines Datenpunkts erfolgt. Die verallgemeinerteForm des Jackknife-Algorithmus ist der delete-d-Jackknife, bei dem bei derBerechnung d Elemente weggelassen werden. Bei großen Datenmengen ver-ringert sich dadurch der Rechenaufwand. In der Praxis wird es oft so ge-handhabt, dass die zugrunde liegende Datenmenge in ungefahr NBlock = 30Blocke eingeteilt wird und bei jedem Jackknife-Schritt ein ganzer Block weg-gelassen wird. In den Gleichungen (3.3), (3.4) und (3.6) muss n durch NBlock
ersetzt werden, da nur noch NBlock Großen θ(i) berechnet werden (weitereDetails in [7, Abschn. 2.3 und Abschn. 5.2]).
Beispiel 3.1 (Stichprobenmittelwert) Der Schatzwert fur den Stichpro-benmittelwert ist
θ =1
n
n∑
i=1
Xi = X .
Die Anwendung des delete-1-Jackknife auf den Stichprobenmittelwert erfor-dert folgende Rechenschritte:
θ(i) =1
n− 1
∑
j 6=i
Xj =1
n− 1(n X −Xi), ∀i = 1, 2 . . . n
θ(·) =1
n
n∑
i=1
θ(i) =1
n (n− 1)
n∑
i=1
(n X −Xi)
=1
n− 1
n∑
i=1
X −1
n− 1
1
n
n∑
i=1
Xi
=n
n− 1X −
1
n− 1X =
1
n− 1(n X − X)
= X
B = (n− 1)(θ(·) − θ) = 0
θ = θ − B = X
ˆvar =n− 1
n
n∑
i=1
(θ(i) − θ(·))2
11
Kapitel 3. Data-Resampling
=n− 1
n
n∑
i=1
(n X −Xi
n− 1− X
)2
=n− 1
n
n∑
i=1
(X −Xi
n− 1
)2=
1
n (n− 1)
n∑
i=1
(X −Xi)2
Die Anwendung des Jackknife-Algorithmus liefert in Beispiel 3.1 keineneuen Informationen. Es ist bereits bekannt, dass der Stichprobenmittelwertkeinen Bias besitzt und dass die Varianz des Stichprobenmittelwerts vonder Form 1
n (n−1)
∑ni=1(X −Xi)
2 ist. Jedoch ist es interessant zu sehen, dassdurch Anwendung von Jackknife die bereits bekannten und guten Losungenfur den Stichprobenmittelwert reproduziert werden.
In komplizierteren Fallen ist es nicht moglich, analytische Losungen fur θ,sowie dessen Bias und Varianz zu finden. Die Berechnung muss dann nume-risch erfolgen. Dabei muss θ(i) fur alle i = 1, 2 . . . n, also n-Mal, berechnetwerden.
Der Jackknife-Algorithmus wird in der Praxis angewendet, um den Biasdes Schatzwerts abzusenken. Das erfolgt oft ohne merkbare Erhohung derVarianz [10, Abschn. 7.1]. Weiters liefert diese Anwendung eine Moglichkeitzur Abschatzung der Varianz des Schatzwerts. Zur Berechnung der Varianzist im Gegensatz zum traditionellen Ansatz auch keine komplizierte Herlei-tung des theoretischen Modells notwendig.
3.2 Bootstrap
Bootstrap generiert wiederholt Datensatze durch Ziehen mit Zurucklegenaus der ursprunglichen Datenmenge. Aus jedem dieser Datensatze kann einSchatzwert berechnet werden. Dadurch entsteht eine Verteilung des Schatz-werts. Aus dieser Verteilung wird die Varianz, der Bias und ein Konfidenz-intervall, das nicht notwendig symmetrisch sein muss, berechnet. Die Mog-lichkeit eine Verteilung des Schatzwerts mit beliebigem Umfang zu erzeugen,macht den Vorteil von Bootstrap gegenuber Jackknife aus.
3.2.1 Der Bootstrap-Algorithmus
Eine weitere Moglichkeit zur Berechnung der Qualitat von Schatzwertenbietet der Bootstrap-Algorithmus. Auch hier existiert eine Stichprobe vomUmfang n, aus der der Schatzwert θ berechnet wird.
Der Bootstrap-Algorithmus:
1. Der Bootstrap-Algorithmus erzeugt eine neue Stichprobe durch Ziehenmit Zurucklegen von n Elementen aus dem ursprunglichen Datensatz.Aus dieser neu erzeugten Stichprobe wird der Schatzwert θ∗ berechnet.
12
Kapitel 3. Data-Resampling
2. Das Erzeugen einer neuen Stichprobe wird B Mal wiederholt. Dar-aus folgen die Schatzwerte θ∗1, θ
∗2, . . . θ
∗B. Diese Schatzwerte sind eine
Approximation der Verteilung von θ.
3. Der Mittelwert der Verteilung ist durch
θ∗(·) =1
B
B∑
i=1
θ∗i (3.7)
gegeben.
4. Wird nun θ∗(·) mit θ verglichen, ergibt sich ein Schatzwert des Biasdurch
B = θ∗(·) − θ . (3.8)
5. Ist der Mittelwert der Verteilung θ∗(·) großer als θ, so wird angenommen,
dass der Schatzwert großer ist als der wahre Wert θ. Ist θ∗(·) kleiner als
θ, so gilt das Umgekehrte. Daher folgt fur den korrigierten Wert
θ∗ = θ − B . (3.9)
6. Die Varianz kann ebenfalls direkt aus der Verteilung berechnet werden.
ˆvar =1
B − 1
B∑
i=1
(θ∗i − θ∗(·))2 (3.10)
Die Wahl der Anzahl der WiederholungenB ist abhangig von der Anwen-dung des Bootstrap-Algorithmus. Soll die Varianz, die wegen der Simulationam Computer entsteht, hochstens 10% der Varianz der Daten betragen, somuss B > 10n gewahlt werden (weitere Details siehe [2, Abschn. 2.5.2]).
Beispiel 3.2 (Stichprobenmittelwert) Die Simulation durch Bootstrapkann ausschließlich uber den numerischen Zugang durchgefuhrt werden. In
Tabelle 3.1: Bootstrap-Simulation zur Berechnung des Stichprobenmittel-werts und dessen Varianz.
Direkte Methode Bootstrap-Simulation
θ −1,037 e−002 −1,042 e−002
ˆvar +0,011 e−002 +0,012 e−002
ˆstd +1,061 e−002 +1,103 e−002
B – +0,004 e−002
diesem Beispiel wird eine Stichprobe vom Umfang n = 100 erzeugt. Die
13
Kapitel 3. Data-Resampling
Stichprobe soll normalverteilte Daten um den Mittelwert θ = 0 mit einerStandardabweichung σ = 0,1 enthalten. Der Bootstrap-Algorithmus wird furB = 10n− 1 = 999 Wiederholungen durchgefuhrt.
Das Ergebnis der Simulation ist in Tabelle 3.1 dargestellt. Es ist ersicht-lich, dass Mittelwert und Varianz aus der Bootstrap-Simulation gut mit demdirekt erhaltenen Wert ubereinstimmen. Außerdem ist der Bias B viel klei-ner als die Standardabweichung ˆstd und damit nicht signifikant. Es ist nichtnotwendig den MSE (mittleren quadratischen Fehler) zu berechnen, da sichdieser kaum von der Varianz unterscheidet.
Die Verteilung von θ, aus der Mittelwert und Varianz berechnet werden,ist in Abbildung 3.1 dargestellt.
−0.05 −0.04 −0.03 −0.02 −0.01 0 0.01 0.02 0.030
5
10
15
20
25
30
35
40
θi*
f(θ i* )
Abbildung 3.1: Verteilung des Stichprobenmittelwerts.
Der Matlab-Quellcode fur dieses Beispiel ist im Anhang A.1 (Seite 41)angefuhrt.
Eine Stichprobe vom Umfang n besitzt 2n − 1 nichtleere Untermengen.Der delete-1-Jackknife-Algorithmus verwendet nur n von diesen. Der Grund-gedanke bei der Einfuhrung von Bootstrap war, eine großere Anzahl an Un-termengen oder sogar alle 2n − 1 Untermengen zur Berechnung des Schatz-werts heranzuziehen (siehe [7, S. 9]). Da es in der Praxis bei Stichprobengroßeren Umfangs nicht moglich ist alle 2n−1 Untermengen zu untersuchen,wird die Monte-Carlo-Methode zur Approximation des Schatzwerts ver-wendet. Unter Monte-Carlo-Methode ist das Ziehen einer Stichprobe vomUmfang B aus der Population mit 2n − 1 Elementen zu verstehen.
14
Kapitel 3. Data-Resampling
3.2.2 Konfidenzintervalle
Die Berechnung von Konfidenzintervallen ist eine sehr wichtige Anwendungdes Bootstrap-Algorithmus. Dabei muss nicht davon ausgegangen werden,dass der Schatzwert normalverteilt ist, da der Bootstrap-Algorithmus selbereine geschatzte Verteilung von θ liefert.
θ↑i seien die Werte θ∗i nach aufsteigender Reihenfolge sortiert. Das (1−α)Konfidenzintervall ist durch
[θ↑((B+1)α2), θ
↑((B+1)(1−α
2))] (3.11)
definiert. Damit ist es moglich, die Konfidenzintervalle von verschieden ver-teilten Schatzwerten zu berechnen.
Beispiel 3.3 (Konfidenzintervall) In Abbildung 3.2 sind die Konfidenz-intervalle des Stichprobenmittelwerts normalverteilter Daten in Abhangigkeitvon der Anzahl von Berechnungen B aufgetragen. Der Matlab-Quellcode fur
0 200 400 600 800 1000−0.01
0
0.01
0.02
0.03
0.04
0.05
B
Kon
fiden
zint
erva
ll vo
n θ (⋅)*
Bootstraptheoretisch
Abbildung 3.2: 90% Konfidenzintervall des Mittelwerts normalverteilterDaten.
dieses Beispiel ist im Anhang A.1 (Seite 42) angefuhrt.
Die Genauigkeit der Konfidenzintervalle ist sehr stark abhangig von derAnzahl B der durchgefuhrten Berechnungen (siehe Beispiel 3.3). B mussgroß genug sein, damit die Anwendung von (3.11) sinnvoll ist (weitere Detailssiehe [2, Abschn. 2.5.2]).
15
Kapitel 3. Data-Resampling
3.3 Theoretische Hintergrunde
Eine Anwendung von Data-Resampling darf nur erfolgen, falls ein konsis-tenter Schatzwert mit hinreichender asymptotischer Genauigkeit vorliegt.
Jackknife und Bootstrap sind intuitiv anzuwenden. Trotzdem ist bei ihrerNutzung eine theoretische Bestatigung wichtig. Dies ist vor allem bei kom-plizierten Modellen notwendig. Dabei wird in der Statistik oft vom asym-ptotischen Verhalten einer Stichprobe vom Umfang n mit limn→∞ ausgegan-gen. Asymptotisch korrekte Schatzwerte sind naturlich immer asymptotischunkorrekten vorzuziehen. Bei der theoretischen Betrachtung von Modellensollen folgende Punkte behandelt werden:
1. Konsistenz
2. Asymptotische Genauigkeit
Konsistenz bedeutet, dass im asymptotischen Verhalten der Schatzwertgegen den wahren Wert strebt. Die asymptotische Genauigkeit gibt an, mitwelcher Ordnung der Fehler verschwindet (z. B. 1
n, 1
n12
, 1
n14
, . . .). Eine aus-
fuhrliche theoretische Behandlung von Jackknife und Bootstrap ist in [7] zufinden.
3.4 Fehlerfortpflanzung
Es wird die Durchfuhrung der Fehlerfortpflanzung mit Data-Resampling er-klart und aufgezeigt welche Vorteile dadurch entstehen.
Der traditionelle Zugang der Fehlerfortpflanzung durch Berechnung derKovarianzmatrix und der Jacobi-Matrix wurde bereits im Kap. 2 erlautert.Nach dieser Methode wird die Fehlerfortpflanzung aufwendig, falls sie uberviele Zwischenschritte wiederholt erfolgt. Das angegebene Beispiel fur Feh-lerfortpflanzung uber mehrere Schritte war
α = f(θ, a)
β = g(α, b). . .
ω = z(ψ, x) .
Damit die Unsicherheit von ω angegeben werden kann, mussen in jedemSchritt die Ableitungen nach den abhangigen Variablen berechnet werden.Zusatzlich muss jedes Mal die Kovarianzmatrix berechnet werden.
Ganz anders ist der Zugang uber Data-Resampling. Dabei wird die ge-wunschte Große nicht nur auf Basis des ursprunglichen Datensatzes berech-
16
Kapitel 3. Data-Resampling
net, sondern es werden alle kunstlich generierten Datensatze verwendet. Da-durch entstehen die neuen Werte
θ(1), θ(2), . . . θ(n)
α(1), α(2), . . . α(n)
β(1), β(2), . . . β(n)
. . .
ω(1), ω(2), . . . ω(n) ,
wobei α(1) = f(θ(1), a), α(2) = f(θ(2), a), . . . sind. Diese Werte sind in der
Jackknife-Notation angegeben. Fur Bootstrap gilt naturlich das Aquivalente.Auf den Datenvorrat der jeweiligen Variable konnen die Formeln fur Varianzund Bias des verwendeten Algorithmus direkt angewendet werden. Dadurchist eine Fehlerfortpflanzung ohne langwierige Berechnung von Ableitungenund Rucksichtnahme auf Kovarianzmatrizen moglich.
3.5 Jackknife vs. Bootstrap
Die Vorteile von Data-Resampling gegenuber dem traditionellen Zugang wer-den aufgezeigt. Außerdem erfolgt eine Gegenuberstellung von Jackknife undBootstrap.
1. Jackknife und Bootstrap benotigen zur Fehlerfortpflanzung im Ge-gensatz zur traditionellen Methode keine Kovarianzmatrix und keineJacobi-Matrix.
2. Beide haben gegenuber dem traditionellen Zugang den Vorteil, dasssie keine Herleitungen des theoretischen Modells benotigen.
3. Der nicht-parametrische Bootstrap und Jackknife basieren nicht auftheoretischen Modellen und daher mussen keine Annahmen bezug-lich des theoretischen Modells getroffen werden. Dies wird als Robust-heit bezeichnet. Beispielsweise kann gezeigt werden, dass der JackknifeVarianz-Schatzwert in linearen Regressionsmodellen robust gegenuberunterschiedlichen Varianzen der Datenpunkte (Heteroskedastizitat) ist(siehe [4] und [9]).
4. Der Jackknife Varianz-Schatzwert (3.6) ist fur viele Statistiken konsis-tent, jedoch benotigt er eine strengere Glattheit von θ als der Boot-strap Varianz-Schatzwert (3.10). Beispielsweise ist der Schatzwert derJackknife Varianz des Stichprobenmedians nicht konsistent. Bootstraphingegen darf fur die Berechnung der Varianz des Stichprobenmediansverwendet werden (siehe [7, S. 69]).
17
Kapitel 3. Data-Resampling
5. Der delete-d-Jackknife-Algorithmus ist sehr effizient, da bei jedem aus-zufuhrenden Resampling-Schritt d Elemente weggelassen werden. Dasverringert die Rechenzeit erheblich.
6. Fur die Berechnung der Varianz und des Bias ist Jackknife aufgrunddes geringeren Rechenaufwands Bootstrap vorzuziehen.
7. Der Bootstrap-Algorithmus funktioniert nicht fur zeitlich abhangige(d. h. zeitlich korrelierte) Daten, falls die Korrelation in den Schatz-wert einfließt. Durch willkurliche Auswahl der Daten wurden die Kor-relationen zerstort werden.
8. Der Bootstrap-Algorithmus bietet den Vorteil, dass eine Approxima-tion der Verteilung von θ berechnet werden kann. Dadurch ist es mog-lich Konfidenzintervalle aufzustellen, die auch fur den Fall, dass θ nichtnormalverteilt ist, gultig sind.
18
Kapitel 4
Anwendungen
4.1 Lineare Regression
4.1.1 Mikroskop
Diese Anwendung bietet einen Einstieg in Jackknife und Bootstrap. MitJackknife wird die Varianz der Regressionsparameter bestimmt und mit Boot-strap ein 95% Konfidenzband der Regressionsgeraden aufgestellt. Alle Ergeb-nisse konnen auch analytisch berechnet werden.
Die Daten fur dieses Beispiel wurden aus der Laborubung Bestimmungder Vergroßerung und der Brennweiten eines Mikroskops der Lehrveran-staltung Grundpraktikum 1 (Mechanik, Warme, Optik) entnommen. DieMesswerte fur die Gesamtvergroßerung Vges in Abhangigkeit von der me-chanischen Tubuslange t sind in Tabelle 4.1 aufgestellt.
Tabelle 4.1: Gesamtvergroßerung Vges in Abhangigkeit von der mecha-nischen Tubuslange t/mm.
Nr. t/mm Vges
1 190 59
2 185 57
3 180 55
4 175 54
5 170 52
6 165 51
7 160 49
8 155 47
9 150 46
10 145 44
11 140 43
19
Kapitel 4. Anwendungen
Die Gesamtvergroßerung steht in linearem Zusammenhang mit der me-chanischen Tubuslange. Ziel ist es, eine lineare Regression der Daten durch-zufuhren und ein 95% Konfidenzband der Regressionsgeraden zu erzeugen.Das bedeutet, dass zu jedem Datenpunkt der x-Achse ein 95% Konfidenzin-tervall des zugehorigen y-Werts berechnet wird. Dadurch entsteht ein Bandin der die Regressionsgerade mit 95%iger Wahrscheinlichkeit verlauft. DieBerechnung der Varianz der Regressionsparameter erfolgt durch Anwendungvon delete-1-Jackknife und Bootstrap. Die Berechnung des Konfidenzbandserfolgt nur mit Bootstrap.
Dieses Beispiel ist als Einstieg fur die Anwendung von Jackknife undBootstrap zu betrachten. Alle Ergebnisse konnen auch uber direkte Metho-den, wenn auch mit komplizierterer Rechnung, erhalten werden.
Die lineare Regression basiert auf dem Geradenmodell
y = f(x) = a1 + a2 x . (4.1)
Im Geradenmodell (4.1) entspricht a1 dem Achsenabschnitt und a2 der Stei-gung. Die Durchfuhrung der linearen Regression und Berechnung der Vari-anzen ist z. B. in [1, Kap. 6] erklart. Als Schatzwerte fur a1 und a2 werdendie Beziehungen
a1 = 〈y〉 − 〈x〉 a2 (4.2)
a2 =〈xy〉 − 〈x〉〈y〉
〈x2〉 − 〈x〉2(4.3)
verwendet. Diese Beziehungen gelten, falls alle Datenpunkte yi die gleichenUnsicherheiten σ besitzen. Fur die Varianzen der Parameter gilt
var(a1) =σ2
n
〈x2〉
〈x2〉 − 〈x〉2
var(a2) =σ2
n
1
〈x2〉 − 〈x〉2,
wobei n die Anzahl der Datenpunkte und σ die Unsicherheit der y-Werte ist.σ wurde durch den Mittelwert der quadrierten Residuen berechnet ( [2, S.258]).
σ2 =1
n− 2
n∑
i=1
(
yi − (a1 + a2 xi))2
Jackknife und Bootstrap konnen auch zur Berechnung der Varianz desAchsenabschnitts a1 und der Varianz der Steigung a2 verwendet werden. Beider Anwendung des Jackknife-Algorithmus mussen jeweils n = 11 Berech-nungen des Achsenabschnitts und der Steigung durchgefuhrt werden, wobeifur die Auswertung von a1,i und a2,i der i-te Datenpunkt der mechanischen
20
Kapitel 4. Anwendungen
Tubuslange t und der i-te Datenpunkt der Gesamtvergroßerung Vges wegge-lassen werden. Mit den noch vorhandenen n−1 Daten werden die Parameterdes Geradenmodells nach (4.2) und (4.3) berechnet.
Bei Verwendung von Bootstrap wird die neue Stichprobe durch Ziehenmit Zurucklegen von n Elementen aus der ursprunglichen Datenmenge er-zeugt. Wird das i-te Element der Tubuslange t gezogen, so muss auch dasi-te Element der Gesamtvergroßerung Vges gezogen werden. Aufgrund die-ser Eigenschaft wird dieses Verfahren auch als Paired Bootstrap bezeich-net. Mathematisch bedeutet das, dass die Stichprobe durch eine empirischeVerteilungsfunktion F erzeugt wird, die ein zusammengehorendes Paar derTubuslange und der Gesamtvergroßerung liefert.
Die Ergebnisse der linearen Regression sind in Tabelle 4.2 zusammenge-fasst.
Tabelle 4.2: Ergebnisse der linearen Regression der Daten aus Tabelle 4.1mit dem Modell (4.1).
Direkte Methode Jackknife Bootstrap
a1 −2,164 −2,157 −2,165
ˆvar +1,050 +1,637 +1,118ˆstd +1,024 +1,279 +1,057
B – −0,007 +0,001
a2 +3,2000 e−001 +3,2004 e−001 +3,2005 e−001
ˆvar +0,0004 e−001 +0,0006 e−001 +0,0004 e−001
ˆstd +0,0618 e−001 +0,0773 e−001 +0,0637 e−001
B – −0,0004 e−001 −0,0005 e−001
Aus den Daten der Tabelle 4.2 kann abgelesen werden, dass sich dieSchatzwerte fur a1 und a2 in allen drei Fallen unter Beachtung der Stan-dardabweichung ˆstd ahneln. Außerdem ist der Bias B des Bootstrap- undJackknife-Schatzwerts klein im Vergleich zur Standardabweichung ˆstd undtragt damit nicht wesentlich zum Ergebnis bei.
Die Erstellung des Konfidenzbands der Regressionsgeraden erfolgt durchBerechnung der Verteilung der Gesamtvergroßerung Vges zu jeder Tubus-lange t. Dazu wird das Geradenmodell (4.1) fur die Tubuslange t mit allenWerten des Achsenabschnitts a∗1,i und der Steigung a∗2,i ausgewertet und dieerhaltenen Gesamtvergroßerungen aufsteigend sortiert. Die Berechnung desKonfidenzintervalls der Gesamtvergroßerung fur eine bestimmte mechani-sche Tubuslange t erfolgt durch Anwendung von (3.11). Wird das fur allegewunschten Tubuslangen (hier t ∈ [0, 200]) wiederholt, entsteht das Konfi-denzband der Regressionsgeraden.
Die grafische Darstellung der Daten erfolgt in Abbildung 4.1. Es sind dieDatenpunkte der Messungen, die Regressionsgerade auf Basis des Jackknife-
21
Kapitel 4. Anwendungen
Algorithmus und das 95% Konfidenzband, welches durch den Bootstrap-Algorithmus ausgewertet wurde, eingetragen. Die Darstellung der Geraden
0 50 100 150 200−10
0
10
20
30
40
50
60
70
Tubuslänge t/mm
Ges
amtv
ergr
ößer
ung
Vge
s
Gemessene DatenLineare Regression − Jackknife95% Konfidenzband − Bootstrap
Abbildung 4.1: Geradenmodell der Gesamtvergroßerung Vges in Abhangig-keit von der mechanischen Tubuslange t.
und des Konfidenzbands ist fur weitere Berechnungen bis t = 0 gewahlt.
Der Matlab-Quellcode fur dieses Beispiel ist im Anhang A.2 (Seite 43)angefuhrt.
4.2 Regression mit nichtlinearen Parametern
4.2.1 Modellfunktion y = a xbe−c x
Eine Modellfunktion wird durch Regression mit nichtlinearen Paramerternnach dem Levenberg-Marquardt-Algorithmus an kunstlich generierteDatenpunkte angepasst. Jackknife liefert fur die Varianz der Modellparame-ter vergleichbare Werte zur direkten Methode. Eine Berechnung des Konfi-denzbands der Regressionsfunktion erfolgt mit Bootstrap. Abschließend wirdder Vorteil von Resampling Methoden bei der Fehlerfortpflanzung aufgezeigt,der vor allem darin besteht, dass keine Ableitung ausgewertet werden muss,
22
Kapitel 4. Anwendungen
um die Varianz zu bestimmen.
In diesem Abschnitt wird die Regression mit nichtlinearen Parameternanhand der Modellfunktion
y = axb e−c x (4.4)
durchgefuhrt.Die Modellfunktion (4.4) erzeugt im Bereich x ∈ [1, 4] in aquidistantem
Abstand n = 150 Datenpunkte mit den Parametern a = 5, b = 3 undc = 2. Diese Datenpunkte werden durch Addition von zufalligen Wertender Normalverteilung (Mittelwert µ = 0, sehr große Standardabweichungσ = 0,2) verrauscht. Es folgt die Regression der als bekannt vorausgesetztenModellfunktion (4.4) mit den erzeugten Datenpunkten.
Da dieses Modell nicht linear ist, muss der Fit mit einem iterativenVerfahren durchgefuhrt werden. Es wurde der Levenberg-Marquardt-Algorithmus gewahlt, da dieser in vielen Fallen konvergiert und zudemschnell zur Konvergenz fuhrt. Die Parameter, die mit Hilfe des Levenberg-
Marquardt-Algorithmus berechnet werden, sind a, b und c. Zusatzlich lie-fert dieser Algorithmus einen Schatzwert der Kovarianzmatrix, die notwen-dig ist, um eine korrekte Fehlerfortpflanzung durchzufuhren (siehe Kap. 2).
Tabelle 4.3: Ergebnisse der nichtlinearen Regression der kunstlich erzeugtenDaten mit dem Modell (4.4).
Direkte Methode Jackknife Bootstrap
a +5,803 +5,595 +5,660
ˆvar +2,182 +1,975 +2,068ˆstd +1,477 +1,405 +1,438
B – +0,207 +0,142
b +3,155 +3,142 +3,157
ˆvar +0,343 +0,322 +0,313ˆstd +0,586 +0,568 +0,559
B – +0,014 −0,016
c +2,138 +2,129 +2,140
ˆvar +0,097 +0,088 +0,086ˆstd +0,311 +0,297 +0,293
B – +0,008 −0,002
Die Anwendung von Jackknife bzw. Bootstrap erfolgt analog zur linea-ren Regression. Zuerst werden Daten nach den jeweiligen Methoden gene-riert. Anschließend wird mit diesen Daten die nichtlineare Regression nachLevenberg-Marquardt durchgefuhrt. Es entstehen die Schatzwerte a(i),
23
Kapitel 4. Anwendungen
b(i) und c(i) aus denen die Varianz der Regressionsparameter nach dem je-weiligen Algorithmus berechnet werden konnen.
Das Ergebnis der Regression mit nichtlinearen Parametern ist in Ta-belle 4.3 dargestellt.
Aus Tabelle 4.3 ist ersichtlich, dass die Parameter a, b und c innerhalbder Standardabweichungen ˆstd bei den ursprunglichen Parametern a, b undc liegen, obwohl die Daten mit σ = 0,2 verrauscht wurden.
Die grafische Darstellung der nichtlinearen Regression erfolgt in Abbil-dung 4.2.
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
x
y
Generierte DatenExakte FunktionNichtlineare Regression − Jackknife95% Konfidenzband − Bootstrap
Abbildung 4.2: Darstellung der nichtlinearen Regression der Modellfunk-tion (4.4).
Abbildung 4.2 enthalt die erzeugten Daten, die exakte Funktion, die an-gepasste Modellfunktion – berechnet nach dem Jackknife-Algorithmus – unddas 95% Konfidenzband auf Basis des Bootstrap-Algorithmus. Die exakteFunktion ist die Modellfunktion (4.4) mit den Parametern a = 5, b = 3 undc = 2. Die durch Jackknife angepasste Funktion kann nicht genau der exak-ten Funktion entsprechen, da eine endliche Zahl an Datenpunkten mit derNormalverteilung verrauscht wurde. Ziel ist es aber, diese exakte Funktionbestmoglich anzunahern. In diesem Beispiel ist ersichtlich, dass das Konfi-denzband von den Daten gestutzt wird, da im Datenzentrum die Breite desKonfidenzbands minimal ist und nach außen hin zunimmt. In den Auslaufen
24
Kapitel 4. Anwendungen
der Funktion nimmt die Breite abermals ab, da durch die spezielle Form derModellfunktion (4.4) Einschrankungen der Funktionswerte entstehen.
Der Vergleich der Fehlerfortpflanzung mit der direkten Berechnung derVarianz uber Data-Resampling soll am Beispiel
θ = f(b, c) =c
b(4.5)
durchgefuhrt werden. Die theoretischen Grundlagen dazu wurden in Kap. 2und in Abschn. 3.4 eingefuhrt. Die Berechnung der Varianz von θ ist aufdrei Arten moglich:
1. Fehlerfortpflanzung ohne Berucksichtigung der Kovarianzen
2. Fehlerfortpflanzung mit Berucksichtigung der Kovarianzen
3. Berechnung der Varianz uber eine Data-Resampling Methode
Wird angenommen, dass die einzelnen Parameter, nach denen entwickeltwird, unabhangig sind, so gilt Gleichung (2.11). Angewandt auf die Funktionθ = f(b, c) gilt (4.6). Diese Beziehung spiegelt den Fall ohne Berucksichti-gung der Kovarianzen wider.
σ2θ ≈
(
∂f
∂b
)2
σ2b +
(
∂f
∂c
)2
σ2c (4.6)
Die Berucksichtigung der Kovarianzen erfolgt durch Hinzunahme gemischterAbleitungen, was allgemein in (2.10) beschrieben wurde. Diese Gleichungenthalt den Term der Kovarianz zwischen b und c, der mit σ2
bc bezeichnetwird. Fur θ = f(b, c) folgt
σ2θ ≈
(
∂f
∂b
)2
σ2b +
(
∂f
∂c
)2
σ2c + 2
(
∂f
∂b
)(
∂f
∂c
)
σ2bc . (4.7)
Gleichung (4.7) kann, wie in (2.12) beschrieben, in kompakter Matrixnota-tion geschrieben werden.
σ2θ ≈ JTCJ (4.8)
J =
(
∂f∂b∂f∂c
)
C =
(
σ2b σ2
bc
σ2bc σ2
c
)
J entspricht der Jacobi-Matrix und C der Kovarianzmatrix.Die Berechnung der Varianz mit Hilfe von Data-Resampling erfolgt auf
direktem Weg. Die Schatzwerte θ∗i werden aus den zuvor berechneten c∗i undb∗i durch
θ∗i =c∗ib∗i
∀i = 1, 2 . . . n
25
Kapitel 4. Anwendungen
festgelegt. Daraus konnen die Varianz ˆvar und der Bias B von θ abgeschatztwerden.
Die Ergebnisse der Fehlerfortpflanzung von θ sind in Tabelle 4.4 darge-stellt. Dabei entspricht σb und σc in (4.6), (4.7) und (4.8) der geschatztenStandardabweichung ˆstd der jeweiligen Variablen aus Tabelle 4.3.
Tabelle 4.4: Ergebnisse der Berechnung der Varianz von θ.
ohne Kov. mit Kov. Jackknife Bootstrap
θ +0,678 +0,678 +0,673 +0,672
ˆvar +0,027 +0,001 +0,001 +0,002ˆstd +0,163 +0,034 +0,033 +0,039
B – – +0,005 +0,006
Die Varianzen der Data-Resampling-Verfahren sind in der gleichen Gro-ßenordnung wie die Varianz der Fehlerfortpflanzung mit Kovarianz. Der Vor-teil der Data-Resampling-Methoden ist jedoch, dass keine Ableitungen desModells (4.5) notwendig sind. Die Varianz der Fehlerfortpflanzung ohne Be-rucksichtigung der Kovarianz ist jedoch um mehr als eine Zehnerpotenz gro-ßer als die anderen Varianzen. Dies liegt daran, dass die Parameter b undc voneinander abhangig sind, was durch die schrage Lage der Konturen inAbbildung 4.3, die den χ2-Wert von (4.5) bei Variation der Parameter b undc darstellt, ersichtlich ist.
b
c
3 3.05 3.1 3.15 3.2 3.25 3.3 3.35
1.95
2
2.05
2.1
2.15
2.2
2.25
2.3
Abbildung 4.3: χ2-Wert von (4.5) bei Variation der Parameter b und cdargestellt als Konturplot.
Mit Hilfe von Bootstrap kann auch die Verteilung von θ und somit einkorrektes Konfidenzintervall angegeben werden. Die Verteilung ist in Abbil-
26
Kapitel 4. Anwendungen
dung 4.4 grafisch dargestellt. Zum Vergleich ist auch die Normalverteilungmit der geschatzten Standardabweichung σθ (= ˆstd aus Tabelle 4.4) einge-zeichnet.
0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.90
5
10
15
θi*
f(θ i* )
Normalverteilung von θVerteilung von θ − Bootstrap
Abbildung 4.4: Vergleich der Verteilung von θ, berechnet durch Bootstrap,mit der Normalverteilung.
Der Matlab-Quellcode fur dieses Beispiel ist im Anhang A.3 (Seite 48)angefuhrt. Die Implementierung des Levenberg-Marquardt Algorithmuswurde von Ao.Univ.-Prof. Dr. Sormann durchgefuhrt. Eine theoretische Be-schreibung des Algorithmus findet sich in [8, Abschn. 4.5].
4.2.2 Fernrohr
Auch diese Aufgabe behandelt Regression mit nichtlinearen Parametern. Je-doch ist die Varianz der einzelnen Datenpunkte verschieden. Das bedeutet,dass Datenpunkte mit unterschiedlichem Gewicht in die Berechnung einflie-ßen. Zu groß abgeschatzte Varianzen der Datenpunkte, die zum Gewichtenangegeben werden aber nicht im Einklang mit den tatsachlichen Varianzender Datenpunkte stehen, werden von den Data-Resampling Methoden korri-giert.
Die Daten fur dieses Beispiel wurden aus der Laborubung Vergroße-rung eines Fernrohrs der Lehrveranstaltung Grundpraktikum 1 (Mechanik,Warme, Optik) entnommen. Die Vergroßerung des Fernrohrs V in Abhan-gigkeit von der Entfernung des Fernrohrs l/m vom beobachteten Objekt istin Tabelle 4.5 zusammengefasst.
Die Vergroßerung V wurde nach
V = V0M0
m0 l0
ml
M(4.9)
27
Kapitel 4. Anwendungen
Tabelle 4.5: Vergroßerung V in Abhangigkeit von der Entfernung l/m.
l . . . Abstand Messlatte-Okular; ∆l = 0,01m
M . . . Uberdeckende Einheiten auf der Messlatte; ∆M = 0,05m . . . Uberdeckende Einheiten des Okularmaßstabs; ∆m = 0,05V . . . Vergroßerung nach (4.9)
∆V . . . Unsicherheit der Vergroßerung nach (4.12)
Nr. l/m M m V ∆V
1 2, 00 31, 00 5, 00 29, 57 0, 77
2 2, 50 34, 00 4, 00 26, 96 0, 72
3 3, 00 33, 00 3, 00 25, 00 0, 72
4 3, 50 40, 00 3, 00 24, 06 0, 70
5 4, 00 38, 00 2, 50 24, 12 0, 75
6 4, 50 45, 00 2, 50 22, 92 0, 71
7 5, 00 47, 00 2, 30 22, 43 0, 72
8 5, 50 45, 00 2, 00 22, 40 0, 77
9 6, 00 42, 00 1, 70 22, 26 0, 84
10 6, 50 41, 00 1, 50 21, 80 0, 89
11 7, 00 47, 00 1, 60 21, 84 0, 85
12 7, 50 31, 00 1, 00 22, 18 1, 23
13 8, 00 41, 00 1, 20 21, 46 1, 03
14 8, 50 44, 00 1, 20 21, 25 1, 02
15 9, 00 46, 00 1, 20 21, 52 1, 03
16 9, 50 41, 00 1, 00 21, 24 1, 17
17 10, 00 44, 00 1, 00 20, 83 1, 15
berechnet, wobei V0, M0, m0 und l0 aus einer anderen Messmethode stam-men und folgende Werte besitzen
V0 = 25 ± 0,4
M0 = 33 ± 0,05
m0 = 3 ± 0,05
l0 = 3 ± 0,01 . (4.10)
Die Vergroßerung in Abhangigkeit von der Entfernung kann mit demModell (4.11) beschrieben werden.
V (l) = V∞l
l − (L+ f1)(4.11)
In (4.11) entspricht V∞ der Vergroßerung eines Objekts, das unendlich weitvom Fernrohr entfernt ist. L ist die Lange des Fernrohrs und f1 die Brenn-weite des Objektivs.
28
Kapitel 4. Anwendungen
Fur die Auswertung des Modells wird der Levenberg-Marquardt-Algorithmus verwendet, der die Parameter V∞ und ˆL+ f1 berechnet. In denLevenberg-Marquardt-Algorithmus fließt die Unsicherheit der Vergro-ßerung V ein. Fur die Fehlerrechnung von V wurde angenommen, dass dieUnsicherheiten der Parameter von (4.9) unabhangig und klein sind. Daherist es moglich, die Fehlerfortpflanzungsformel ohne Betrachtung der Kovari-anzen (2.11) zu verwenden und es folgt fur (4.9)
σV =
[(
∂V
∂V0σV0
)2
+
(
∂V
∂M0σM0
)2
+
(
∂V
∂m0σm0
)2
+
(
∂V
∂l0σl0
)2
+
(
∂V
∂MσM
)2
+
(
∂V
∂mσm
)2
+
(
∂V
∂lσl
)2] 1
2
. (4.12)
Die Standardabweichungen werden aus (4.10) entnommen und furM , m undl analog zu M0, m0 und l0 geschatzt. Die partiellen Ableitungen werden vonder Funktion (4.9) durchgefuhrt. Die Standardabweichungen σV sind fur dieeinzelnen Datenpunkte verschieden und in der letzten Spalte in Tabelle 4.5angegeben.
Nachdem die Daten nach dem Jackknife- bzw. Bootstrapverfahren selek-tiert wurden, muss fur jeden Datensatz eine Regression durchgefuhrt werden.Dadurch entstehen die Werte V∞,(i) und ( ˆL+ f1)(i) aus denen die Varianzder Parameter berechnet werden kann. Mit Bootstrap kann, wie bereits be-schrieben, zusatzlich das Konfidenzband der Regressionsfunktion berechnetwerden.
Das Ergebnis der Regression mit nichtlinearen Parametern ist in Ta-belle 4.6 angegeben.
Tabelle 4.6: Ergebnisse der nichtlinearen Regression der Daten aus Ta-belle 4.5 mit dem Modell (4.11).
Direkte Methode Jackknife Bootstrap
V∞ +19,673 +19,649 +19,654
ˆvar +0,111 +0,007 +0,010ˆstd +0,333 +0,086 +0,098
B – +0,024 +0,019
ˆL+ f1 +0,6664 +0,6712 +0,6707
ˆvar +0,0020 +0,0001 +0,0003ˆstd +0,0443 +0,0104 +0,0183
B – −0,0048 −0,0043
Auffallend am Ergebnis ist, dass die Varianz der Parameter aus den Data-
29
Kapitel 4. Anwendungen
Resampling Methoden weit geringer ist, als die Varianz der Parameter ausdem direkten Fit. Der Unterschied zwischen den Varianzen entsteht dadurch,dass die Fehler von V0, M0, m0 und l0 im Vergleich zur Streuung der Da-ten zu groß abgeschatzt wurden. Deshalb sind die nach (4.12) berechnetenUnsicherheiten der Vergroßerung σV zu groß. Die zu großen Unsicherhei-ten σV werden der Fit-Routine mitgeteilt und fließen in die Berechnung deszu minimierenden Werts χ2 ein, was wiederum eine großere Varianz der Re-gressionsparameter zur Folge hat. Bei der Anwendung von Data-Resamplingtragt nur die Streuung von V∞,(i) und ( ˆL+ f1)(i) zur Varianz von V∞ und
ˆL+ f1 bei. Das bedeutet, dass die Varianz der Streuung der Daten ange-passt ist. Daher ist die Data-Resampling Varianz kleiner als die Varianz, diemit der direkten Methode und den zu groß abgeschatzten Unsicherheiten σV
berechnet wurde. Die direkte Methode und die Data-Resampling Methodenwurden einen ahnlichen Schatzwert der Varianz liefern, falls die Unsicher-heiten σV im Einklang mit der Steuung der Daten waren. Das ware derFall, wenn ca. ein Drittel der Datenpunkte mit ihrer StandardabweichungσV außerhalb der Regressionskurve liegen wurden.
Abbildung 4.5 enthalt die grafische Auswertung der Ergebnisse von Ta-belle 4.6. In Abbildung 4.5 sind die berechneten Vergroßerungen mit Fehler-
2 4 6 8 10 12 1418
20
22
24
26
28
30
32
Entfernung l/m
Ver
größ
erun
g V
Gemessene DatenNichtlineare Regression − Jackknife95% Konfidenzband − Bootstrap
Abbildung 4.5: Darstellung der Vergroßerung V in Abhangigkeit von derEntfernung l/m mit Modell (4.11).
30
Kapitel 4. Anwendungen
balken, die Regressionskurve nach der Modellfunktion (4.11) und das 95%Konfidenzband eingetragen. Es ist deutlich erkennbar, dass das Konfidenz-band nicht symmetrisch um die Regressionskurve liegt. Daher ist es in diesemFall notwendig, zur Berechnung des Konfidenzbands Bootstrap zu verwen-den, um korrekte Ergebnisse zu erhalten.
Der Matlab-Quellcode fur dieses Beispiel ist im Anhang A.3 (Seite 51)angefuhrt.
4.3 Integrierte Autokorrelationszeit1
Dieser Abschnitt liefert eine Einfuhrung in die Behandlung korrelierter Da-ten. Mit Hilfe von Jackknife wird ein Bias-reduzierter Schatzwert der inte-grierten Autokorrelationszeit τO,int ermittelt. Die integrierte Autokorrelati-onszeit gibt den Grad der Abhangigkeit der Daten an. Außerdem wird einSchatzwert ihrer Varianz berechnet. Erst durch die Selektion der Daten nachdem Jackknife-Algorithmus ist es moglich, die Varianz der integrierten Au-tokorrelationszeit zu berechnen.
Bei korrelierten Daten ist die Berechnung der Varianz σ2O
einer Obser-vablen O nicht trivial durchzufuhren. Dies erklart sich dadurch, dass in kor-relierten Daten weniger Informationen stecken als in unkorrelierten Daten.
Korrelierte Daten werden durch Zeitreihen beschrieben. Eine endlicheZeitreihe, wie sie hier betrachtet wird, besteht aus einer Folge von Zeiten
ti = i∆t mit i = 1, 2, . . . , N und ∆t > 0 ,
wobei jeder Zeit ti eine Zufallsvariable Xi zugeordnet wird [3, Kap. 1]. Dasist ein stochastischer Prozess. Jede Realisierung des stochastischen Prozessesist durch eine Folge Oi mit i = 1, 2, . . . , N gegeben, wobei Oi = O(Xi) deri-ten Messung der Observable O zur Zeit ti entspricht. Eine Zeitserie konntebeispielsweise so aussehen wie in Abbildung 4.6.
Ein Schatzwert des Erwartungswerts ist auch fur korrelierte Daten durchden arithmetischen Mittelwert
O =1
N
N∑
i=1
Oi (4.13)
gegeben. Das ist eine Eigenschaft von Zeitreihen, die aus der Markov-Ketten-Monte-Carlo (MCMC – Markov-Chain-Monte-Carlo) Methode resultiert [3,Abschn. 1.7.3]. Die Varianz von unkorrelierten Daten kann durch
σ2Oi
=1
N2
N∑
i=1
(〈O2i 〉 − 〈Oi〉
2) (4.14)
1vgl. [5, Kap. 4]
31
Kapitel 4. Anwendungen
0 5 000 10 000 15 000 20 000 25 000 30 000−2.5
−2
−1.5
−1
−0.5
0
0.5
1
ti
Oi
Abbildung 4.6: Darstellung einer Zeitserie nach dem Modell (4.18) mitτO,int ≈ 5 000.
berechnet werden. Die Varianz (4.14) muss fur korrelierte Daten um einenFaktor c vergroßert werden, wobei c ∝
τO,int
N.N entspricht der Gesamtanzahl
der Messungen. τO,int heißt integrierte Autokorrelationszeit und ist ein Maßfur die Zeitdauer einer zusammengehorenden Region in einer Zeitserie. InAbbildung 4.6 ist τO,int ≈ 5 000. In dieser Aufgabe soll τO,int berechnetwerden. Die Varianz der korrelierten Daten ist durch
σ2O
= 2τO,int
Nσ2Oi
(4.15)
gegeben. Dabei wird die integrierte Autokorrelationszeit durch
τO,int =1
2+
N∑
k=1
A(k) (4.16)
mit der Autokorrelationsfunktion
A(k) =〈OiOi+k〉 − 〈Oi〉〈Oi+k〉
〈O2i 〉 − 〈Oi〉2
(4.17)
berechnet.Die Herleitung der Gleichungen (4.15) – (4.17) ist in [5, S. 431] durchge-
fuhrt. Ziel fur weiterfuhrende Berechnungen ist die integrierte Autokorrela-tionszeit – und wenn moglich auch ihre Unsicherheit – zu bestimmen.
Fur ein Beispiel wird eine korrelierte Zeitserie ei der Form
e0 = e′0
ei = ρ ei−1 +√
1 − ρ2 e′i, i ≥ 1 (4.18)
konstruiert, wobei ρ ein frei wahlbarer so genannter Korrelationskoeffizientist. Dabei sind e′i normalverteilte Zufallsvariablen mit Erwartungswert
〈e′i〉 = 0 (4.19)
32
Kapitel 4. Anwendungen
und Kovarianzmatrix
〈e′i e′j〉 = δij . (4.20)
Nach Auflosung der Iteration in (4.18) folgt
ek = ρk e0 +√
1 − ρ2
k∑
l=1
ρk−l e′l = ρk e′0 +√
1 − ρ2
k∑
l=1
ρk−l e′l . (4.21)
Unter diesen Annahmen ist die exakte Berechnung der Autokorrelations-funktion nach (4.17) moglich.
〈ek〉(4.21)= ρk 〈e′0〉 +
√
1 − ρ2
k∑
l=1
ρk−l 〈e′l〉
(4.19)= 0 (4.22)
〈e2k〉(4.21)= 〈
[
ρk e′0 +√
1 − ρ2
k∑
l=1
ρk−l e′l
]2〉
= ρ2k 〈e′02〉 + 2
√
1 − ρ2
k∑
l=1
ρk−l 〈e′l e′0〉 ρ
k
+(1 − ρ2)〈[
k∑
l=1
ρk−l e′l
]2〉
(4.20)= ρ2k + (1 − ρ2)
k∑
l=1
k∑
m=1
ρk−l ρk−m 〈e′l e′m〉
(4.20)= ρ2k + (1 − ρ2)
k∑
l=1
(ρ2)k−l = ρ2k + (1 − ρ2)k−1∑
l′=0
(ρ2)l′
= ρ2k + (1 − ρ2)1 − (ρ2)k
1 − ρ2= ρ2k + 1 − (ρ2)k
= 1 (4.23)
〈ei ei+k〉 = 〈e0 ek〉
= 〈e′0 (ρk e′0 +√
1 − ρ2
k∑
l=1
ρk−l e′l)〉
= ρk 〈e′02〉 +
√
1 − ρ2
k∑
l=1
ρk−l 〈e′0 e′l〉
(4.20)= ρk (4.24)
33
Kapitel 4. Anwendungen
Die Umformung im ersten Schritt in (4.24) darf durchgefuhrt werden,da der Erwartungswert 〈eiei+k〉 nur vom Abstand zwischen ei und ei+k ab-hangt und alle ei den selben Erwartungswert besitzen. Damit gilt fur dieAutokorrelationsfunktion (4.17) die Beziehung
A(k) = ρk ≡(
e− 1
τexp
)k
= e− k
τexp . (4.25)
Die Autokorrelationsfunktion besteht in diesem einfachen Fall nur auseiner Exponentialfunktion. In realistischen Fallen existiert immer eine Uber-lagerung von mehreren Exponentialfunktionen. Fur große k kann aber ange-nommen werden, dass die Exponentialfunktion mit dem geringsten Gefallegegenuber den anderen sehr stark dominiert, so dass fur große k nur mehrdiese betrachtet werden muss.
0 200 400 600 800 1000−4
−3
−2
−1
0
1
2
3
t
e
unkorrelierte Daten
0 200 400 600 800 1000−3
−2
−1
0
1
2
3
t
e
τexp
=10
(a) (b)
0 200 400 600 800 1000−2
−1.5
−1
−0.5
0
0.5
1
1.5
2
2.5
3
t
e
τexp
=50
−5 0 50
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45unkorrelierte Datenτexp
=10
τexp
=50
Normalverteilung
(c) (d)
Abbildung 4.7: In (a) sind die ersten 1 000 Werte der Simulation von un-korrelierten Daten dargestellt. Die Abbildungen (b) und (c) enthalten dieersten 1 000 Datenpunkte der korrelierten Zeitreihe (4.21) fur τexp = 10 undτexp = 50. Abbildung (d) zeigt das Histogramm uber alle 100 000 Messwertefur die Falle (a), (b) und (c) im Vergleich zur standardisierten Normalver-teilung.
In diesem Beispiel werden 100 000 aufeinanderfolgende Werte fur ek nachGleichung (4.21) berechnet. Dabei wird einerseits τexp = 10 und andererseits
34
Kapitel 4. Anwendungen
τexp = 50 verwendet. Die ersten 1 000 Werte dieser Zeitreihe werden in Ab-bildung 4.7 (b) und (c) mit der unkorrelierten Zeitreihe (a) verglichen. Ausdem Histogramm in Abbildung 4.7 (d) ist ersichtlich, dass in allen drei Fallendie Daten, wie es in (4.22) und (4.23) berechnet wurde, der standardisiertenNormalverteilung genugen.
Im nachsten Schritt wird die Autokorrelationsfunktion aus den Datenberechnet und mit (4.25) verglichen. Ein Schatzwert fur die Autokorrela-tionsfunktion entsteht durch das Ersetzen der Erwartungswerte in (4.17)durch Mittelwerte.
A(k) =1
N ′
∑N ′
i=1 Oi Oi+k − 1N ′
∑N ′
i=1 Oi1
N ′
∑N ′
i=1 Oi+k
1N
∑Ni=1 O
2i −
(
1N
∑Ni=1 Oi
)2 mit N ′ = N − k
(4.26)Die Autokorrelationsfunktion und ihre Varianz ist fur τexp = 10 in Ab-
bildung 4.8 (a) dargestellt. Es ist darauf zu achten, dass die Werte der Au-tokorrelationsfunktion logarithmisch aufgetragen sind.
0 20 40 60 80 10010
−5
10−4
10−3
10−2
10−1
100
k
log(
A(k
))
τexp
=10 − Schätzwert
τexp
=10 − Exakt
0 20 40 60 80 1000
2
4
6
8
10
12
kmax
τ int(k
max
)
τexp
=10 − Schätzwert
τexp
=10 − Exakt
(a) (b)
Abbildung 4.8: In (a) ist die theoretisch berechnete und die aus den Da-ten geschatzte Autokorrelationsfunktion im halblogarithmischen Diagrammdargestellt. Abbildung (b) enthalt den Vergleich zwischen der theoretischberechneten und der aus den Daten geschatzen integrierten Autokorrelati-onszeit.
Die theoretische Funktion stimmt fur kleine k gut mit den aus den Da-ten berechneten Schatzwerten uberein. Fur k > 30 steigt die Varianz jedochdrastisch an, so dass die Datenpunkte unbrauchbar werden. Die Varianzder Autokorrelationsfunktion kann mit delete-d-Jackknife einfach berechnetwerden. Sonst ist eine Berechnung der Varianz kaum moglich. Bei der An-wendung von delete-d-Jackknife werden die durch (4.21) simulierten Datenin NBlock Blocke unterteilt, wobei d = N
NBlock� τint gelten muss. Hier wurde
NBlock = 20 gewahlt und damit gilt d = 100 00020 = 5000 � τint (siehe (4.27)).
Durch Weglassen eines Blocks werden N − d Datenpunkte ausgewahlt, ausdenen die Autokorrelationsfunktion nach (4.26) berechnet werden kann. Die
35
Kapitel 4. Anwendungen
Daten des Blocks werden bei der Implementierung nicht tatsachlich ausge-lassen, sondern auf Null gesetzt, damit keine falschen Korrelationen berech-net werden. Es entstehen i = 1, . . . NBlock unterschiedliche Werte A(i)(k) zujedem k. Die Varianz der Autokorrelationsfunktion kann mit dem Jackknife-Varianzschatzwert (3.6) fur ein bestimmtes k berechnet werden. Wird dieBerechnung der Varianz fur jedes k wiederholt, entsteht ein Konfidenzbandder Autokorrelationsfunktion, das zwei Standardabweichungen breit ist.
Mit (4.16) kann die integrierte Autokorrelationszeit einerseits aus derexakten Autokorrelationsfunktion und andererseits aus der geschatzten Au-tokorrelationsfunktion berechnet werden.
τint =1
2+
kmax∑
k=1
A(k) =1
2+
kmax∑
k=1
ρk
(4.25)= −
1
2+
kmax∑
k=0
e− k
τexplimkmax→∞
= −1
2+
1
1 − e− 1
τexp
τexp�1≈ τexp (4.27)
τint =1
2+
kmax∑
k=1
A(k) (4.28)
Abbildung 4.8 (b) enthalt den grafischen Vergleich von (4.27) und (4.28).Dabei ist ersichtlich, dass die zweite Gleichung fur die Berechnung der inte-grierten Autokorrelationszeit aufgrund der großen Unsicherheit der geschatz-ten Autokorrelationsfunktion, die bei großen k vorliegt, nicht brauchbar ist.
Gleichung (4.28) konnte verwendet werden, falls A(k) bei großen k ge-nauer bestimmt ware. Um das zu erreichen, wird eine Regression durch-gefuhrt, bei der die zuverlassigen Daten einfließen (in diesem Fall bis k ≈20). Mit dieser Regressionsfunktion konnen die Autokorrelationskoeffizien-ten A(k) fur große k besser bestimmt werden. Wird der Logarithmus derDaten betrachtet, kann ein einfaches Geradenmodell als Regressionsfunk-tion verwendet werden.
log(Afit(k)) = a1 + a2 k (4.29)
Die Funktion Afit beschreibt die Autokorrelationsfunktion fur große k
besser als A(k) und es gilt fur die integrierte Autokorrelationszeit
τint =1
2+
∞∑
k=1
Afit(k)(4.29)= −
1
2+
∞∑
k=0
ea1+a2 k
= −1
2+ ea1
∞∑
k=0
(
ea2
)k
= −1
2+ ea1
1
1 − ea2
. (4.30)
36
Kapitel 4. Anwendungen
In vielen Fallen besteht die Autokorrelationsfunktion aus einer Uberla-gerung mehrerer Exponentialfunktionen. Dann liefert der Logarithmus derDaten keine einfache Gerade. Jedoch gibt es ein k1, ab dem die Autokorrela-tionsfunktion annahernd linear verlauft. Das kann angenommen werden, dajene Exponentialfunktion aus der Uberlagerung aller Exponentialfunktionen,die das geringste Gefalle besitzt, ab einem k1 gegenuber den anderen domi-niert. Dadurch verhalt sich die Autokorrelationsfunktion ab k1 so, als obnur noch eine Exponentialfunktion vorliegen wurde. Gleichung (4.30) musssomit durch die Summe
∑k1−1k=1 A(k) erweitert werden. Die Regression darf
nur mit den zuverlassigen Daten ab dem k1-ten Datenpunkt erfolgen und zuτint beitragen. Es ergibt sich
τint =1
2+
k1−1∑
k=1
A(k) +
∞∑
k=k1
Afit(k) .
Ein Schatzwert der Varianz der integrierten Autokorrelationszeit kannwiederum durch Anwendung von Jackknife erhalten werden. Dazu werdendie i = 1 , . . . NBlock Autokorrelationsfunktionen A(i)(k), die zur Berechnungder Varianz der Autokorrelationsfunktion generiert wurden, verwendet. Mitden zuverlassigen Autokorrelationskoeffizienten A(i)(k) bis k ≈ 20 der i-tenAutokorrelationsfunktion wird eine Regression durchgefuhrt. Dadurch kanndie i-te Autokorrelationsfunktion fur große k besser bestimmt werden undes kann mit (4.30) ein Schatzwert τint,(i) berechnet werden. Nachdem τint,(i)
fur alle i = 1, . . . NBlock berechnet wurde, ist es moglich die Unsicherheitvon τint durch Anwenden des Jackknife-Varianzschatzwerts zu bestimmen.
In Tabelle 4.7 ist ein Vergleich zwischen dem mit (4.27) theoretisch be-rechneten Wert und dem uber den Jackknife-Algorithmus berechneten Wertder integrierten Autokorrelationszeit gegeben. Es ist zu erkennen, dass diebeiden Werte außerordentlich gut ubereinstimmen.
Tabelle 4.7: Ergebnisse der Berechnung der integrierten Autokorrelations-zeit.
theoretischer Wert Jackknife
τint +10,008 +10,024
ˆvar – +0,102ˆstd – +0,320
B – +0,020
Der Matlab-Quellcode fur dieses Beispiel ist im Anhang A.4 (Seite 54)angefuhrt.
37
38
Kapitel 5
Zusammenfassung
Die Anwendung von Data-Resampling Methoden erleichtert in vielen Fal-len eine Auswertung von Daten. Standardabweichungen und Konfidenzin-tervalle von Schatzwerten lassen sich ohne komplizierte Herleitung des zu-grunde liegenden Modells berechnen. Der nicht-parametrische Fall von Data-Resampling zeigt, dass es nicht notwendig ist, Annahmen uber die Verteilungder Daten zu treffen. Es besteht auch die Moglichkeit, halb-parametrischeund parametrische Modelle mit Data-Resampling zu analysieren.
Eine weitere Anwendung von Data-Resampling besteht in der Durchfuh-rung von Hypothesentests, die ausfuhrlich in [2] und [7] beschrieben wird.
Data-Resampling kann in verschiedenen Anwendungsgebieten verwendetwerden. Diese Gebiete umfassen lineare Regression und nichtlineare Regres-sion. Außerdem ist eine Erweiterung auf die Untersuchung von Zeitserienund von anderen abhangigen Daten moglich.
In dieser Arbeit wurde Jackknife zur Berechnung von Standardabwei-chungen gegenuber Bootstrap vorgezogen, da dieser Algorithmus wenigerRechenschritte benotigt, um zum aquivalenten Ergebnis zu gelangen. Boot-strap wurde hingegen zur Berechnung von Konfidenzintervallen, die nichtnotwendig normalverteilt sein mussen, verwendet. Ein Vorteil von Boot-strap gegenuber Jackknife ist, dass Bootstrap eine Verteilung des Schatz-werts generiert, aus der durch herkommliche Vorgehensweise Varianz, Biasund Konfidenzintervalle berechnet werden konnen.
Ganz besonders nutzlich ist die Anwendung von Data-Resampling Me-thoden, wenn die Daten in vielen aufeinanderfolgenden Schritten manipuliertund zusammengefasst werden. In solchen Fallen kann direkt auf die Varianzdes Schatzwerts geschlossen werden, indem der Jackknife- oder Bootstrap-Algorithmus auf die zugrundeliegenden Daten angewendet wird. Ein Beispieldafur ist die Fehlerfortpflanzung.
39
40
Anhang A
Matlab-Quellcode
A.1 Stichprobenmittelwert – Bootstrap
Listing A.1: Stichprobenmittelwert – Bootstrap
1 % Estimates the mean and variance of
2 % a normal distributed sample.
3 %
4 % Author: Knap Michael
5
6 % reset
7 clc;
8 clear all;
9 close all;
10
11 % generate data
12 n=100;
13 sigma =0.1;
14 x=normrnd (0,sigma ,1,n);
15
16 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
17 %% mean (X)
18 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
19
20 % maximum likelihood estimator
21 xQuer=mean (x);
22
23 % calculate variance from residuals
24 varxQuerResidual =1/(n*(n -1))*sum((x-xQuer ).^2);
25
26 % apply Bootstrap
27 % choose randomly (uniform distributed ) from x
28 B=999;
29 index=ceil (n*rand (B,n));
30 xRand=x(index);
31
32 % compute mean
33 xQuerAst =mean (xRand ,2);
41
Anhang A. Matlab-Quellcode
34 xQuerHat =mean (xQuerAst )
35
36 % compute bias -> approximately 0 is expected
37 biasxQuerPB =xQuerHat -xQuer
38 xQuerAdjustedPB =xQuer - biasxQuerPB
39
40 % compute var -> approximately var/n
41 varxQuerPB =1/(B -1)* sum (( xQuerAst -xQuerHat ).^2)
42 % calculate the mse
43 msexQuerPB =varxQuerPB + biasxQuerPB ^2
Listing A.2: Konfidenzintervalle des Stichprobenmittelwerts – Bootstrap
1 % Estimates the 90% confidence interval of the mean
2 % of a normal distributed sample.
3 %
4 % Author: Knap Michael
5
6 % reset
7 clc;
8 clear all;
9 close all;
10
11 % generate data
12 n=100;
13 sigma =0.1;
14 x=normrnd (0,sigma ,1,n);
15
16 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
17 %% quantiles q1 , q2
18 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
19
20 % maximum likelihood estimator
21 xQuer=mean (x);
22
23 % calculate variance from residuals
24 varxQuerResidual =1/(n*(n -1))*sum((x-xQuer ).^2);
25
26 % Bootstrap
27 % choose randomly (uniform distributed ) from x
28 B=[19 ,39 ,99 ,199 ,499 ,699 ,899 ,999];
29 numberOfRepetitions =4;
30 q1=zeros(length(B), numberOfRepetitions );
31 q2=zeros(length(B), numberOfRepetitions );
32 for (k=1: length(B))
33 for (l=1: numberOfRepetitions )
34
35 % realise Bootstrap algorithm
36 for m=1:B(k)
37 index=ceil (length(x)* rand (1, length(x)));
38 XStern=x(index);
39 tQ(m)=mean (XStern );
42
Anhang A. Matlab-Quellcode
40 end
41
42 tQ=sort (tQ);
43 q1(k,l)=tQ((B(k )+1)*0.05);
44 q2(k,l)=tQ((B(k )+1)*0.95);
45 end
46 end
A.2 Lineare Regression
Listing A.3: Lineare Regression – Mikroskop
1 % Linear regression for microscope data .
2 %
3 % Author: Knap Michael
4
5 % reset
6 clc;
7 clear all;
8 close all;
9
10 % measured data
11 x=[190 ,185 ,180 ,175 ,170 ,165 ,160 ,155 ,150 ,145 ,140];
12 y=[59 ,57 ,55 ,54 ,52 ,51 ,49 ,47 ,46 ,44 ,43];
13
14 % data size
15 n=length(x);
16
17 xModel =[0 ,200];
18 xPrediction = [xModel (1):0.1: xModel(length(xModel ))];
19
20 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
21 %%% linear regression
22 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
23 Cxx =1/n*sum(x.^2) -(1/n*sum(x))^2;
24 Cyy =1/n*sum(y.^2) -(1/n*sum(y))^2;
25 Cxy =1/n*sum(x.*y)-(1/n*sum(x)*1/ n*sum(y));
26
27 yQuer=mean (y);
28 xQuer=mean (x);
29
30 a2=Cxy/Cxx;
31 a1=yQuer -xQuer*a2;
32
33 % calculate variance from residuals
34 sSquare =1/(n -2)* sum((y-(a1+a2*x)).^2);
35 varA1Residual =sSquare *(1/ n+xQuer ^2/( n*Cxx ));
36 varA2Residual =sSquare /(n*Cxx);
37 covarA1A2Residual =-sSquare /n*sum(x)/(n*Cxx);
38
39 y_model=a2*xModel+a1;
43
Anhang A. Matlab-Quellcode
40
41 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
42 %%% Jackknife
43 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
44 xQueri =1/(n -1)*(n*xQuer -x);
45 yQueri =1/(n -1)*(n*yQuer -y);
46
47 Cxxi =zeros(1,n);
48 Cxyi =zeros(1,n);
49 for k=1:n
50 cxTemp =0;
51 cyTemp =0;
52 for (l=1:n)
53 if (k~=l)
54 cxTemp=cxTemp+x(l)^2;
55 cyTemp=cyTemp+x(l)*y(l);
56 end;
57 end;
58 Cxxi (k)=1/(n -1)*( cxTemp)-( xQueri(k))^2;
59 Cxyi (k)=1/(n -1)*( cyTemp)-( xQueri(k)* yQueri(k));
60 end;
61
62 a2i=Cxyi ./ Cxxi;
63 a1i=yQueri -xQueri .*a2i;
64
65 a1Dot=mean (a1i );
66 a2Dot=mean (a2i );
67
68 BiasA1=(n -1)*(a1Dot -a1)
69 BiasA2=(n -1)*(a2Dot -a2)
70
71 a1Adjusted =a1 -BiasA1
72 a2Adjusted =a2 -BiasA2
73
74 y_modelAdjusted = a2Adjusted *xModel+a1Adjusted ;
75
76 % calculate variance
77 varA1Jack =(n -1)/n*sum ((a1i -a1Dot ).^2)
78 varA2Jack =(n -1)/n*sum ((a2i -a2Dot ).^2)
79
80
81 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
82 %%% Bootstrap
83 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
84 B=999;
85 % choose randomly (uniform distributed ) from x and y
86 index=ceil (n*rand (B,n));
87 xRand=x(index);
88 yRand=y(index);
89
90 % compute slope and intercept
91 Cxxi =1/n.* sum(xRand .^2 ,2) -(1/ n.* sum(xRand ,2)).^2;
92 Cxyi =1/n.* sum(xRand.*yRand ,2) -...
93 (1/n.* sum(xRand ,2).*1/ n.* sum(yRand ,2));
44
Anhang A. Matlab-Quellcode
94
95 yQueri=mean (yRand ,2);
96 xQueri=mean (xRand ,2);
97
98 a2Ast=Cxyi ./ Cxxi;
99 a1Ast=yQueri -xQueri .* a2Ast;
100 xTubusAst =-a1Ast./ a2Ast;
101
102 a1Hat=mean (a1Ast)
103 a2Hat=mean (a2Ast)
104
105 % compute bias -> approximately 0 is expected
106 biasA1PB =a1Hat -a1
107 biasA2PB =a2Hat -a2
108
109 a1AdjustedPB =a1 -biasA1PB
110 a2AdjustedPB =a2 -biasA2PB
111
112 % compute variance
113 varA1PB =1/(B -1)* sum(( a1Ast -a1Hat ).^2)
114 varA2PB =1/(B -1)* sum(( a2Ast -a2Hat ).^2)
115
116 y_modelBootstrap = a2AdjustedPB *xModel+a1AdjustedPB ;
117
118 % calculate the mse
119 mseJackA1 =varA1Jack +BiasA1 ^2
120 msePBA1=varA1PB+biasA1PB ^2
121
122 mseJackA2 =varA2Jack +BiasA2 ^2
123 msePBA2=varA2PB+biasA2PB ^2
124
125 % prediction interval - Bootstrap
126 yConfidencePB_i =repmat(a1Ast ,1, length(xPrediction ))+...
127 repmat(a2Ast ,1, length(xPrediction )).*...
128 repmat(xPrediction , length(a2Ast ),1);
129 yConfidencePB_i =sort ( yConfidencePB_i );
130 yConfidencePBTop = yConfidencePB_i ((B+1)*0.025 ,:);
131 yConfidencePBBottom = yConfidencePB_i ((B+1)*0.975 ,:);
132 y_area=[ yConfidencePBBottom ,fliplr(yConfidencePBTop )];
133 x_area=[ xPrediction ,fliplr(xPrediction )];
A.3 Regression mit nichtlinearen Parametern
Listing A.4: Implementierung des Levenberg-Marquardt-Algorithmus
1 function [a,sd_a ,var ,covar ,c,chisq] = ...
2 mrqmin(funcs ,x,y,sig ,ndata ,ma ,a,eps);
3
4 % Routine zur Durchfuehrung einer Gauss -Newton -Marquardt -
5 % Iteration nach dem Skriptum , Abschnitt 4.5.5
6
45
Anhang A. Matlab-Quellcode
7 % H. Sormann November 2006
8
9 % BESCHREIBUNG DER PARAMETER :
10
11 % Eingabe -Parameter :
12 % funcs Filename fuer die externe Modellfunktion
13 % x Zeilenvektor der x-Komponenten der geg. Punkte
14 % y Zeilenvektor der y-Komponenten der geg. Punkte
15 % sig Zeilenvektor der Standardabweichungen der geg. Punkte
16 % ndata Zahl der geg. Punkte
17 % ma Zahl der Fit -Parameter
18 % a Zeilenvektor der guessed values fuer die
19 % Fit -Paramter
20 % eps relativer Fehler der Fit -Parameter , der
21 % zu unterschreiten ist
22 %
23 % Ausgabe -Parameter :
24 % a Zeilenvektor der gefitteten Parameter
25 % sd_a Zeilenvektor der Standardabweichungen der Parameter
26 % var Varianz des Fits Anfang und Ende des Sollbereichs
27 % covar Kovarianz - Matrix (nicht normiert )
28
29 % Dieses Programm ruft die folgende Routine mrqcof auf.
30 % Diese Routine arbeitet genau wie im Skriptum (S. 125f)
31 % beschrieben .
32
33 % Die Routine mrqcof ruft ihrerseits die Routine funcs auf
34
35 alamda =0.0003;
36 [alpha ,beta ,chisq] = mrqcof(funcs ,x,y,sig ,ndata ,ma ,a);
37
38 beta =beta ’;
39 ochisq=chisq;
40
41 aalt =a;
42 for iter =1:100
43 normal=alpha;
44 for i=1: ma
45 normal(i,i)= alpha(i,i)*(1+ alamda);
46 end
47
48 augmented =[ normal beta ];
49 cr=rref (augmented );
50 da=cr(:,ma +1);
51
52 atry =a+da ’;
53 [normal ,vektor ,chisq] = ...
54 mrqcof(funcs ,x,y,sig ,ndata ,ma ,atry );
55
56 if(chisq <= ochisq)
57 alamda=alamda /5;
58 ochisq=chisq;
59 alpha=normal;
60 a=atry ;
46
Anhang A. Matlab-Quellcode
61 beta =vektor ’;
62
63
64 % Fehlerdiagnostik :
65
66 delmax =0;
67 for j=1: ma
68 del=abs((a(j)-aalt (j))/a(j));
69 if(del > delmax)
70 delmax=del;
71 end
72 end
73 if(delmax < eps)
74 fprintf(’Genauigkeit wurde erreicht !!\n’);
75 % Berechnung der Varianz :
76 var =[ chisq/(ndata -ma) 1-sqrt (2/( ndata -ma)) ...
77 1+ sqrt (2/( ndata -ma ))];
78 c=( ndata -ma)/ chisq;
79 % Berechnung der Kovarianz -Matrix:
80 alamda =0.0;
81 [alpha ,beta ,chisq] = mrqcof(funcs ,x,y,sig ,ndata ,ma ,a);
82
83 covar=zeros(ma ,ma);
84 unitmat=eye(ma ,ma);
85 for i=1: ma
86 inhom=unitmat (:,i);
87 augmented =[ normal inhom];
88 cr=rref ( augmented );
89 covar(:,i)=cr(:,ma +1);
90 end
91
92 diag_covar =diag (covar);
93 sd_a =sqrt ( diag_covar )’;
94 return
95 else
96 aalt =a;
97 end
98 else
99 alamda=alamda *5;
100 end
101 end
1 function [alpha ,beta ,chisq] = ...
2 mrqcof(funcs ,x,y,sig ,ndata ,ma ,a);
3
4 alpha=zeros(ma ,ma);
5 beta =zeros(1,ma);
6 g=zeros(1, ndata);
7 g=1.0./ sig .^2;
8
9 chisq=0;
10
11 [ymod ,dyda ] = funcs(x,a,ma ,ndata);
12
47
Anhang A. Matlab-Quellcode
13 dy=y-ymod ;
14 for k=1: ndata
15 for i=1: ma
16 for j=1: ma
17 alpha(i,j)= alpha(i,j)+ g(k)* dyda (i,k)* dyda(j,k);
18 end
19 beta (i)=beta (i)+g(k)*dy(k)* dyda (i,k);
20 end
21 chisq=chisq+g(k)*dy(k)^2;
22 end
Listing A.5: Regression mit nichtlinearen Parametern – Modellfunktion:Regressionsfunktion
1 function [z,dzda] = expFunc (x,a,ma ,ndata)
2
3 if (nargin <=3)
4 ndata=length(x);
5 end
6 if (nargin <=2)
7 ma=length(a);
8 end
9
10 dzda =zeros(ma ,ndata);
11 z=zeros(1, ndata);
12
13 z=expFunc2 (a,x);
14
15 if (nargin >2)
16 fac1 =x.^( ones (1, ndata)*a(2));
17 fac2 =exp(-a(3)*x);
18
19 dzda (1 ,:)= fac1 .* fac2 ;
20 dzda (2 ,:)=a(1).* fac1 .* log(x).* fac2 ;
21 dzda (3,:)=-a(1)*x.^( ones (1, ndata )*(a(2)+1)).* fac2;
22 end
Listing A.6: Regression mit nichtlinearen Parametern – Modellfunktion:Hauptprogramm
1 % Nonlinear regression
2 % for y(x)=a*x^b*exp(-c*x).
3 %
4 % Author: Knap Michael
5
6 % reset
7 clear all;
8 close all
9 clc;
48
Anhang A. Matlab-Quellcode
10
11 % create function and mock data
12 x=1:0.02:4 -0.02;
13 xPrediction =0:0.02:5 -0.02;
14 n=length(x);
15 a=5; b=3; c=2;
16 yExact=a*x.^( ones (1,n)*b).*exp(-c*x);
17
18 sigma =0.2;
19 y=yExact+normrnd (0, sigma*ones (1, length(yExact )));
20
21 % guessed values for parameters
22 p0 =[1 ,1 ,1];
23 nParameters =length(p0);
24
25 % apply Levenberg -Marquardt
26 sig=ones (size (y));
27 eps =10^ -6;
28 [p,sd_p ,var ,CResidual ]= mrqmin(@expFunc ,x,y,sig ,n ,...
29 nParameters ,p0 ,eps );
30 yModel = expFunc(x,p);
31 % once again for correct covariance matrix
32 sSquare =1/(n- nParameters )* sum((y-yModel ).^2);
33 sig=sig*sqrt (sSquare );
34 [p,sd_p ,var ,CResidual ,cT , chiSqFit ]= mrqmin(@expFunc ,...
35 x,y,sig ,n,nParameters ,p0 ,eps);
36
37 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
38 %%% Jackknife
39 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
40 p_i=zeros(n,nParameters );
41 m=n;
42 xPrediction =0:0.02:5 -0.02;
43 yConfidence_i =zeros(n, length(xPrediction ));
44 for k=1:m
45 xTmp =[x(1:k-1),x(k+1:n)]; yTmp =[y(1:k-1),y(k+1:n)];
46 [p_i(k ,:)]= mrqmin(@expFunc ,xTmp ,yTmp ,sig (1:n-1) ,...
47 n-1, nParameters ,p0 ,eps);
48 theta_i(k)=p_i(k ,3)/ p_i(k ,2);
49 yConfidence_i (k ,:)= expFunc (xPrediction ,p_i(k ,:));
50 end
51 pDot =mean (p_i );
52 biasJack =(n -1)*(pDot -p);
53 pJack=pDot ;
54 pAdjustedJack =p-biasJack ;
55
56 varJack=zeros(1, nParameters ); for k=1: nParameters
57 varJack(k)=(n -1)/n*sum(( p_i(:,k)-p(k)).^2);
58 end
59
60 mseJack=varJack+biasJack .^2
61
62 yModelAdjustedJack = expFunc (x,pAdjustedJack );
63
49
Anhang A. Matlab-Quellcode
64 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
65 %%% Bootstrap
66 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
67
68 B=999; % estimated
69 % choose randomly (uniform distributed ) from x and y
70 index=ceil (n*rand (B,n));
71 xRand=x(index);
72 yRand=y(index);
73
74 % compute parameters
75 pPB_i=zeros(B, nParameters );
76 yConfidencePB_i =zeros(B,length(xPrediction ));
77 exitflags =zeros(B ,1);
78 for k=1:B
79 [pPB_i(k ,:)]= mrqmin(@expFunc ,xRand(k,:), yRand(k ,:) ,...
80 sig ,n,nParameters ,p0 ,eps);
81 thetaPB_i (k)= pPB_i(k ,3)/ pPB_i(k,2);
82 yConfidencePB_i (k ,:)= expFunc(xPrediction ,pPB_i(k ,:));
83 end
84
85 % compute bias
86 pPBHat=mean (pPB_i);
87 biasPB=pPBHat -p;
88 pAdjustedPB =p-biasPB;
89
90 % compute variance
91 varPB=zeros(1, nParameters );
92 for k=1: nParameters
93 varPB(k)=1/B*sum(( pPB_i(:,k)-p(k)).^2);
94 end
95
96 msePB=varPB+biasPB .^2
97
98 yModelAdjustedPB =expFunc (x,pAdjustedPB );
99
100 % prediction interval - Bootstrap
101 yConfidencePB_i =sort ( yConfidencePB_i );
102 yConfidencePBTop = yConfidencePB_i ((B+1)*0.025 ,:);
103 yConfidencePBBottom = yConfidencePB_i ((B+1)*0.975 ,:);
104 y_area=[ yConfidencePBBottom ,fliplr(yConfidencePBTop )];
105 x_area=[ xPrediction ,fliplr(xPrediction )];
106
107 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
108 %%% errorpropagation
109 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
110 disp (’Fehlerfortpflanzung ’);
111 disp (’theta=c/b’);
112 disp (’sigmaTheta =sqrt ((-c/b^2* sigmaB )^2+(1/b*sigmaC )^2) ’);
113 disp (’ ’);
114 theta=p(3)/ p(2);
115 sigmaTheta =sqrt ((-p(3)/ p(2)^2* sqrt (CResidual (2 ,2)))^2...
116 +(1/ b*sqrt (CResidual (3 ,3)))^2);
117
50
Anhang A. Matlab-Quellcode
118 % with covariances
119 % Jacobian
120 J=[0,-p(3)/ p(2)^2 ,1/p(2)];
121 varThetaCo =J*CResidual *J’;
122
123 % Jackknife
124 thetaDot =mean (theta_i );
125 biasThetaJack =(n -1)*( thetaDot -theta);
126 thetaAdjustedJack =theta -biasThetaJack ;
127 varThetaJack =(n -1)/n*sum(( theta_i -theta ).^2);
128
129 mseThetaJack = varThetaJack +biasThetaJack .^2;
130
131 % Bootstrap
132 thetaPBHat =mean (thetaPB_i );
133
134 biasThetaPB =thetaPBHat -theta;
135 thetaAdjustedPB =theta - biasThetaPB ;
136 varThetaPB =1/B*sum(( thetaPB_i -theta ).^2);
137
138 mseThetaPB =varThetaPB + biasThetaPB .^2;
Listing A.7: Regression mit nichtlinearen Parametern – Teleskop: Regres-sionsfunktion
1 function [z,dzda] = magnification (x,a, ma , ndata)
2
3 if (nargin <=3)
4 ndata=length(x);
5 end
6 if (nargin <=2)
7 ma=length(a);
8 end
9
10 dzda =zeros(ma ,ndata);
11 dzda (1 ,:)=x./(x-a(2));
12 dzda (2 ,:)=a(1)*x./(x-a (2)).^2;
13
14 z=zeros(1, ndata);
15 z=a(1)* x./(x-a(2));
Listing A.8: Regression mit nichtlinearen Parametern – Teleskop: Haupt-programm
1 % Nonlinear regression for telescope data .
2
3 % Author: Michael Knap
4
5 close all;
51
Anhang A. Matlab-Quellcode
6 clear all;
7 clc;
8
9 % measured data
10 xData =[2.0 ,2.5 ,3.0 ,3.5 ,4.0 ,4.5 ,5.0 ,5.5 ,6.0 ,6.5 ,7.0 ,7.5 ,8.0 ,...
11 8.5 ,9.0 ,9.5 ,10.0];
12 yData =[29.57 ,26.96 ,25.00 ,24.06 ,24.12 ,22.92 ,22.43 ,22.40 ,...
13 22.26 ,21.80 ,21.84 ,22.18 ,21.46 ,21.25 ,21.52 ,21.24 ,20.83];
14
15 v0 =25; deltav0 =0.4;
16 M0 =33; deltaM =0.05;
17 m0 =3; deltam =0.05;
18 l0 =3; deltal =0.01;
19
20 M=[31 ,34 ,33 ,40 ,38 ,45 ,47 ,45 ,42 ,41 ,47 ,31 ,41 ,44 ,46 ,41 ,44];
21 m=[5 ,4 ,3 ,3 ,2.5 ,2.5 ,2.3 ,2 ,1.7 ,1.5 ,1.6 ,1.0 ,1.2 ,1.2 ,1.2 ,...
22 1.0 ,1.0];
23
24 n=length(xData);
25
26 % guessed values for parameters
27 p0 =[1 ,1];
28 nParameters =length(p0);
29 xModel = [2:0.01:15];
30
31 % errorpropagation
32 fac1 =v0*M0/( m0*l0)*(m.* xData )./M;
33 dvdv0=fac1 /v0;
34 dvdM0=fac1 /M0;
35 dvdm0=-fac1 /m0;
36 dvdl0=-fac1 /l0;
37
38 dvdm =fac1 ./m;
39 dvdl =fac1 ./ xData;
40 dvdM =-fac1 ./M;
41
42 sig=sqrt (( dvdv0*deltav0 ).^2+( dvdM0*deltaM ).^2+...
43 (dvdm0*deltam ).^2+( dvdl0*deltal ).^2+...
44 (dvdm*deltam )^2+( dvdM *deltaM ).^2+...
45 (dvdl*deltal ).^2);
46 [p,sd_p ,var ,CResidual ,cT]= mrqmin(@magnification ,...
47 xData ,yData ,sig ,n,nParameters ,p0 ,eps);
48 yModel = magnification (xModel ,p);
49
50 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
51 %%% Jackknife
52 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
53 p_i=zeros(n ,2);
54
55 % calculate parameters
56 lF=xModel;
57 v_i=zeros(n,length(lF));
58 for k=1:n
59 xDataTmp =[ xData(1:k-1), xData(k+1:n)];
52
Anhang A. Matlab-Quellcode
60 yDataTmp =[ yData(1:k-1), yData(k+1:n)];
61 sigTmp=[ sig (1:k-1),sig(k+1:n)];
62
63 [p_i(k ,:)]= mrqmin(@magnification ,xDataTmp ,yDataTmp ,...
64 sigTmp ,n-1, nParameters ,p0 ,eps);
65 v_i(k,:)= p_i(k ,1)* lF ./(lF -p_i(k,2));
66 end
67 pDot =mean (p_i );
68 biasJack =(n -1)*(pDot -p);
69 pAdjustedJack =p-biasJack ;
70
71 varJack (1)=(n -1)/n*sum (( p_i(:,1)-p(1)).^2);
72 varJack (2)=(n -1)/n*sum (( p_i(:,2)-p(2)).^2);
73
74 mseJack=varJack+biasJack .^2;
75
76 yModelAdjustedJack = magnification (xModel ,pAdjustedJack );
77
78 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
79 %%% Bootstrap
80 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
81
82 B=999;
83 % choose randomly (uniform distributed ) from X
84 index=ceil (n*rand (B,n));
85 xRand=xData(index);
86 yRand=yData(index);
87 sigRand=sig(index);
88
89 % compute parameters
90 p_i=zeros(B ,2);
91 yConfidencePB_i =zeros(B,length(lF ));
92 for k=1:B
93 [p_i(k ,:)]= mrqmin(@magnification ,xRand(k,:), yRand(k ,:) ,...
94 sigRand(k,:),n,nParameters ,p0 ,eps);
95 yConfidencePB_i (k ,:)= magnification (lF ,p_i(k ,:));
96 end
97
98 pHat =mean (p_i );
99
100 biasPB=pHat -p;
101 pAdjustedPB =p-biasPB;
102
103 varPB (1)=1/B*sum(( p_i (:,1)-p(1)).^2);
104 varPB (2)=1/B*sum(( p_i (:,2)-p(2)).^2);
105
106 msePB=varPB+biasPB .^2;
107
108 yModelAdjustedPB = magnification (xModel , pAdjustedPB );
109
110 % prediction interval - Bootstrap
111 yConfidencePB_i =sort ( yConfidencePB_i );
112 yConfidencePBTop = yConfidencePB_i ((B+1)*0.025 ,:);
113 yConfidencePBBottom = yConfidencePB_i ((B+1)*0.975 ,:);
53
Anhang A. Matlab-Quellcode
114 y_area=[ yConfidencePBBottom ,fliplr(yConfidencePBTop )];
115 x_area=[lF ,fliplr(lF)];
A.4 Integrierte Autokorrelationszeit
Listing A.9: Erzeugung von korrelierten Daten
1 % Correlated data
2 % Generate correlated data .
3 %
4 % Author: Knap Michael
5
6 % reset
7 close all
8 clear all
9 clc
10
11 % create independent (uncorrelated ) data
12 my =0;
13 sigma=1;
14 n=100000;
15 nDisplay =1000;
16
17 x=normrnd(my*ones (1,n),sigma*ones (1,n));
18 figure (1)
19 plot ([1: nDisplay ],x(1: nDisplay ),’.’);
20
21 % calculate correlated data
22 tau =[10 ,50];
23 rho=exp (-1./tau );
24
25 for l=1: length(tau)
26 e(l,1)= x(1);
27 for k=2:n
28 fac1 =rho(l)^(k -1)* e(l ,1);
29 fac2 =sqrt (1-rho(l)^2);
30 e(l,k)= fac1+fac2 *sum(rho(l).^(k-[2: k]).* x(2:k));
31 end
32 figure(l+1)
33 plot ([1: nDisplay ],e(l ,1: nDisplay ),’.’);
34 disp ([’finished correlation ’, num2str(l),’ with tau=’, ...
35 num2str(tau(l))]);
36 end
37
38 % write to file
39 save (’data2.dat ’,’x’,’-ascii’);
40 save (’data2.dat ’,’e’,’-ascii’,’-append’);
54
Anhang A. Matlab-Quellcode
Listing A.10: Berechnung der Autokorrelationsfunktion und der integrier-ten Autokorrelationszeit
1 % Correlated data
2 % Calculate autocorrelationfunction and integrated
3 % autocorrelation time .
4 %
5 % Author: Knap Michael
6
7 % reset
8 clear all
9 close all
10 clc;
11
12 % read data
13 my =0;
14 sigma =1;
15 tau_exp =[10 ,50] ’;
16 rho=exp (-1./ tau_exp );
17
18 disp (’Reading data ... ’);
19 tic
20 load (’data .dat ’)
21 x=data (1 ,:);
22 e=data (2: length(tau_exp )+1 ,:);
23 n=length(x);
24 toc
25 disp (’ ’);
26
27 % calculate autocorrelation -coefficient
28 disp (’Autokorrelationsfunktion (theoretisch )’);
29 tic
30 kMaxAtA =100;
31 k=[0:0.1: kMaxAtA ];
32 ATheor=( repmat(rho ,1, length(k))).^( repmat(k,...
33 length(tau_exp ) ,1));
34 toc
35
36 for m=1: length(tau_exp)
37 disp ([’Autokorrelationsfunktion fur Tau_exp=’ ,...
38 num2str(tau_exp (m))]);
39 tic
40 aDenominator =1/n*sum(e(m ,:).^2) -(1/n*sum(e(m ,:)))^2;
41 for l=[1:1: kMaxAtA]
42 nDash=n-l;
43 aEnumerator (l)=1/ nDash*sum(e(m,1: nDash).*e(m,1+l:n))...
44 -1/nDash*sum(e(m ,1: nDash ))*1/ nDash*sum(e(m,1+l:n));
45 end
46 A(m,:)= aEnumerator / aDenominator ;
47 toc
48 end
49 disp (’ ’);
50
51 % calculate integrated autocorrelation time
52 disp (’integrierte Autokorrelationszeit (theoretisch )’);
55
Anhang A. Matlab-Quellcode
53 tic
54 k_max =[0:0.1: kMaxAtA ];
55 tau_intTheor =tau_exp .*(1+1./(12.* tau_exp .^2));
56 fac1 =2* tau_exp ./(2* tau_exp +1);
57 fac2 =exp(-repmat(k_max , length(tau_exp ) ,1)./...
58 repmat(tau_exp ,1, length(k_max )));
59 fac3 =1- repmat(fac1 ,1, length(k_max )).* fac2 ;
60 tau_intTheorOfk =repmat(tau_intTheor ,1, length(k_max )).* fac3 ;
61 toc
62
63 disp (’integrierte Autokorrelationszeit’);
64 tic
65 tau_intOfk =cumsum(A ,2);
66 toc
67 disp (’ ’);
68
69 % calculate tau_int with restriction : k_max >=6* tau_int(kmax )
70 for m=1: length(tau_exp)
71 for l=[1:1: kMaxAtA]
72 if (l >=6* tau_intOfk (m,l))
73 sig= tau_intOfk (m,l)*sqrt (2*(2*l+1)/ n);
74 disp ([ ’Autokorrelationszeit mit Bedingung ...
75 k_max >=6* tau_int(kmax ) fur Tau_exp=’ ,...
76 num2str(tau_exp (m)),’: ’, num2str (...
77 tau_intOfk (m,l)),’(+/-)’,num2str(sig )]);
78 break;
79
80 end
81 end
82 end
Listing A.11: Berechnung der integrierten Autokorrelationszeit mit Jack-knife
1 % Correlated data
2 % Calculate integrated correlation time with Jackknife .
3 %
4 % Author: Knap Michael
5
6 % reset
7 clear all
8 close all
9 clc;
10
11 % read data
12 my =0;
13 sigma=1;
14 tau_exp =10;
15 rho=exp (-1./ tau_exp );
16 kMax =20;
17
18 disp (’Reading data ... ’);
56
Anhang A. Matlab-Quellcode
19 tic
20 load (’data .dat ’)
21 % read only tau =10
22 e=data (2 ,:);
23 n=length(e);
24 toc
25 disp (’ ’);
26
27 % calculate autocorrelation function
28 aDenominator =1/n*sum(e.^2) -(1/ n*sum(e))^2;
29 for l=[1:1: kMax]
30 nDash=n-l;
31 aEnumerator (l)=1/ nDash*sum(e(1: nDash ).*e(1+l:n))...
32 -1/nDash*sum(e(1: nDash ))*1/ nDash*sum(e(1+l:n));
33
34 end
35 A=aEnumerator / aDenominator ;
36
37 % linear regression with log(A)
38 logA =log(A);
39
40 x=([1:1: kMax ]);
41 % linear regression
42 Cxx =1/ kMax *sum(x.^2) -(1/ kMax *sum(x))^2;
43 Cyy =1/ kMax *sum(logA .^2) -(1/ kMax *sum(logA ))^2;
44 Cxy =1/ kMax *sum(x.* logA )-(1/ kMax *sum(x)*1/ kMax *sum(logA ));
45
46 logAQuer =mean (logA );
47 xQuer=mean (x);
48
49 p(2)= Cxy/Cxx;
50 p(1)= logAQuer -xQuer*p(2);
51
52 % calculate areas
53 area = -1/2+ exp(p(1))*1/(1- exp(p(2)));
54
55 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
56 %%% Jackknife
57 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
58 numberOfBlocks =20;
59
60 for m=1: numberOfBlocks
61 tic
62 nBlock=n/ numberOfBlocks ;
63 % save original data
64 eTemp=e;
65
66 e(nBlock *(m -1)+1: nBlock*m)=0;
67 eFullData =[e(1: nBlock *(m-1)),e(nBlock*m+1: end )];
68
69 % calculate autocorrelation function
70 aDenominator =1/(n-nBlock )* sum(e.^2) -...
71 (1/(n-nBlock )* sum(e))^2;
72 for l=[1:1: kMax]
57
Anhang A. Matlab-Quellcode
73 nDash=n-l;
74 if (m==1)
75 aEnumerator (l)=1/(nDash -nBlock)*sum(e(1: nDash ).*...
76 e(1+l:n)) -1/( nDash -nBlock)* sum(e(1: nDash ))*1/...
77 (nDash -nBlock+l)* sum(e(1+l:n));
78 end
79 if (m== numberOfBlocks )
80 aEnumerator (l)=1/(nDash -nBlock)*sum(e(1: nDash ).*...
81 e(1+l:n)) -1/( nDash -nBlock+l)* sum(e(1: nDash ))*1/...
82 (nDash -nBlock )* sum(e(1+l:n));
83 else
84 aEnumerator (l)=1/(nDash -nBlock -l)* sum(e(1: nDash ...
85 ).*e(1+l:n)) -1/( nDash -nBlock)*sum(e(1: nDash ))*1/...
86 (nDash -nBlock )* sum(e(1+l:n));
87 end
88 end
89 AJack(m ,:)= aEnumerator /aDenominator ;
90
91 % restore original data
92 e=eTemp;
93 disp ([’Autokorrelationsfunktion fur Block ’,num2str(m ,...
94 ’%2.2 d’),’ - Dauer: ’, num2str(toc , ’%5.3 f’)]);
95 end
96
97 % linear regression with log(A)
98 logAJack =log(AJack);
99
100 x=([1:1: kMax ]);
101 for m=1: numberOfBlocks
102 % linear regression
103 Cxx =1/ kMax*sum(x.^2) -(1/ kMax *sum(x))^2;
104 Cyy =1/ kMax*sum(logAJack (m ,:).^2) -...
105 (1/ kMax*sum(logAJack (m ,:)))^2;
106 Cxy =1/ kMax*sum(x.* logAJack (m ,:)) -...
107 (1/ kMax*sum(x)*1/ kMax *sum(logAJack (m ,:)));
108
109 logAJackQuer =mean (logAJack (m ,:));
110 xQuer=mean (x);
111
112 pJack(m ,2)= Cxy/Cxx;
113 pJack(m ,1)= logAJackQuer -xQuer*pJack(m,2);
114
115 % calculate area
116 areaJack (m)= -1/2+ exp(pJack(m ,1))*1/(1- exp(pJack(m ,2)));
117 end
118
119 tauint=mean (areaJack );
120 biasTauintJack =( numberOfBlocks -1)*( tauint -area );
121 tauintAdjustedJack =tauint -biasTauintJack ;
122 % calculate variance
123 varTauintJack =( numberOfBlocks -1)/( numberOfBlocks )*...
124 sum(( areaJack -tauint ).^2);
125 % calculate the mse
126 mseTauintJack = varTauintJack +biasTauintJack ^2;
58
Literaturverzeichnis
[1] Bevington, P. R. und D. K. Robinson: Data Reduction and ErrorAnalysis for the Physical Sciences. McGraw-Hill, New York, 3 Aufl.,2003.
[2] Davison, A. C. und D. V. Hinkley: Bootstrap Methods and TheirApplication, Bd. 1 d. Reihe Cambridge Series in Statistical and Proba-bilistic Mathematics. Cambridge University Press, Cambridge, 1997.
[3] Evertz, H. G.: Computersimulationen (515.819), 2007.
[4] Hinkley, D. V.: Jackknifing in Unbalanced Situations. Technometrics,19(3):285–292, August 1977.
[5] Janke, W.: Quantum Simulations of Complex Many-Body Systems:From Theory to Algorithms, Bd. 10 d. Reihe NIC Series, Kap. Statisti-cal Analysis of Simulations: Data Correlations and Error Estimation,S. 423–445. John von Neumann Institute for Computing, Julich, 2002.
[6] Papoulis, A.: Probability, Random Variables, and Stochastic Proces-ses. McGraw-Hill, Kogakusha, International Student Edition Aufl.,1965.
[7] Shao, J. und D. Tu: The Jackknife and Bootstrap. Springer-Verlag,New York, 1995.
[8] Sormann, H.: Numerische Methoden in der Physik (515.421), 2006.
[9] Wu, C. F. J.: Jackknife, Bootstrap and Other Resampling Methods inRegression Analysis. The Annals of Statistics, 14(4):1261–1295, De-cember 1986.
[10] Yang, M. C. K. und D. H. Robinson: Understanding and LearningStatistics by Computer , Bd. 4 d. Reihe Series in Computer Science.World Scientific Publishing Co Pte Ltd., Singapore, 1986.
59
60
Messbox zur Druckkontrolle
— Druckgroße kontrollieren! —
Breite = 100 mmHohe = 50 mm
— Diese Seite nach dem Druck entfernen! —
61
top related