Mainz, 25. Juni 2015

Statistics, Data Analysis, and Simulation – SS 2015 08.128.730 Statistik, Datenanalyse und Simulation

Dr. Michael O. Distler

Dr. Michael O. Distler

Statistics, Data Analysis, and Simulation – SS 2015

1 / 28

8. Einführung in die Bayes-Statistik Wiederholung von: 1.1 Theory of probability Probability theory, mathematics:

−→ Kolmogorov axioms Classical interpretation, frequentist probability: Pragmatical definition of probability: n N→∞ N

p(E) = lim

n(E) = number of events E N = number of trials (experiments) Experiments have to be repeatable (in principle). Disadvantage: Strictly speaking one cannot make statements on the probability of any true value. Only upper and lower limits are possible given a certain confidence level. Dr. Michael O. Distler

Statistics, Data Analysis, and Simulation – SS 2015

2 / 28

Theory of probability Probability theory, mathematics Classical interpretation, frequentist probability Bayesian statistics, subjective probability: Prior subjective assumptions enter into the calculation of probabilities of a hypotheses H. p(H) = degree of belief that H is true Metaphorically speaking: Probabilities are the ratio of the (maximum) wager and the anticipated prize in a bet.

Dr. Michael O. Distler

Statistics, Data Analysis, and Simulation – SS 2015

3 / 28

Theory of probability Probability theory, mathematics Classical interpretation, frequentist probability Bayesian statistics, subjective probability:

Disadvantage: Prior hypotheses influence the probability. Advantages for rare and one-time events, like noisy signals or catastrophe modeling.

Dr. Michael O. Distler

Statistics, Data Analysis, and Simulation – SS 2015

3 / 28

Theory of probability Probability theory, mathematics Classical interpretation, frequentist probability Bayesian statistics, subjective probability:

Disadvantage: Prior hypotheses influence the probability. Advantages for rare and one-time events, like noisy signals or catastrophe modeling. In this lecture we focused on the classical statistics, e.g. error estimates have to be understood as confidence regions.

Dr. Michael O. Distler

Statistics, Data Analysis, and Simulation – SS 2015

3 / 28

Theory of probability Probability theory, mathematics Classical interpretation, frequentist probability Bayesian statistics, subjective probability:

Disadvantage: Prior hypotheses influence the probability. Advantages for rare and one-time events, like noisy signals or catastrophe modeling. In this chapter we will discuss the Bayesian statistics. Bayesian inference (e.g. mean, variance) is derived solely from posterior distributions (probabilities).

Dr. Michael O. Distler

Statistics, Data Analysis, and Simulation – SS 2015

3 / 28

Bayes–Theorem

Aus der Gleichung p(A und B) = p(A) · p(B|A) = p(B) · p(A|B) erhält man das Bayes-Theorem: p(A|B) = p(B|A) ·

Dr. Michael O. Distler

p(A) p(B)

Statistics, Data Analysis, and Simulation – SS 2015

4 / 28

Bayes–Theorem für diskrete Ereignisse p(A|B) = p(B|A) ·

p(A) p(B)

Beispiel: In einem Experiment soll der leptonische Zerfall der ˇ K0 -Mesonen studiert werden. Es ist geplant, einen CerenkovDetektor zu verwenden, um die Leptonen nachzuweisen. Dazu muss untersucht werden, ob ein Detektor ausreicht, um die leptonischen Ereignisse von dem kleinen Untergrund abzutrennen, der ebenfalls den Detektor auslösen kann. p(B) ≡ Wahrscheinlichkeit, dass ein Ereignis den ˇ Cerenkov-Detektor auslöst. p(A) ≡ Wahrscheinlichkeit, dass sich ein echter leptonischer Zerfall ereignet. p(B|A) ≡ Wahrscheinlichkeit, dass ein echtes leptonisches ˇ Ereignis den Cerenkov-Detektor auslöst p(A|B) ≡ Wahrscheinlichkeit, dass es sich bei einem Ereignis um einen echten leptonischen Zerfall handelt, unter der ˇ Voraussetzung, dass der Cerenkov-Detektor auslöst. Dr. Michael O. Distler

Statistics, Data Analysis, and Simulation – SS 2015

5 / 28

Bayes–Theorem für diskrete Ereignisse

p(A|B) = p(B|A) · p(B)

p(A) p(B)

≡

Wahrscheinlichkeit, dass ein Ereignis den ˇ Cerenkov-Detektor auslöst. p(A) ≡ Wahrscheinlichkeit, dass sich ein echter leptonischer Zerfall ereignet. p(B|A) ≡ Wahrscheinlichkeit, dass ein echtes leptonisches ˇ Ereignis den Cerenkov-Detektor auslöst p(A|B) ≡ Wahrscheinlichkeit, dass es sich bei einem Ereignis um einen echten leptonischen Zerfall handelt, unter der ˇ Voraussetzung, dass der Cerenkov-Detektor auslöst. p(B) kann gemessen werden. p(A) ergibt sich aus früheren Messungen bzw. Theorie. p(B|A) wird aus einer Simulation bestimmt. ,→ p(A|B)

Dr. Michael O. Distler

Statistics, Data Analysis, and Simulation – SS 2015

6 / 28

Bayes-Theorem für Bayesianer Wenn es sich bei A und B nicht um Ereignisklassen sondern um Hypothesen handelt, dann wird der Unterschied zwischen den beiden Statistik-Schulen offensichtlich. Als frequentist kann man einer Hypothese keine Wahrscheinlichkeit zuweisen. Der bayesian interpretiert p(H) als Grad des Vertrauens in die Hypothese. p(H|E) = p(H)

≡

p(H|E) p(E|H) p(E)

≡ ≡ ≡

p(H) · p(E|H) p(E)

prior Wahrscheinlichkeit Wissen (Grad des Vertrauens) vor der Datennahme posterior Wahrscheinlichkeit likelihood Normalisierungsfaktor

Das Ergebnis (Erwartungswert, Varianz, . . . ) einer Bayes-Analyse wird allein dem posterior entnommen. Dr. Michael O. Distler

Statistics, Data Analysis, and Simulation – SS 2015

7 / 28

Cookie-Problem Stellen Sie sich 2 gefüllte Keksdosen vor. Dose 1 enthält 10 Chocolate Chip Cookies und 30 Plain Cookies. Bei Dose 2 ist des Verhältnis 20/20. Unser Freund Fred wählt zunächst zufällig eine Keksdose aus und entnimmt dann zufällig einen Cookie. Es ist ein Plain Cookie. Mit welcher Wahrscheinlichkeit stammt er aus Dose 1? Hypothesen: H1 : der Keks stammt aus Dose 1. H2 : der Keks stammt aus Dose 2. Prior: p(H1 ) = p(H2 ) = 1/2 Ereignis: E: der Keks ist ein Plain Cookie. Likelihood: p(E|H1 ) = 3/4 p(E|H2 ) = 1/2 Bayes-Theorem: p(H1 |E) =

3 p(H1 ) × p(E|H1 ) = p(H1 ) · p(E|H1 ) + p(H2 ) · p(E|H2 ) 5

Dr. Michael O. Distler

Statistics, Data Analysis, and Simulation – SS 2015

8 / 28

Bayessche Inferenz∗ für Binomialexperimente Häufig hat man es mit großen Populationen zu tun, von der ein Anteil p ein gewisse Eigenschaft aufweist. Beispiel: Für die Bevölkerung einer Stadt könnte die Eigenschaft lauten „plant Kandidat A bei der Bürgermeisterwahl zu wählen“. Wir zählen die Anzahl von “Erfolgen” in n unabhängigen Versuchen, wobei jeder Versuch nur zwei mögliche Ergebnisse haben kann: Erfolg oder Miserfolg. Erfolg bedeutet, dass bei dem i-ten Versuch die geforderte Eigenschaft auftrat. Die Anzahl von Erfolgen in n Versuchen, die Erfolgswahrscheinlichkeit bei jedem einzelnen Versuch sei p, ist binomialverteilt. Die bedingte Wahrscheinlichkeit ist:   n f (k |p) = pk (1 − p)n−k k ∈ [1, n] k ∗

Dr. Michael O. Distler

Inferenz = Schlussfolgerung

Statistics, Data Analysis, and Simulation – SS 2015

9 / 28

Betaverteilung Wahrscheinlichkeitsdichte der Betaverteilung f (x; a, b) =

1 x a−1 (1 − x)b−1 B(a, b)

x ∈ [0, 1]

mit der Eulerschen Betafunktion Z 1 Γ(a)Γ(b) B(a, b) = = u a−1 (1 − u)b−1 du Γ(a + b) 0 Extremum: xextrem

  b − 1 −1 a−1 = 1+ = a−1 a+b−2

Erwartungswert und Varianz: E(X ) =

a a+b

Var(X ) =

Dr. Michael O. Distler

a·b (a + b + 1)(a + b)2

Statistics, Data Analysis, and Simulation – SS 2015

10 / 28

Betaverteilung beta(.5,.5)

0

0.5

beta(.5,1.)

1

0

beta(1.,.5)

0

0.5

0.5

1

0

0.5

0

0.5

1

0

0.5

1

0

0

0.5

Dr. Michael O. Distler

1

0

0.5

1

0

0.5

1

0

0

0.5

1

0.5

1

beta(2.,3.)

1

0

beta(3.,2.)

1

0.5 beta(1.,3.)

beta(2.,2.)

beta(3.,1.)

1

0.5

beta(.5,3.)

beta(1.,2.)

beta(2.,1.)

beta(3.,.5)

0

1

beta(1.,1.)

beta(2.,.5)

0

0.5

beta(.5,2.)

0.5

1

beta(3.,3.)

1

0

0.5

1

Statistics, Data Analysis, and Simulation – SS 2015

11 / 28

Referendum Beispiel: Es steht ein Referendum wegen eines Bauvorhabens in Ihrer Stadt aus. Da Sie im Bekanntenkreis das Thema bereits diskutiert haben, gehen Sie von einer knappen Entscheidung aus, wobei Sie sich sicher sind (C.L.: 95%), dass weder Befürworter noch Gegner mehr als 60% der Stimmen erreichen werden. Aus Ihrem Vorwissen konstruieren Sie den Prior: a E(X ) = = 0.5 ⇒ a = b a+b Var(X ) =

(a +

1 ab = (0.05)2 = 4(2a + 1) + b + 1)

b)2 (a

Nähert man also die Betaverteilung mit der Normalverteilung und setzt 95%c.l. ' 2σ ' ±10% so ergibt sich a = b ' 50 (die exakte Rechnung ergibt a = b = 47.2998). Dr. Michael O. Distler

Statistics, Data Analysis, and Simulation – SS 2015

12 / 28

Referendum (exakte Berechnung der Betaverteilung) zu lösen ist: Z

0.6

f (x; a, a) dx = 0.95 0.4

Mathematica: f [α_] := NIntegrate[PDF[BetaDistribution[α, α], x], {x, 0.4, 0.6}] FindRoot[f [α] == 0.95, {α, 50}]

Python: from scipy.stats import beta from scipy.optimize import newton def f(x): return (beta.cdf(0.6,x,x) -beta.cdf(0.4,x,x)-0.95) print(newton(f, 50))

Dr. Michael O. Distler

Statistics, Data Analysis, and Simulation – SS 2015

13 / 28

Referendum (Teil 2) In einer repräsentativen Umfrage haben sich von N = 1500 Betroffenen nur k = 720 Personen für das Bauvorhaben ausgesprochen. Ermitteln Sie die Wahrscheinlichkeit, dass im Referendum die Gegner eine Mehrheit erzielen. Die posterior Dichte g(x) ergibt sich aus: prior × likelihood Betaverteilung(x; a, b) × Binomialverteilung(x; N, k ) Γ(a+b) a−1 N! k N−k (1 − x)b−1 × k ! (N−k Γ(a)Γ(b) x )! x (1 − x) Im Falle eines Beta-Priors ergibt sich einfach: g(x) = Beta(x; a + k , b + N − k ) = Beta(x; 770, 830)

Dr. Michael O. Distler

Statistics, Data Analysis, and Simulation – SS 2015

14 / 28

Referendum (Teil 3) prior 12

posterior 45

Beta(x,50,50) Beta(x,30,70) Beta(x,0.5,0.5) Beta(x,1.,1.)

10

40 35

8

30 25

6

20

4

15 10

2 0

Beta(x,770,830) Beta(x,750,850) Beta(x,720.5,780.5) Beta(x,721,781)

5 0

0.2

0.4

0.6

0.8

1

0 0.42 0.44 0.46 0.48

0.5

0.52 0.54

Aus der roten Kurve (rechts) ermitteln wir unser Ergebnis: Z 0.5 Beta(x; 770, 830)dx = 0.933 0

Das Referendum wird also mit einer Wahrscheinlichkeit von 93.3% abgelehnt. Dr. Michael O. Distler

Statistics, Data Analysis, and Simulation – SS 2015

15 / 28

Referendum (Teil 4) 35

Beta(x,50,50) Beta(x,30,70) Beta(x,0.5,0.5) Beta(x,1.,1.) Beta(x,770,830) Beta(x,750,850) Beta(x,720.5,780.5) Beta(x,721,781)

30

25

20

15

10

5

0

0.2

0.3

0.4

0.5

0.6

0.7

Das Bild zeigt noch einmal deutlich, dass das Ergebnis nur schwach von der Wahl des Priors abhängt. Dr. Michael O. Distler

Statistics, Data Analysis, and Simulation – SS 2015

16 / 28

Vorsichtsmaßnahmen bei Verwendung eines konjugierten Priors z.B. Beta-Prior bei Binomialverteilungen 1 Plotten Sie Ihren Beta(a, b)-Prior. Passen Sie notfalls Mittelwert π0 und Varianz σ02 an, bis diese Ihren Vorstellungen entsprechen. 2 Berechnen Sie die äquivalente Stichprobengröße. Für den Fall, dass diese unrealistisch groß ist, vergrößern Sie die Varianz Ihres Priors und berechnen diesen neu. Für eine Binomialverteilung mit Trefferwahrscheinlichkeit π und Versuchsanzahl n ist die Varianz π(1 − π)/n. Dies setzen wir der Prior-Varianz gleich: π0 (1 − π0 ) ab = neq (a + b + 1)(a + b)2 Mit π0 =

a a+b

und (1 − π0 ) =

b a+b

ergibt sich

neq = a + b + 1 Dr. Michael O. Distler

Statistics, Data Analysis, and Simulation – SS 2015

17 / 28

Nichtinformative Priori-Dichten

Ein nichtinformativer (engl. uninformative or objective) Prior drückt eine vage bzw. unbestimmte Kenntnis der gesuchten Größe aus. Die einfachste und älteste Methode einen nichtinformativen Prior zu konstruieren stellt das Indifferenzprinzip dar. Demnach wird allen Möglichkeiten die gleiche Wahrscheinlichkeit zugewiesen. Dabei kann leicht ein uneigentlicher (engl. improper) Prior entstehen, d.h. der Prior ist nicht normiert und damit auch keine Wahrscheinlichkeitsdichte. Das stellt jedoch im allgemeinen kein Problem dar, da sich die Posterior-Dichte meist normieren lässt.

Dr. Michael O. Distler

Statistics, Data Analysis, and Simulation – SS 2015

18 / 28

Nichtinformative Priori-Dichten

Der “flache” Prior ist jedoch nicht wirklich “objektiv”, wovon man sich leicht überzeugen kann, wenn man eine (nicht-lineare) Variablentransformation durchführt. Nach der Transformation ist der flache Prior nicht mehr flach. Bessere Eigenschaften besitzt der Jeffreys Prior, der ebenfalls als nichtinformativer Prior bezeichnet wird. Eine Bayes-Analyse mit einem nichtinformativen Prior liefert meist ähnliche oder identische Ergebnisse wie die klassische Maximum Likelihood Methode.

Dr. Michael O. Distler

Statistics, Data Analysis, and Simulation – SS 2015

19 / 28

Bayessche Inferenz für Poisson

Die Poisson-Verteilung wird verwendet, um das Auftreten von seltenen Ereignissen zu zählen. Die Ereignisse treten zufällig in Zeit (oder Raum) auf, jedoch mit einer konstanten mittleren Rate. Die Poisson-Verteilung kann etwa verwendet werden, um die Anzahl der Unfälle auf einer Autobahn innerhalb eines Monats zu modellieren. Allerdings kann es nicht verwendet werden, um die Zahl der Todesopfer auf der Autobahn zu modellieren, da einige Unfälle mehrere Todesopfer aufweisen können.

Dr. Michael O. Distler

Statistics, Data Analysis, and Simulation – SS 2015

20 / 28

Die Gammaverteilung Wahrscheinlichkeitsdichte der Gammaverteilung  p  b x p−1 e−bx x > 0 f (x) = Γ(p) 0 x ≤0 Maximum (für p > 1): xmax =

p−1 b

Erwartungswert: E(X ) = Varianz: Var(X ) =

Dr. Michael O. Distler

p b p b2

Statistics, Data Analysis, and Simulation – SS 2015

21 / 28

1.0

Gammaverteilung

0.6

0.5, 0.5, 1, b 1, b 2, b 2, b

b=2 b=1 =2 =1 =2 =1

0.4

= = = = = =

0.0

0.2

Ɣ(p,b)

0.8

p p p p p p

0

1

2

3

4

5

x

Dr. Michael O. Distler

Statistics, Data Analysis, and Simulation – SS 2015

22 / 28

Bayes-Theorem für Poisson Parameter

Wir betrachten eine Stichprobe y1 , . . . , yn aus einer Poisson(µ) Verteilung. Die Proportionalitätsform des Bayes-Theorems lautet: posterior ∝ prior × likelihood g(µ|y1 , . . . , yn ) ∝ g(µ) × f (y1 , . . . , yn |µ) Durch Normierung erhalten wir die tatsächliche Posterior Dichte: g(µ|y1 , . . . , yn ) = R ∞ 0

Dr. Michael O. Distler

g(µ) × f (y1 , . . . , yn |µ) g(µ) × f (y1 , . . . , yn |µ) dµ

Statistics, Data Analysis, and Simulation – SS 2015

23 / 28

Likelihood für Poisson Parameter Die Likelihood für eine einmalige Ziehung von einer Poisson-Verteilung ist bekannt: f (y |µ) =

µy e−µ y!

Die Form wird dabei festgelegt durch f (y |µ) ∝ µy e−µ Für eine größere Stichprobe werden die ursprünglichen Likelihoods multipliziert: f (y1 , . . . , yn |µ) =

n Y i=1 P

∝ µ Dr. Michael O. Distler

f (yi |µ) yi

e−nµ

Statistics, Data Analysis, and Simulation – SS 2015

24 / 28

Gleichverteilte Prior Dichte

Wenn wir keine Information über µ haben, bevor wir die Daten betrachten, dann wäre ein gleichverteilter Prior eine mögliche Wahl: g(µ) = 1 für µ > 0 Dies ist ein uneigentlicher (improper) Prior! g(µ|y1 , . . . , yn ) ∝ g(µ) × f (y1 , . . . , yn |µ) ∝ 1×µ

P

yi

e−nµ

P Dies entspricht einer gamma(p, b) Verteilung mit p = y + 1 und b = n. Somit erhalten wir einen normierten Posterior, obwohl wir mit einem improper Prior gestartet waren.

Dr. Michael O. Distler

Statistics, Data Analysis, and Simulation – SS 2015

25 / 28

Jeffreys Prior für Poisson Ein Jeffreys Prior ist objektiv in dem Sinne, dass er invariant ist unter bestimmten Transformationen des Parameters. Der Jeffreys Prior für Poisson lautet: 1 g(µ) ∝ √ µ

für µ > 0

Dies ist ebenfalls ein uneigentlicher (improper) Prior! g(µ|y1 , . . . , yn ) ∝ g(µ) × f (y1 , . . . , yn |µ) P 1 ∝ √ × µ yi e−nµ µ P

∝ µ

yi −1/2

e−nµ

P Dies entspricht einer gamma(p, b) Verteilung mit p = y + 21 und b = n. Wiederum erhalten wir einen normierten Posterior, obwohl wir mit einem improper Prior gestartet waren. Dr. Michael O. Distler

Statistics, Data Analysis, and Simulation – SS 2015

26 / 28

Konjugierte Priors für Poisson Die Gammaverteilung bildet die Familie von konjugierten Priors für Poisson, d.h. sowohl Prior als auch Posterior stammen aus der gleichen Familie. Für eine Stichprobe y1 , . . . , yn aus einer Poissonverteilung und einer Prior gamma(p, b) ergibt sich der Posterior: X gamma(p0 , b0 ) mit p0 = p + y , b0 = b + n Der Prior lässt sich leicht aus einer Kenntnis von Mittelwert µ und Varianz s2 konstruieren. Aus µ=

p b

und s2 =

p b2

p=

µ2 s2

und b =

µ s2

folgt

Dr. Michael O. Distler

Statistics, Data Analysis, and Simulation – SS 2015

27 / 28

press any key

Dr. Michael O. Distler

Statistics, Data Analysis, and Simulation – SS 2015

28 / 28