python stats

Parameter Estimation

Javier Rico, IFAE

Master of Multidisciplinary Research in Experimental Sciences, UPF, 2018-2019

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Bibliography ★ G. Cowan, “Statistical Data Analysis”, 1998, Oxford

Science Publications

★ M. Tanabashi et al. (Particle Data Group), Phys. Rev. D 98, 030001 (2018).

2

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Basic concepts

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Description of the problem ★ Consider a random variable x, with probability density function (pdf) f(x)

✦ E.g.: height of human beings

✦ A set of independent observations of x: x1,…,xn is called a sample of size n or n independent and identically distributed values

✦ The whole dataset can be considered as one multi-dimensional random variable, with joint pdf:

fsample(x1,…,xn) = f(x1) f(x2)… f(xn)

(since observations are independent from each other and AND relation means multiplying probabilities)

★ The problem: we have n measurements of x, based on which we want to infer properties of the unknown pdf f(x): ✦ We may be interested in estimating general pdf parameters like the mean or the variance

✦ Also often we have a parametric expression for f(x;θ), which depends on unknown parameters θ, the goal is then to estimate those parameters

4

x

f(x)

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Example ★ We measure the height of 30

human beings, randomly selected

★ xi ={170.8, 177.5, 155.1, 180.4, 163.7, 146.2, 160.3, 179.0, 175.0, 152.8, 179.3, 168.7, 176.4, 175.0, 152.6, 175.3, 160.9, 193.2, 180.9, 170.7, 155.3, 170.8, 172.4, 151.9, 203.4, 175.9, 195.4, 153.3, 194.5, 149.6}

★ if:

★ 𝜇 and 𝜎 are unknown parameters

★ We want to build estimators, i.e.:

• •

5

d�PP

dE =

1

mdm⌧dm

dN�

dE (5)

J(⌦) = 1

4⇡

Z

⌦

Z

los ⇢(l) dl d⌦ (6)

⇢(r) = ⇢s

(r/rs)(1 + r/rs)2 (7)

L(s, b|n,m) = (s+ b)n

n! e�(s+b) (⌧b)

m

m! e�⌧b (8)

L(M(✓)|S,B) = Y

i2B PB(Ei)

Y

j2S PS(M(✓)|Ej) (9)

PB(E) =

Z 1

0 ⌧ �B(E

0)RB(E 0, E) dE0 (10)

PS(M(✓)|E) =

Z 1

0 �B(E

0)RB(E 0, E) dE0 + ✓0

Z 1

0 �G(E

0|✓)RG(E 0, E) dE0 (11)

LT (M(✓)|ST , BT ) = Y

i

Li(M(✓)|Si, Bi) (12)

� < 2 p

Non TEG/Tobs

Ae↵ TEG (13)

�(l, b) = �CR + �e+e� + �EG�� + �Gal��(l, b) + �EG-DM + �Gal-DM('(l, b)) (14)

�(l1, b1)� �(l2, b2) = �Gal��(l1, b1) + �Gal-DM('1)� �Gal��(l2, b2)� �Gal-DM('2) (15)

�Gal-DM('1) = �(l1, b1)� �(l2, b2)� �Gal��(l1, b1) + �Gal��(l2, b2)

1� �Gal-DM('2) �Gal-DM('1)

(16)

f(x) = 1p 2⇡�

e� 1 2(

x�µ � )2 (17)

2

d�PP

dE =

1

mdm⌧dm

dN�

dE (5)

J(⌦) = 1

4⇡

Z

⌦

Z

los ⇢(l) dl d⌦ (6)

⇢(r) = ⇢s

(r/rs)(1 + r/rs)2 (7)

L(s, b|n,m) = (s+ b)n

n! e�(s+b) (⌧b)

m

m! e�⌧b (8)

L(M(✓)|S,B) = Y

i2B PB(Ei)

Y

j2S PS(M(✓)|Ej) (9)

PB(E) =

Z 1

0 ⌧ �B(E

0)RB(E 0, E) dE0 (10)

PS(M(✓)|E) =

Z 1

0 �B(E

0)RB(E 0, E) dE0 + ✓0

Z 1

0 �G(E

0|✓)RG(E 0, E) dE0 (11)

LT (M(✓)|ST , BT ) = Y

i

Li(M(✓)|Si, Bi) (12)

� < 2 p Non TEG/Tobs

Ae↵ TEG (13)

�(l, b) = �CR + �e+e� + �EG�� + �Gal��(l, b) + �EG-DM + �Gal-DM('(l, b)) (14)

�(l1, b1)� �(l2, b2) = �Gal��(l1, b1) + �Gal-DM('1)� �Gal��(l2, b2)� �Gal-DM('2) (15)

�Gal-DM('1) = �(l1, b1)� �(l2, b2)� �Gal��(l1, b1) + �Gal��(l2, b2)

1� �Gal-DM('2) �Gal-DM('1)

(16)

f(x) = 1p 2⇡�

e� 1 2(

x�µ � )2 (17)

µ̂(x1, ..., xn) (18)

�̂(x1, ..., xn) (19)

2

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Estimators ★ Statistic = a function of the measurements (x1,…,xn) which contains no

unknown parameter (i.e. given the data we can compute the statistic, e.g. the sum of all measurements).

★ Estimator = a statistic used to estimate some property of a pdf like its mean, variance or others (e.g. 𝜇, 𝜎) noted in general as 𝜃 to distinguish it from the true (and forever unknown) parameter 𝜃.

★ An estimator 𝜃 is said to be consistent if it converges to 𝜃 in the limit of large n, i.e. for any 𝜖>0 arbitrarily small

which is a minimum requirement for a useful estimator. ✦ The limit n→∞ is referred to as the “large sample” or “asymptotic” limit

★ The procedure of estimating a parameter given the data (x1,…,xn) is called parameter fitting

6

^^^

^

d�PP

dE =

1

mdm⌧dm

dN�

dE (5)

J(⌦) = 1

4⇡

Z

⌦

Z

los ⇢(l) dl d⌦ (6)

⇢(r) = ⇢s

(r/rs)(1 + r/rs)2 (7)

L(s, b|n,m) = (s+ b)n

n! e�(s+b) (⌧b)

m

m! e�⌧b (8)

L(M(✓)|S,B) = Y

i2B PB(Ei)

Y

j2S PS(M(✓)|Ej) (9)

PB(E) =

Z 1

0 ⌧ �B(E

0)RB(E 0, E) dE0 (10)

PS(M(✓)|E) =

Z 1

0 �B(E

0)RB(E 0, E) dE0 + ✓0

Z 1

0 �G(E

0|✓)RG(E 0, E) dE0 (11)

LT (M(✓)|ST , BT ) = Y

i

Li(M(✓)|Si, Bi) (12)

� < 2 p Non TEG/Tobs

Ae↵ TEG (13)

�(l, b) = �CR + �e+e� + �EG�� + �Gal��(l, b) + �EG-DM + �Gal-DM('(l, b)) (14)

�(l1, b1)� �(l2, b2) = �Gal��(l1, b1) + �Gal-DM('1)� �Gal��(l2, b2)� �Gal-DM('2) (15)

�Gal-DM('1) = �(l1, b1)� �(l2, b2)� �Gal��(l1, b1) + �Gal��(l2, b2)

1� �Gal-DM('2) �Gal-DM('1)

(16)

f(x) = 1p 2⇡�

e� 1 2(

x�µ � )2 (17)

µ̂(x1, ..., xn) (18)

�̂(x1, ..., xn) (19)

lim n!1

P (|✓̂ � ✓| > ✏) = 0 (20)

2

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Expectation value and bias ★ Estimators, as function of random variables (the measurements), are also

random variables themselves

★ I.e. if the entire experiment was repeated m times we would obtain m values (𝜃1, …,𝜃m), following some pdf g(𝜃;𝜃), which depends in general on the true value 𝜃.

★ The pdf of an estimator (g in this case) is called sampling distribution

★ Expectation value of an estimator 𝜃 with sampling pdf g(𝜃;𝜃):

★ Bias of the estimator 𝜃 is defined as:

★ Note that: ✦ The bias does not depend on the measured values

✦ In general, it depends on the size of the sample, the functional form of the estimator and the true properties of the estimator (the true value 𝜃 among them)

7

^

d�PP

dE =

1

mdm⌧dm

dN�

dE (5)

J(⌦) = 1

4⇡

Z

⌦

Z

los ⇢(l) dl d⌦ (6)

⇢(r) = ⇢s

(r/rs)(1 + r/rs)2 (7)

L(s, b|n,m) = (s+ b)n

n! e�(s+b) (⌧b)

m

m! e�⌧b (8)

L(M(✓)|S,B) = Y

i2B PB(Ei)

Y

j2S PS(M(✓)|Ej) (9)

PB(E) =

Z 1

0 ⌧ �B(E

0)RB(E 0, E) dE0 (10)

PS(M(✓)|E) =

Z 1

0 �B(E

0)RB(E 0, E) dE0 + ✓0

Z 1

0 �G(E

0|✓)RG(E 0, E) dE0 (11)

LT (M(✓)|ST , BT ) = Y

i

Li(M(✓)|Si, Bi) (12)

� < 2 p

Non TEG/Tobs

Ae↵ TEG (13)

�(l, b) = �CR + �e+e� + �EG�� + �Gal��(l, b) + �EG-DM + �Gal-DM('(l, b)) (14)

�(l1, b1)� �(l2, b2) = �Gal��(l1, b1) + �Gal-DM('1)� �Gal��(l2, b2)� �Gal-DM('2) (15)

�Gal-DM('1) = �(l1, b1)� �(l2, b2)� �Gal��(l1, b1) + �Gal��(l2, b2)

1� �Gal-DM('2) �Gal-DM('1)

(16)

f(x) = 1p 2⇡�

e� 1 2(

x�µ � )2 (17)

µ̂(x1, ..., xn) (18)

�̂(x1, ..., xn) (19)

lim n!1

P (|✓̂ � ✓| > ✏) = 0 (20)

E[✓̂(x)] =

Z d✓̂ ✓̂ g(✓̂; ✓) =

Z dx1...

Z dxn ✓̂(x) f(x1; ✓) ... f(xn; ✓) (21)

2

d�PP

dE =

1

mdm⌧dm

dN�

dE (5)

J(⌦) = 1

4⇡

Z

⌦

Z

los ⇢(l) dl d⌦ (6)

⇢(r) = ⇢s

(r/rs)(1 + r/rs)2 (7)

L(s, b|n,m) = (s+ b)n

n! e�(s+b) (⌧b)

m

m! e�⌧b (8)

L(M(✓)|S,B) = Y

i2B PB(Ei)

Y

j2S PS(M(✓)|Ej) (9)

PB(E) =

Z 1

0 ⌧ �B(E

0)RB(E 0, E) dE0 (10)

PS(M(✓)|E) =

Z 1

0 �B(E

0)RB(E 0, E) dE0 + ✓0

Z 1

0 �G(E

0|✓)RG(E 0, E) dE0 (11)

LT (M(✓)|ST , BT ) = Y

i

Li(M(✓)|Si, Bi) (12)

� < 2 p Non TEG/Tobs

Ae↵ TEG (13)

�(l, b) = �CR + �e+e� + �EG�� + �Gal��(l, b) + �EG-DM + �Gal-DM('(l, b)) (14)

�(l1, b1)� �(l2, b2) = �Gal��(l1, b1) + �Gal-DM('1)� �Gal��(l2, b2)� �Gal-DM('2) (15)

�Gal-DM('1) = �(l1, b1)� �(l2, b2)� �Gal��(l1, b1) + �Gal��(l2, b2)

1� �Gal-DM('2) �Gal-DM('1)

(16)

f(x) = 1p 2⇡�

e� 1 2(

x�µ � )2 (17)

µ̂(x1, ..., xn) (18)

�̂(x1, ..., xn) (19)

lim n!1

P (|✓̂ � ✓| > ✏) = 0 (20)

E[✓̂(x)] =

Z d✓̂ ✓̂ g(✓̂; ✓) =

Z dx1...

Z dxn ✓̂(x) f(x1; ✓) ... f(xn; ✓) (21)

b = E[✓̂]� ✓ (22)

2

^ ^

^^

^

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Classification of estimator according to bias

★ If b=0 we call it an unbiased estimator

★ If b=0 only for n→∞ we call the estimator asymptotically unbiased

★ An estimator can be consistent AND biased (e.g as estimator of the mean)

★ An estimator can be unbiased and non-consistent (e.g. x1 for {x1,… xn} as estimator of the mean)

★ Unbiased estimators are particularly useful when combining several measurements

★ Normally a bias smaller than statistical error is accepted

★ Another measure of the quality of an estimator is provided by the mean squared error (MSE) [✪]:

★ Classical statistics does not provide unique method to define an estimator, but once defined we can characterized by variance, bias and/or MSE

8

onN

0 10 20 30 40 50 60 70 80 90 100

310×

=0 )

sy st

τ ∆

Se ns

iti vi

ty /

Se ns

iti vi

ty (

1−10

1

=3.0systτ∆

=2.0systτ∆

=1.0systτ∆

=0.0systτ∆

=1τ

onN

0 10 20 30 40 50 60 70 80 90 100

310×

=0 )

sy st

τ ∆

Se ns

iti vi

ty /

Se ns

iti vi

ty (

1−10

1

=3τ

onN

0 10 20 30 40 50 60 70 80 90 100

310×

=0 )

sy st

τ ∆

Se ns

iti vi

ty /

Se ns

iti vi

ty (

1−10

1

=10τ

τ

0 1 2 3 4 5 6 7 8 9 10

=0 )

sy st

τ ∆

=1 ,

τ Se

ns iti

vi ty

/ Se

ns iti

vi ty

(

1−10

1

10

=10.0systτ∆

=5.0systτ∆

=3.0systτ∆ =1.5systτ∆

=0.0systτ∆

=10000onN

τ

0 1 2 3 4 5 6 7 8 9 10 =0

) sy

st τ

∆ =1

, τ

Se ns

iti vi

ty /

Se ns

iti vi

ty (

1−10

1

10

=50000onN

τ

0 1 2 3 4 5 6 7 8 9 10

=0 )

sy st

τ ∆

=1 ,

τ Se

ns iti

vi ty

/ Se

ns iti

vi ty

(

1−10

1

10

=100000onN

b = E[✓̂]� ✓ (22)

MSE = E[(✓̂ � ✓)2] = E[(✓̂ � E[✓̂])2] + (E[✓̂ � ✓])2 = V [✓̂] + b2 (23)

3

S2 = 1

n

nX

i=1

(xi � µ)2 = x2 � µ2 (29)

1

n

nX

i=1

xi + 1

n (30)

5

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Sample mean as estimator of population mean

★ Suppose a sample of size n of a random variable x: x1,…,xn

★ x is distributed according to an unknown pdf f(x)

★ We want to build an estimator of the expectation value of x (𝜇≡E[x]), also known as the population or distribution mean)

★ One possibility is to use the arithmetic mean of the xi sample (called sample mean), i.e:

★ Weak law of large numbers: if V[x] exists (true if x represents a physical quantity), then x is a consistent estimator of 𝜇 (i.e: n→∞ ⇒ x→𝜇) [✪]

★ Expectation value of x:

i.e., x is an unbiased estimator of for the mean of any pdf ★ Variance of x [✪]: (√V(𝜃) is typically referred to as the error of 𝜃)

9

onN

0 10 20 30 40 50 60 70 80 90 100

310×

=0 )

sy st

τ ∆

Se ns

iti vi

ty /

Se ns

iti vi

ty (

1−10

1

=3.0systτ∆

=2.0systτ∆

=1.0systτ∆

=0.0systτ∆

=1τ

onN

0 10 20 30 40 50 60 70 80 90 100

310×

=0 )

sy st

τ ∆

Se ns

iti vi

ty /

Se ns

iti vi

ty (

1−10

1

=3τ

onN

0 10 20 30 40 50 60 70 80 90 100

310×

=0 )

sy st

τ ∆

Se ns

iti vi

ty /

Se ns

iti vi

ty (

1−10

1

=10τ

τ

0 1 2 3 4 5 6 7 8 9 10

=0 )

sy st

τ ∆

=1 ,

τ Se

ns iti

vi ty

/ Se

ns iti

vi ty

(

1−10

1

10

=10.0systτ∆

=5.0systτ∆

=3.0systτ∆ =1.5systτ∆

=0.0systτ∆

=10000onN

τ

0 1 2 3 4 5 6 7 8 9 10

=0 )

sy st

τ ∆

=1 ,

τ Se

ns iti

vi ty

/ Se

ns iti

vi ty

(

1−10

1

10

=50000onN

τ

0 1 2 3 4 5 6 7 8 9 10

=0 )

sy st

τ ∆

=1 ,

τ Se

ns iti

vi ty

/ Se

ns iti

vi ty

(

1−10

1

10

=100000onN

b = E[✓̂]� ✓ (22)

MSE = E[(✓̂ � ✓)2] = E[(✓̂ � E[✓̂])2] + (E[✓̂ � ✓])2 = V [✓̂] + b2 (23)

x̄ = 1

n

nX

i=1

xi (24)

3

onN

0 10 20 30 40 50 60 70 80 90 100

310×

=0 )

sy st

τ ∆

Se ns

iti vi

ty /

Se ns

iti vi

ty (

1−10

1

=3.0systτ∆

=2.0systτ∆

=1.0systτ∆

=0.0systτ∆

=1τ

onN

0 10 20 30 40 50 60 70 80 90 100

310×

=0 )

sy st

τ ∆

Se ns

iti vi

ty /

Se ns

iti vi

ty (

1−10

1

=3τ

onN

0 10 20 30 40 50 60 70 80 90 100

310×

=0 )

sy st

τ ∆

Se ns

iti vi

ty /

Se ns

iti vi

ty (

1−10

1

=10τ

τ

0 1 2 3 4 5 6 7 8 9 10

=0 )

sy st

τ ∆

=1 ,

τ Se

ns iti

vi ty

/ Se

ns iti

vi ty

(

1−10

1

10

=10.0systτ∆

=5.0systτ∆

=3.0systτ∆ =1.5systτ∆

=0.0systτ∆

=10000onN

τ

0 1 2 3 4 5 6 7 8 9 10

=0 )

sy st

τ ∆

=1 ,

τ Se

ns iti

vi ty

/ Se

ns iti

vi ty

(

1−10

1

10

=50000onN

τ

0 1 2 3 4 5 6 7 8 9 10

=0 )

sy st

τ ∆

=1 ,

τ Se

ns iti

vi ty

/ Se

ns iti

vi ty

(

1−10

1

10

=100000onN

b = E[✓̂]� ✓ (22)

MSE = E[(✓̂ � ✓)2] = E[(✓̂ � E[✓̂])2] + (E[✓̂ � ✓])2 = V [✓̂] + b2 (23)

x̄ = 1

n

nX

i=1

xi (24)

E[x̄] = E

" 1

n

nX

i=1

xi

# =

1

n

nX

i=1

E[xi] = 1

n

nX

i=n

µ = µ (25)

3

_

_

_

(%)systτ ∆ 0 0.5 1 1.5 2 2.5 3

=0 )

sy st

τ ∆

Se ns

iti vi

ty /

Se ns

iti vi

ty (

1

1.2

1.4

1.6

1.8

2

2.2

2.4

2.6

2.8

3

3.2

3.4

3.6 = 100 GeVDMm

=5165)on=25h (NobsT

=1 (theoretical)binsN

=1binsN

=2binsN

=10binsN

=0binsN

= 100 GeVDMm

(%)systτ ∆ 0 0.5 1 1.5 2 2.5 3

=0 )

sy st

τ ∆

Se ns

iti vi

ty /

Se ns

iti vi

ty (

1

1.2

1.4

1.6

1.8

2

2.2

2.4

2.6

2.8

3

3.2

3.4

3.6 = 1000 GeVDMm

=5165)on=25h (NobsT

=1 (theoretical)binsN

=1binsN

=2binsN

=10binsN

=0binsN

= 1000 GeVDMm

(%)systτ ∆ 0 0.5 1 1.5 2 2.5 3

=0 )

sy st

τ ∆

Se ns

iti vi

ty /

Se ns

iti vi

ty (

1

1.2

1.4

1.6

1.8

2

2.2

2.4

2.6

2.8

3

3.2

3.4

3.6 = 10000 GeVDMm

=5165)on=25h (NobsT

=1 (theoretical)binsN

=1binsN

=2binsN

=10binsN

=0binsN

= 10000 GeVDMm

(%)systτ ∆ 0 0.5 1 1.5 2 2.5 3

=0 )

sy st

τ ∆

Se ns

iti vi

ty /

Se ns

iti vi

ty (

1

1.2

1.4

1.6

1.8

2

2.2

2.4

2.6

2.8

3

3.2

3.4

3.6 = 100 GeVDMm

=16460)on=80h (NobsT

=1 (theoretical)binsN

=1binsN

=2binsN

=10binsN

=0binsN

= 100 GeVDMm

(%)systτ ∆ 0 0.5 1 1.5 2 2.5 3

=0 )

sy st

τ ∆

Se ns

iti vi

ty /

Se ns

iti vi

ty (

1

1.2

1.4

1.6

1.8

2

2.2

2.4

2.6

2.8

3

3.2

3.4

3.6 = 1000 GeVDMm

=16460)on=80h (NobsT

=1 (theoretical)binsN

=1binsN

=2binsN

=10binsN

=0binsN

= 1000 GeVDMm

(%)systτ ∆ 0 0.5 1 1.5 2 2.5 3

=0 )

sy st

τ ∆

Se ns

iti vi

ty /

Se ns

iti vi

ty (

1

1.2

1.4

1.6

1.8

2

2.2

2.4

2.6

2.8

3

3.2

3.4

3.6 = 10000 GeVDMm

=16460)on=80h (NobsT

=1 (theoretical)binsN

=1binsN

=2binsN

=10binsN

=0binsN

= 10000 GeVDMm

(%)systτ ∆ 0 0.5 1 1.5 2 2.5 3

=0 )

sy st

τ ∆

Se ns

iti vi

ty /

Se ns

iti vi

ty (

1

1.2

1.4

1.6

1.8

2

2.2

2.4

2.6

2.8

3

3.2

3.4

3.6 = 100 GeVDMm

=41067)on=200h (NobsT

=1 (theoretical)binsN

=1binsN

=2binsN

=10binsN

=0binsN

= 100 GeVDMm

(%)systτ ∆ 0 0.5 1 1.5 2 2.5 3

=0 )

sy st

τ ∆

Se ns

iti vi

ty /

Se ns

iti vi

ty (

1

1.2

1.4

1.6

1.8

2

2.2

2.4

2.6

2.8

3

3.2

3.4

3.6 = 1000 GeVDMm

=41067)on=200h (NobsT

=1 (theoretical)binsN

=1binsN

=2binsN

=10binsN

=0binsN

= 1000 GeVDMm

(%)systτ ∆ 0 0.5 1 1.5 2 2.5 3

=0 )

sy st

τ ∆

Se ns

iti vi

ty /

Se ns

iti vi

ty (

1

1.2

1.4

1.6

1.8

2

2.2

2.4

2.6

2.8

3

3.2

3.4

3.6 = 10000 GeVDMm

=41067)on=200h (NobsT

=1 (theoretical)binsN

=1binsN

=2binsN

=10binsN

=0binsN

= 10000 GeVDMm

V [x̄] = �2

n (28)

4

^ ^

_

_

x

f(x)

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Estimators of population variance ★ If the population mean 𝜇 is known, then the statistic S2, defined by

is an unbiased estimator of 𝜎2 [✪] (variance of f(x), also known as population variance 𝜎2 = E[(x-𝜇)2]), i.e. E[s2] = 𝜎2

★ Normally we do not know 𝜇, though. The sample variance s2, defined as:

is an unbiased estimator of 𝜎2 [✪]

★ The variance of the estimator s2 can be computed to be

where 𝜇k is the k-th central moment (e.g. 𝜇2 = 𝜎2), that can be estimated by

10

onN

0 10 20 30 40 50 60 70 80 90 100

310×

=0 )

sy st

τ ∆

Se ns

iti vi

ty /

Se ns

iti vi

ty (

1−10

1

=3.0systτ∆

=2.0systτ∆

=1.0systτ∆

=0.0systτ∆

=1τ

onN

0 10 20 30 40 50 60 70 80 90 100

310×

=0 )

sy st

τ ∆

Se ns

iti vi

ty /

Se ns

iti vi

ty (

1−10

1

=3τ

onN

0 10 20 30 40 50 60 70 80 90 100

310×

=0 )

sy st

τ ∆

Se ns

iti vi

ty /

Se ns

iti vi

ty (

1−10

1

=10τ

τ

0 1 2 3 4 5 6 7 8 9 10

=0 )

sy st

τ ∆

=1 ,

τ Se

ns iti

vi ty

/ Se

ns iti

vi ty

(

1−10

1

10

=10.0systτ∆

=5.0systτ∆

=3.0systτ∆ =1.5systτ∆

=0.0systτ∆

=10000onN

τ

0 1 2 3 4 5 6 7 8 9 10

=0 )

sy st

τ ∆

=1 ,

τ Se

ns iti

vi ty

/ Se

ns iti

vi ty

(

1−10

1

10

=50000onN

τ

0 1 2 3 4 5 6 7 8 9 10

=0 )

sy st

τ ∆

=1 ,

τ Se

ns iti

vi ty

/ Se

ns iti

vi ty

(

1−10

1

10

=100000onN

b = E[✓̂]� ✓ (22)

MSE = E[(✓̂ � ✓)2] = E[(✓̂ � E[✓̂])2] + (E[✓̂ � ✓])2 = V [✓̂] + b2 (23)

x̄ = 1

n

nX

i=1

xi (24)

E[x̄] = E

" 1

n

nX

i=1

xi

# =

1

n

nX

i=1

E[xi] = 1

n

nX

i=n

µ = µ (25)

s2 = 1

n� 1

nX

i=1

(xi � x̄)2 = n

n� 1 (x2 � x̄2) (26)

3

S2 = 1

n

nX

i=1

(xi � µ)2 = x2 � µ2 (29)

1

n

nX

i=1

xi + 1

n (30)

V̂xy = 1

n� 1

nX

i=1

(xi � x̄)(yi � ȳ) = n

n� 1 (xy � x̄ȳ) (31)

rxy = V̂xy

sxsy =

P n

i=1(xi � x̄)(yi � ȳ)qP n

j=1(xj � x̄)2 P

n

k=1(yk � ȳ)2 =

xy � x̄ȳq (x2 � x̄2)(y2 � ȳ2)

(32)

x̄ = nX

i=1

wixi (33)

E[x̄] = µ ) nX

i=1

wiE[xi] = µ ) nX

i=1

wiµ = µ ) µ nX

i=1

wi = µ ) nX

i=1

wi = 1 (34)

wi = 1/�2

iP n

j=1 1/� 2 j

(35)

V [s2] = 1

n

✓ µ4 �

n� 3

n� 1 µ2 2

◆ (36)

5

S2 = 1

n

nX

i=1

(xi � µ)2 = x2 � µ2 (29)

1

n

nX

i=1

xi + 1

n (30)

V̂xy = 1

n� 1

nX

i=1

(xi � x̄)(yi � ȳ) = n

n� 1 (xy � x̄ȳ) (31)

rxy = V̂xy

sxsy =

P n

i=1(xi � x̄)(yi � ȳ)qP n

j=1(xj � x̄)2 P

n

k=1(yk � ȳ)2 =

xy � x̄ȳq (x2 � x̄2)(y2 � ȳ2)

(32)

x̄ = nX

i=1

wixi (33)

E[x̄] = µ ) nX

i=1

wiE[xi] = µ ) nX

i=1

wiµ = µ ) µ nX

i=1

wi = µ ) nX

i=1

wi = 1 (34)

wi = 1/�2

iP n

j=1 1/� 2 j

(35)

V [s2] = 1

n

✓ µ4 �

n� 3

n� 1 µ2 2

◆ (36)

mk = 1

n� 1

nX

i=1

(xi � x̄)k (37)

5

E[x̄] = E

" 1

n

nX

i=1

xi

# =

1

n

nX

i=1

E[xi] = 1

n

nX

i=n

µ = µ (60)

s2 = 1

n� 1

nX

i=1

(xi � x̄)2 = n

n� 1 (x2 � x̄2) (61)

�2 = E[(x� µ)2] (62)

V [x̄] = �2

n (63)

S2 = 1

n

nX

i=1

(xi � µ)2 (64)

1

n

nX

i=1

xi + 1

n (65)

V̂xy = 1

n� 1

nX

i=1

(xi � x̄)(yi � ȳ) = n

n� 1 (xy � x̄ȳ) (66)

rxy = V̂xy

sxsy =

P n

i=1(xi � x̄)(yi � ȳ)qP n

j=1(xj � x̄)2 P

n

k=1(yk � ȳ)2 =

xy � x̄ȳq (x2 � x̄2)(y2 � ȳ2)

(67)

x̄ = nX

i=1

wixi (68)

E[x̄] = µ ) nX

i=1

wiE[xi] = µ ) nX

i=1

wiµ = µ ) µ nX

i=1

wi = µ ) nX

i=1

wi = 1 (69)

wi = 1/�2

iP n

j=1 1/� 2 j

(70)

V [s2] = 1

n

✓ µ4 �

n� 3

n� 1 µ2 2

◆ (71)

mk = 1

n� 1

nX

i=1

(xi � x̄)k (72)

E[⌘i(✓)] ⇡ ⌘(µ) (73)

Uij ⇡ nX

k,l=1

@⌘i @✓k

@⌘j @✓l

���� ✓=µ

Vkl (74)

8

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Estimator of the correlation coefficient

★ In a similar way, one can show that the covariance between two measured sets of variables (x1,…xn) and (y1,…,yn) (with e.g. x=height, y=water consumption)

is an unbiased estimator of the covariance Vxy between two random variables x and y of unknown mean

And then, an estimator r of the correlation coefficient 𝜌 can be written as:

11

S2 = 1

n

nX

i=1

(xi � µ)2 = x2 � µ2 (29)

1

n

nX

i=1

xi + 1

n (30)

V̂xy = 1

n� 1

nX

i=1

(xi � x̄)(yi � ȳ) = n

n� 1 (xy � x̄ȳ) (31)

rxy = V̂xy

sxsy =

P n

i=1(xi � x̄)(yi � ȳ)qP n

j=1(xj � x̄)2 P

n

k=1(yk � ȳ)2 =

xy � x̄ȳq (x2 � x̄2)(y2 � ȳ2)

(32)

5

S2 = 1

n

nX

i=1

(xi � µ)2 = x2 � µ2 (29)

1

n

nX

i=1

xi + 1

n (30)

V̂xy = 1

n� 1

nX

i=1

(xi � x̄)(yi � ȳ) = n

n� 1 (xy � x̄ȳ) (31)

5

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Weighted mean ★ Suppose that we have a set of observations (x1,…,xn) where each xi is distributed according to

a different fi(xi), all with a common mean 𝜇 but different variance 𝜎i2

★ How do we optimally combine the xi into an estimator of 𝜇? ✦ typical problem when we have several experiments measuring the same parameter and want to get the

global result

✦ E.g. the height of a person measured with different instruments

★ We generalize the arithmetic mean to the weighted mean:

with wi called weights

★ For x to be an unbiased estimator, then:

it suffices that the xi estimators are unbiased and that the sum of weights is 1

★ Weights are chosen to minimize the variance of the estimator x, i.e.[✪]:

12

S2 = 1

n

nX

i=1

(xi � µ)2 = x2 � µ2 (29)

1

n

nX

i=1

xi + 1

n (30)

V̂xy = 1

n� 1

nX

i=1

(xi � x̄)(yi � ȳ) = n

n� 1 (xy � x̄ȳ) (31)

rxy = V̂xy

sxsy =

P n

i=1(xi � x̄)(yi � ȳ)qP n

j=1(xj � x̄)2 P

n

k=1(yk � ȳ)2 =

xy � x̄ȳq (x2 � x̄2)(y2 � ȳ2)

(32)

x̄ = NX

i=1

wixi (33)

E[x̄] = µ ) NX

i=1

wiE[xi] = µ ) NX

i=1

wiµ = µ ) µ NX

i=1

wi = µ ) NX

i=1

wi = 1 (34)

wi = 1/�2

iP n

j=1 1/� 2 j

(35)

5

S2 = 1

n

nX

i=1

(xi � µ)2 = x2 � µ2 (29)

1

n

nX

i=1

xi + 1

n (30)

V̂xy = 1

n� 1

nX

i=1

(xi � x̄)(yi � ȳ) = n

n� 1 (xy � x̄ȳ) (31)

rxy = V̂xy

sxsy =

P n

i=1(xi � x̄)(yi � ȳ)qP n

j=1(xj � x̄)2 P

n

k=1(yk � ȳ)2 =

xy � x̄ȳq (x2 � x̄2)(y2 � ȳ2)

(32)

x̄ = nX

i=1

wixi (33)

E[x̄] = µ ) nX

i=1

wiE[xi] = µ ) nX

i=1

wiµ = µ ) µ nX

i=1

wi = µ ) nX

i=1

wi = 1 (34)

wi = 1/�2

iP n

j=1 1/� 2 j

(35)

5

S2 = 1

n

nX

i=1

(xi � µ)2 = x2 � µ2 (29)

1

n

nX

i=1

xi + 1

n (30)

V̂xy = 1

n� 1

nX

i=1

(xi � x̄)(yi � ȳ) = n

n� 1 (xy � x̄ȳ) (31)

rxy = V̂xy

sxsy =

P n

i=1(xi � x̄)(yi � ȳ)qP n

j=1(xj � x̄)2 P

n

k=1(yk � ȳ)2 =

xy � x̄ȳq (x2 � x̄2)(y2 � ȳ2)

(32)

x̄ = nX

i=1

wixi (33)

E[x̄] = µ ) nX

i=1

wiE[xi] = µ ) nX

i=1

wiµ = µ ) µ nX

i=1

wi = µ ) nX

i=1

wi = 1 (34)

wi = 1/�2

iP n

j=1 1/� 2 j

(35)

5

_

_

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Exercises 1. Read x, y values (e.g. x=height, y=water consumption) from input file

ParEst_input_1.txt and estimate:

a. the population mean and variance of the x and y values (considered separately), with errors;

b. the correlation coefficient between x and y

2. Consider a Gaussian pdf with 𝜇=0 and 𝜎=1. Show numerically that: a. the arithmetic mean is a consistent and unbiased (e.g. for n=10) estimator

of 𝜇;

b. s2 (see the slides) is a consistent and unbiased (e.g. for n=10) estimator of 𝜎2, but that it is biased if we change the normalization factor 1/(n-1) by 1/n

3. Read x, 𝜟x values from input file ParEst_input_3.txt, and compute:

a. The weighted mean

b. Its variance

13

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Least-squares

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

The least-squares (LS) method ★ Consider a set of measured values yi, each distributed with Gaussian pdf

centred around 𝜆i (this is true in many practical situations as expected from the Central Limit theorem)

★ Assume that we have N independent yi measurements, which depend on xi (known without error). ✦ E.g. xi is the age in years and yi = 𝜆(xi) the average height of the population with that

age

★ Each yi has different unknown mean 𝜆i and different known variance 𝜎i2

★ Suppose the true value is given by the function of x 𝜆(x;𝜽), with unknown parameters 𝜽=(𝜃1,…,𝜃m), to be estimated

★ This can be done by finding the values of 𝜽 that minimize the quantity

which is called the LS method

15

S2 = 1

n

nX

i=1

(xi � µ)2 = x2 � µ2 (29)

1

n

nX

i=1

xi + 1

n (30)

V̂xy = 1

n� 1

nX

i=1

(xi � x̄)(yi � ȳ) = n

n� 1 (xy � x̄ȳ) (31)

rxy = V̂xy

sxsy =

P n

i=1(xi � x̄)(yi � ȳ)qP n

j=1(xj � x̄)2 P

n

k=1(yk � ȳ)2 =

xy � x̄ȳq (x2 � x̄2)(y2 � ȳ2)

(32)

x̄ = nX

i=1

wixi (33)

E[x̄] = µ ) nX

i=1

wiE[xi] = µ ) nX

i=1

wiµ = µ ) µ nX

i=1

wi = µ ) nX

i=1

wi = 1 (34)

wi = 1/�2

iP n

j=1 1/� 2 j

(35)

�2(✓) = NX

i=1

(yi � �(xi;✓))2

�2 i

(36)

5

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Elements of the LS method ★ Example of LS method

elements: ✦ N=9 independent

measurements y1,…yN

✦ Each measurement with known error 𝜎1,…,𝜎N

✦ Measurements corresponding to the values x1,…, xN, known without errors

✦ True value 𝜆i of yi assumed to be given by function 𝜆i = 𝜆(xi;𝜽)

✦ The LS estimators of 𝜽 are found by minimizing 𝜒2

16

S2 = 1

n

nX

i=1

(xi � µ)2 = x2 � µ2 (29)

1

n

nX

i=1

xi + 1

n (30)

V̂xy = 1

n� 1

nX

i=1

(xi � x̄)(yi � ȳ) = n

n� 1 (xy � x̄ȳ) (31)

rxy = V̂xy

sxsy =

P n

i=1(xi � x̄)(yi � ȳ)qP n

j=1(xj � x̄)2 P

n

k=1(yk � ȳ)2 =

xy � x̄ȳq (x2 � x̄2)(y2 � ȳ2)

(32)

x̄ = nX

i=1

wixi (33)

E[x̄] = µ ) nX

i=1

wiE[xi] = µ ) nX

i=1

wiµ = µ ) µ nX

i=1

wi = µ ) nX

i=1

wi = 1 (34)

wi = 1/�2

iP n

j=1 1/� 2 j

(35)

�2(✓) = NX

i=1

(yi � �(xi;✓))2

�2 i

(36)

�(x;✓) = ✓0 ⇣ e�✓1x cos(2⇡✓2x)

⌘ (37)

yi ± �i (38)

5

L(✓) = nY

i=1

f(xi;✓) (44)

�(x;✓) = ✓0 + ✓1x

(1 + ✓2x+ ✓3x2) (45)

L(µ;n) = NY

i=1

µni

ni! e�µ (46)

6

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

LS for correlated measurements ★ If the measurements yi are not independent from each other, but described

by an N-dimensional Gaussian with covariance matrix V and unknown mean values, the parameters 𝜽 are found by minimizing:

✦ Both expressions used in practice even if the pdf are not Gaussian ✦ Note that we do not need to know the pdf, just the covariance matrix! ✦ Eg. of correlated measurements: x = age; y = average height of the

population over that age

★ The parameters that minimize 𝜒2 are called least-squares estimators, 𝜃1, …,𝜃m

★ The resulting minimum value 𝜒2min follows under certain circumstances the 𝜒2 distribution, justifying its name. In general though, it does not follow the 𝜒2 distribution (more later)

17

S2 = 1

n

nX

i=1

(xi � µ)2 = x2 � µ2 (29)

1

n

nX

i=1

xi + 1

n (30)

V̂xy = 1

n� 1

nX

i=1

(xi � x̄)(yi � ȳ) = n

n� 1 (xy � x̄ȳ) (31)

rxy = V̂xy

sxsy =

P n

i=1(xi � x̄)(yi � ȳ)qP n

j=1(xj � x̄)2 P

n

k=1(yk � ȳ)2 =

xy � x̄ȳq (x2 � x̄2)(y2 � ȳ2)

(32)

x̄ = nX

i=1

wixi (33)

E[x̄] = µ ) nX

i=1

wiE[xi] = µ ) nX

i=1

wiµ = µ ) µ nX

i=1

wi = µ ) nX

i=1

wi = 1 (34)

wi = 1/�2

iP n

j=1 1/� 2 j

(35)

�2(✓) = NX

i=1

(yi � �(xi;✓))2

�2 i

(36)

�(x;✓) = ✓0 ⇣ e�✓1x cos(2⇡✓2x)

⌘ (37)

yi ± �i (38)

�2(✓) = NX

i,j=1

(yi � �(xi;✓))(V �1)ij(yj � �(xj ;✓)) (39)

5

^ ^

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Linear LS fit ★ Assume the function 𝜆(x;𝜽) is linear in the 𝜽 parameters, i.e:

with aj(x) linearly independent functions of x

★ The value of 𝜆 evaluated at xi can be written as

★ And the expression of 𝜒2 in matrix notation is then:

★ To find the estimators 𝜽, we minimize 𝜒2 wrt 𝜽:

which results in:

(i.e. linear functions of the measurements y)

18

S2 = 1

n

nX

i=1

(xi � µ)2 = x2 � µ2 (29)

1

n

nX

i=1

xi + 1

n (30)

V̂xy = 1

n� 1

nX

i=1

(xi � x̄)(yi � ȳ) = n

n� 1 (xy � x̄ȳ) (31)

rxy = V̂xy

sxsy =

P n

i=1(xi � x̄)(yi � ȳ)qP n

j=1(xj � x̄)2 P

n

k=1(yk � ȳ)2 =

xy � x̄ȳq (x2 � x̄2)(y2 � ȳ2)

(32)

x̄ = nX

i=1

wixi (33)

E[x̄] = µ ) nX

i=1

wiE[xi] = µ ) nX

i=1

wiµ = µ ) µ nX

i=1

wi = µ ) nX

i=1

wi = 1 (34)

wi = 1/�2

iP n

j=1 1/� 2 j

(35)

�2(✓) = NX

i=1

(yi � �(xi;✓))2

�2 i

(36)

�(x;✓) = ✓0 ⇣ e�✓1x cos(2⇡✓2x)

⌘ (37)

yi ± �i (38)

�2(✓) = NX

i,j=1

(yi � �(xi;✓))(V �1)ij(yj � �(xj ;✓)) (39)

�(x;✓) = mX

j=1

aj(x)✓j (40)

5

S2 = 1

n

nX

i=1

(xi � µ)2 = x2 � µ2 (29)

1

n

nX

i=1

xi + 1

n (30)

V̂xy = 1

n� 1

nX

i=1

(xi � x̄)(yi � ȳ) = n

n� 1 (xy � x̄ȳ) (31)

rxy = V̂xy

sxsy =

P n

i=1(xi � x̄)(yi � ȳ)qP n

j=1(xj � x̄)2 P

n

k=1(yk � ȳ)2 =

xy � x̄ȳq (x2 � x̄2)(y2 � ȳ2)

(32)

x̄ = nX

i=1

wixi (33)

E[x̄] = µ ) nX

i=1

wiE[xi] = µ ) nX

i=1

wiµ = µ ) µ nX

i=1

wi = µ ) nX

i=1

wi = 1 (34)

wi = 1/�2

iP n

j=1 1/� 2 j

(35)

�2(✓) = NX

i=1

(yi � �(xi;✓))2

�2 i

(36)

�(x;✓) = ✓0 ⇣ e�✓1x cos(2⇡✓2x)

⌘ (37)

yi ± �i (38)

�2(✓) = NX

i,j=1

(yi � �(xi;✓))(V �1)ij(yj � �(xj ;✓)) (39)

�(x;✓) = mX

j=1

aj(x)✓j (40)

�(xi;✓) = mX

j=1

aj(xi)✓j ⌘ mX

j=1

Aij✓j (41)

�2 = (y � �)TV �1(y � �) = (y �A✓)TV �1(y �A✓) (42)

5

S2 = 1

n

nX

i=1

(xi � µ)2 = x2 � µ2 (29)

1

n

nX

i=1

xi + 1

n (30)

V̂xy = 1

n� 1

nX

i=1

(xi � x̄)(yi � ȳ) = n

n� 1 (xy � x̄ȳ) (31)

rxy = V̂xy

sxsy =

P n

i=1(xi � x̄)(yi � ȳ)qP n

j=1(xj � x̄)2 P

n

k=1(yk � ȳ)2 =

xy � x̄ȳq (x2 � x̄2)(y2 � ȳ2)

(32)

x̄ = nX

i=1

wixi (33)

E[x̄] = µ ) nX

i=1

wiE[xi] = µ ) nX

i=1

wiµ = µ ) µ nX

i=1

wi = µ ) nX

i=1

wi = 1 (34)

wi = 1/�2

iP n

j=1 1/� 2 j

(35)

�2(✓) = NX

i=1

(yi � �(xi;✓))2

�2 i

(36)

�(x;✓) = ✓0 ⇣ e�✓1x cos(2⇡✓2x)

⌘ (37)

yi ± �i (38)

�2(✓) = NX

i,j=1

(yi � �(xi;✓))(V �1)ij(yj � �(xj ;✓)) (39)

�(x;✓) = mX

j=1

aj(x)✓j (40)

�(xi;✓) = mX

j=1

aj(xi)✓j ⌘ mX

j=1

Aij✓j (41)

�2 = (y � �)TV �1(y � �) = (y �A✓)TV �1(y �A✓) (42)

r�2 = �2(ATV �1y �ATV �1A✓) = 0 (43)

5

S2 = 1

n

nX

i=1

(xi � µ)2 = x2 � µ2 (29)

1

n

nX

i=1

xi + 1

n (30)

V̂xy = 1

n� 1

nX

i=1

(xi � x̄)(yi � ȳ) = n

n� 1 (xy � x̄ȳ) (31)

rxy = V̂xy

sxsy =

P n

i=1(xi � x̄)(yi � ȳ)qP n

j=1(xj � x̄)2 P

n

k=1(yk � ȳ)2 =

xy � x̄ȳq (x2 � x̄2)(y2 � ȳ2)

(32)

x̄ = nX

i=1

wixi (33)

E[x̄] = µ ) nX

i=1

wiE[xi] = µ ) nX

i=1

wiµ = µ ) µ nX

i=1

wi = µ ) nX

i=1

wi = 1 (34)

wi = 1/�2

iP n

j=1 1/� 2 j

(35)

�2(✓) = NX

i=1

(yi � �(xi;✓))2

�2 i

(36)

�(x;✓) = ✓0 ⇣ e�✓1x cos(2⇡✓2x)

⌘ (37)

yi ± �i (38)

�2(✓) = NX

i,j=1

(yi � �(xi;✓))(V �1)ij(yj � �(xj ;✓)) (39)

�(x;✓) = mX

j=1

aj(x)✓j (40)

�(xi;✓) = mX

j=1

aj(xi)✓j ⌘ mX

j=1

Aij✓j (41)

�2 = (y � �)TV �1(y � �) = (y �A✓)TV �1(y �A✓) (42)

r�2 = �2(ATV �1y �ATV �1A✓) = 0 (43)

✓̂ = (ATV �1A)�1ATV �1y ⌘ By (44)

5

^

E[⌘i(✓)] ⇡ ⌘(µ) (75)

Uij ⇡ nX

k,l=1

@⌘i @✓k

@⌘j @✓l

���� ✓=µ

Vkl (76)

U ⇡ AV AT (77)

Aij = @⌘i @✓j

���� ✓=µ

(78)

p =

Z 1

�2 f(z;nd) dz (79)

V =

✓ �2 a cov[a, b]

cov[a, b] �2 b

◆ (80)

probability that xi in [xi, xi + dxi] for all i = nY

i=1

f(xi;✓)dxi (81)

L(✓) = nY

i=1

f(xi;✓) (82)

�(x;✓) = ✓0 + ✓1x

(1 + ✓2x+ ✓3x2) (83)

L(µ;n) = NY

i=1

µni

ni! e�µ ) logL = �Nµ+

NX

i=1

(ni log µ� log ni!) (84)

@ logL

@µ

���� µ=µ̂ML

= 0 ) N � 1

µ̂ML

NX

i=1

ni = 0 ) µ̂ML = 1

N

NX

i=1

ni = n̄ (85)

f(t; ⌧) = 1

⌧ e�

t ⌧ (86)

logL(⌧) = nX

i=1

✓ log

1

⌧ � ti

⌧

◆ (87)

f(x;µ,�) = 1p 2⇡�

e �(x�µ)2

2�2 (88)

�(xi;✓) = mX

j=1

aj(xi)✓j ⌘ mX

j=1

Aij✓j ) � = A✓ (89)

9

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Bias and variance of linear LS estimators

★ Estimators 𝜽 from linear LS fit are unbiased:

★ The covariance matrix for the estimators 𝜽(U) can be computed with error propagation (see later):

★ The inverse of the covariance matrix is given by [✪]:

which coincides with the Rao-Cramér-Frechet limit if measurements y are Gaussian distributed, in which case 𝜒2 = -2logL (see later). Then the 𝜽 are the unbiased estimators with the lowest possible variance

19

^

✓̂ = (ATV �1A)�1ATV �1y ⌘ By (44)

E[✓̂] = E[(ATV �1A)�1ATV �1y] = (ATV �1A)�1ATV �1E[y]

= (ATV �1A)�1ATV �1A✓ = ✓

6

✓̂ = (ATV �1A)�1ATV �1y ⌘ By (44)

E[✓̂] = E[(ATV �1A)�1ATV �1y] = (ATV �1A)�1ATV �1E[y]

= (ATV �1A)�1ATV �1A✓ = ✓

U = B V BT = (ATV �1A)�1 (45)

6

^

✓̂ = (ATV �1A)�1ATV �1y ⌘ By (44)

E[✓̂] = E[(ATV �1A)�1ATV �1y] = (ATV �1A)�1ATV �1E[y]

= (ATV �1A)�1ATV �1A✓ = ✓

U = B V BT = (ATV �1A)�1 (45)

(U�1)ij = 1

2

 @2�2

@✓i@✓j

�

✓=✓̂

(46)

6

^

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Interlude: error propagation ★ Consider a set of n random quantities 𝜽=(𝜃1,…,𝜃n) distributed according to an

unknown pdf f(𝜽), for which we know the expectation value (mean), noted by µ=(µ1, …,µn) (µ=E[𝜽]), and the covariance matrix Vij=cov[𝜃i,𝜃j].

★ Now consider m functions of 𝜽: 𝜼(𝜽)=(𝜂1(𝜽),…,𝜂m(𝜽)).

★ With error propagation we determine the expectation value of 𝜼 and covariance matrix Uij = cov[𝜂i,𝜂j]. For that we Taylor-expand the functions 𝜼(𝜽) about µ, getting [✪]:

The latter can be written in matrix notation as U ≈ A V AT, with

★ The approximation is exact if 𝜼(𝜽) are linear on 𝜽. Otherwise it will break down if 𝜼(𝜽) is significantly non-linear close to µ in a region of size comparable to the standard deviations of 𝜽.

20

E[x̄] = E

" 1

n

nX

i=1

xi

# =

1

n

nX

i=1

E[xi] = 1

n

nX

i=n

µ = µ (60)

s2 = 1

n� 1

nX

i=1

(xi � x̄)2 = n

n� 1 (x2 � x̄2) (61)

�2 = E[(x� µ)2] (62)

V [x̄] = �2

n (63)

S2 = 1

n

nX

i=1

(xi � µ)2 = x2 � µ2 (64)

1

n

nX

i=1

xi + 1

n (65)

V̂xy = 1

n� 1

nX

i=1

(xi � x̄)(yi � ȳ) = n

n� 1 (xy � x̄ȳ) (66)

rxy = V̂xy

sxsy =

P n

i=1(xi � x̄)(yi � ȳ)qP n

j=1(xj � x̄)2 P

n

k=1(yk � ȳ)2 =

xy � x̄ȳq (x2 � x̄2)(y2 � ȳ2)

(67)

x̄ = nX

i=1

wixi (68)

E[x̄] = µ ) nX

i=1

wiE[xi] = µ ) nX

i=1

wiµ = µ ) µ nX

i=1

wi = µ ) nX

i=1

wi = 1 (69)

wi = 1/�2

iP n

j=1 1/� 2 j

(70)

V [s2] = 1

n

✓ µ4 �

n� 3

n� 1 µ2 2

◆ (71)

mk = 1

n� 1

nX

i=1

(xi � x̄)k (72)

E[⌘i(✓)] ⇡ ⌘(µ) (73)

Uij ⇡ nX

k,l=1

@⌘i @✓k

@⌘j @✓l

���� ✓=µ

Vkl (74)

8

U ⇡ AV AT (75)

Aij = @⌘i @✓j

���� ✓=µ

(76)

p =

Z 1

�2 f(z;nd) dz (77)

V =

✓ �2 a cov[a, b]

cov[a, b] �2 b

◆ (78)

probability that xi in [xi, xi + dxi] for all i = nY

i=1

f(xi;✓)dxi (79)

L(✓) = nY

i=1

f(xi;✓) (80)

�(x;✓) = ✓0 + ✓1x

(1 + ✓2x+ ✓3x2) (81)

L(µ;n) = NY

i=1

µni

ni! e�µ ) logL = �Nµ+

NX

i=1

(ni logµ� log ni!) (82)

@ logL

@µ

���� µ=µ̂ML

= 0 ) N � 1

µ̂ML

NX

i=1

ni = 0 ) µ̂ML = 1

N

NX

i=1

ni = n̄ (83)

f(t; ⌧) = 1

⌧ e�

t ⌧ (84)

logL(⌧) = NX

i=1

✓ log

1

⌧ � ti

⌧

◆ (85)

9

E[⌘i(✓)] ⇡ ⌘(E[✓]) = ⌘(µ) (75)

Uij ⇡ nX

k,l=1

@⌘i @✓k

@⌘j @✓l

���� ✓=µ

Vkl (76)

U ⇡ AV AT (77)

Aij = @⌘i @✓j

���� ✓=µ

(78)

p =

Z 1

�2 f(z;nd) dz (79)

V =

✓ �2 a cov[a, b]

cov[a, b] �2 b

◆ (80)

probability that xi in [xi, xi + dxi] for all i = nY

i=1

f(xi;✓)dxi (81)

L(✓) = nY

i=1

f(xi;✓) (82)

�(x;✓) = ✓0 + ✓1x

(1 + ✓2x+ ✓3x2) (83)

L(µ;n) = NY

i=1

µni

ni! e�µ ) logL = �Nµ+

NX

i=1

(ni log µ� log ni!) (84)

@ logL

@µ

���� µ=µ̂ML

= 0 ) N � 1

µ̂ML

NX

i=1

ni = 0 ) µ̂ML = 1

N

NX

i=1

ni = n̄ (85)

f(t; ⌧) = 1

⌧ e�

t ⌧ (86)

logL(⌧) = nX

i=1

✓ log

1

⌧ � ti

⌧

◆ (87)

f(x;µ,�) = 1p 2⇡�

e �(x�µ)2

2�2 (88)

�(xi;✓) = mX

j=1

aj(xi)✓j ⌘ mX

j=1

Aij✓j ) � = A✓ (89)

9

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

𝜒2+1 contours in linear LS fit ★ For the case of 𝜆(x;𝜽) linear in the parameters 𝜽, the 𝜒2 is quadratic in 𝜽 (see

before) and Taylor expansion around the minimum (𝜒2(𝜽)) up to 2nd order is exact:

★ If we evaluate 𝜒2 at 𝜃i=𝜃i+𝜎i (with 𝜃j=𝜃j (∀j≠i)):

★ So 𝛥𝜒2 = 1 gives the 1𝜎 contour of our measurement. That is, the contours in parameter space whose tangents are at one standard deviation away from the LS estimates.

★ If 𝜆(x;𝜽) is NOT linear in the parameters 𝜽, the contour is not elliptical and the standard deviations cannot be obtained from the tangents. Still, it defines a region in the parameter space as a confidence region (see later).

21

✓̂ = (ATV �1A)�1ATV �1y ⌘ By (44)

E[✓̂] = E[(ATV �1A)�1ATV �1y] = (ATV �1A)�1ATV �1E[y]

= (ATV �1A)�1ATV �1A✓ = ✓

U = B V BT = (ATV �1A)�1 (45)

(U�1)ij = 1

2

 @2�2

@✓i@✓j

�

✓=✓̂

(46)

�2(✓) = �2(✓̂) + 1

2

mX

i,j=1

 @2�2

@✓i@✓j

�

✓=✓̂

(✓i � ✓̂i)(✓j � ✓̂j)

= �2(✓̂) + mX

i,j=1

(✓i � ✓̂i)(U �1)ij(✓j � ✓̂j) (47)

6

^ ^^

✓̂ = (ATV �1A)�1ATV �1y ⌘ By (44)

E[✓̂] = E[(ATV �1A)�1ATV �1y] = (ATV �1A)�1ATV �1E[y]

= (ATV �1A)�1ATV �1A✓ = ✓

U = B V BT = (ATV �1A)�1 (45)

(U�1)ij = 1

2

 @2�2

@✓i@✓j

�

✓=✓̂

(46)

�2(✓) = �2(✓̂) + 1

2

mX

i,j=1

 @2�2

@✓i@✓j

�

✓=✓̂

(✓i � ✓̂i)(✓j � ✓̂j)

= �2(✓̂) + mX

i,j=1

(✓i � ✓̂i)(U �1)ij(✓j � ✓̂j) (47)

��2 = U�1 ii

�̂i 2 =

1

�̂i 2 �̂i

2 = 1 (48)

6

^

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

𝜒2+1 contours in linear LS fit

22

100 The method of least squares

(a) x2 = X2

mln + 1 46.5

0.8 1-----1-- ············································ .. T· .. ······ ......................................... .

46

0.6

45.5 ........•.......... : ......................... . LS estimate

0.4 L-_---L"--_...l..-_--'-__ -L.l._---' 2.5 2.6 2.7 2.8 2.9 0.4 0.6 0.8 1.2 1.4

Fig. 7.3 (a) The X2 as a function of eo for the zero-order polynomial fit shown in Fig. 7.2. The horizontal lines indicate X;'in and X;'in + 1. The corresponding eo values (vertical lines) are the LS estimate 00 and 00 ± 0-eo' (b) The LS estimates Bo and [h for the first-order polynomial

fit in Fig. 7.2. The tangents to the contour X2 (Bo, B1 ) = X;'in + 1 correspond to 00 ± 0-eo and

Ol±0-81'

(;-9 0

yr;: = 0.30

(;-9 1

;0:; = 0.10

cov[eo, e1 ] = U0 1 = -0.028,

corresponding to a correlation coefficient of r = -0.90. As in the case of maxi- mum likelihood, the standard deviations correspond to the tangents of the ellipse, and the correlation coefficient to its width and angle of inclination (see equations {6.31} and (6.32)).

Since the two estimators eo and e1 have a strong negative correlation, it is important to include the covariance, or equivalently the correlation coefficient, when reporting the results of the fit. Recall from Section 1.7 that one can always define two new quantities, i}o and i}1, from the original eo and e1 by means of an orthogonal transformation such that cov[i}o, 7h] = O. However, although it is generally easier to deal with uncorrelated quantities, the transformed parameters may not have as direct an interpretation as the original ones.

7.4 Least squares with binned data In the previous examples, the function relating the 'true' values). to the variable x was not necessarily a p.d.f. for x, but an arbitrary function. It can, however, be a p.d.f., or it can be proportional to one. Suppose, for example, one has n

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Confidence intervals with LS estimators

★ (Will be extended/systematized later)

★ A confidence interval of an estimator is an interval that covers the true value of a parameter with a specified probability (called coverage or confidence level (CL))

★ If the pdf of the estimator is Gaussian, 𝜃+/-𝜎𝜃 will contain the true value 𝜃 in 68.27% of the experiments (see proof later).

★ In LS, this corresponds to the interval defined by 𝛥𝜒2=1

★ In 2D, the ellipse defined by 𝛥𝜒2=1 has a coverage of 39%

23

^ ^

39. Statistics 29

parameter estimates and their correlation, one sometimes reports the value of θi, which minimizes χ2 at a fixed value of θj , such as the PDG best value. This θi value lies along the dotted line between the points where the ellipse becomes tangent to vertical, and has statistical error σinner as shown on the figure, where σinner = (1 − ρ2

ij) 1/2σi. Instead of

the correlation ρij , one reports the dependency d!θi/dθj , which is the slope of the dotted

line. This slope is related to the correlation coefficient by d!θi/dθj = ρij × σi σj

.

θ i

φ

θ i

jσ

θj

iσ

jσ

iσ

^

θ j ^

ij iρ σ

innerσ

�

Figure 39.5: Standard error ellipse for the estimators !θi and !θj . In the case shown the correlation is negative.

As in the single-variable case, because of the symmetry of the Gaussian function between θ and !θ, one finds that contours of constant lnL or χ2 cover the true values with a certain, fixed probability. That is, the confidence region is determined by

lnL(θ) ≥ lnLmax − ∆ lnL , (39.68)

or where a χ2 has been defined for use with the method of least-squares,

χ2(θ) ≤ χ2 min + ∆χ2 . (39.69)

Values of ∆χ2 or 2∆ lnL are given in Table 39.2 for several values of the coverage probability 1 − α and number of fitted parameters m. For Gaussian distributed data, these are related by ∆χ2 = 2∆ lnL = F−1

χ2 m

(1 − α), where F−1 χ2

m is the chi-square quantile

(inverse of the cumulative distribution) for m degrees of freedom.

For non-Gaussian data samples, the probability for the regions determined by Equations (39.68) or (39.69) to cover the true value of θ becomes independent of θ only in the large-sample limit. So for a finite data sample these are not exact confidence regions according to our previous definition. Nevertheless, they can still have a coverage probability only weakly dependent on the true parameter, and approximately as given in Table 39.2. In any case, the coverage probability of the intervals or regions obtained according to this procedure can in principle be determined as a function of the true parameter(s), for example, using a Monte Carlo calculation.

One of the practical advantages of intervals that can be constructed from the log-likelihood function or χ2 is that it is relatively simple to produce the interval

October 1, 2016 19:59

30 39. Statistics

Table 39.2: Values of ∆χ2 or 2∆ lnL corresponding to a coverage probability 1 − α in the large data sample limit, for joint estimation of m parameters.

(1 − α) (%) m = 1 m = 2 m = 3

68.27 1.00 2.30 3.53

90. 2.71 4.61 6.25

95. 3.84 5.99 7.82

95.45 4.00 6.18 8.03

99. 6.63 9.21 11.34

99.73 9.00 11.83 14.16

for the combination of several experiments. If N independent measurements result in log-likelihood functions lnLi(θ), then the combined log-likelihood function is simply the sum,

lnL(θ) = N!

i=1

lnLi(θ) . (39.70)

This can then be used to determine an approximate confidence interval or region with Eq. (39.68), just as with a single experiment.

39.4.2.3. Poisson or binomial data:

Another important class of measurements consists of counting a certain number of events, n. In this section, we will assume these are all events of the desired type, i.e., there is no background. If n represents the number of events produced in a reaction with cross section σ, say, in a fixed integrated luminosity L, then it follows a Poisson distribution with mean µ = σL. If, on the other hand, one has selected a larger sample of N events and found n of them to have a particular property, then n follows a binomial distribution where the parameter p gives the probability for the event to possess the property in question. This is appropriate, e.g., for estimates of branching ratios or selection efficiencies based on a given total number of events.

For the case of Poisson distributed n, limits on the mean value µ can be found from the Neyman procedure as discussed in Section 39.4.2.1 with n used directly as the statistic x . The upper and lower limits are found to be

µlo = 1 2F−1

χ2 (αlo; 2n) , (39.71a)

µup = 1 2F−1

χ2 (1 − αup; 2(n + 1)) , (39.71b)

where confidence levels of 1 − αlo and 1 − αup, refer separately to the corresponding

intervals µ ≥ µlo and µ ≤ µup, and F−1 χ2 is the quantile of the χ2 distribution (inverse of

the cumulative distribution). The quantiles F−1 χ2 can be obtained from standard tables or

from the ROOT routine TMath::ChisquareQuantile. For central confidence intervals at confidence level 1 − α, set αlo = αup = α/2.

October 1, 2016 19:59

“1𝜎”

“2𝜎”

“3𝜎”

{

68 .2

7% C

L

{68.27% CL

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Goodness of fit of LS estimators ★ (Will be extended/systematized

in Hypothesis Testing part) ★ Goodness of fit tells us how well

the data are described by 𝜆(x;𝜽) ★ For least-squares, that can be

quantified by the probability that, with 𝜆(x;𝜽) being true, the value of 𝜒2 was as high (at least) as the one obtained, i.e:

with f(z;nd) the 𝜒2 pdf for nd degrees of freedom (nd=N-m, number of data points-number of free parameters)

24

S2 = 1

n

nX

i=1

(xi � µ)2 = x2 � µ2 (29)

1

n

nX

i=1

xi + 1

n (30)

V̂xy = 1

n� 1

nX

i=1

(xi � x̄)(yi � ȳ) = n

n� 1 (xy � x̄ȳ) (31)

rxy = V̂xy

sxsy =

P n

i=1(xi � x̄)(yi � ȳ)qP n

j=1(xj � x̄)2 P

n

k=1(yk � ȳ)2 =

xy � x̄ȳq (x2 � x̄2)(y2 � ȳ2)

(32)

x̄ = nX

i=1

wixi (33)

E[x̄] = µ ) nX

i=1

wiE[xi] = µ ) nX

i=1

wiµ = µ ) µ nX

i=1

wi = µ ) nX

i=1

wi = 1 (34)

wi = 1/�2

iP n

j=1 1/� 2 j

(35)

V [s2] = 1

n

✓ µ4 �

n� 3

n� 1 µ2 2

◆ (36)

mk = 1

n� 1

nX

i=1

(xi � x̄)k (37)

Uij ⇡ nX

k,l=1

@⌘i @✓k

@⌘j @✓l

���� ✓̂

Vkl (38)

U ⇡ AV AT (39)

Aij = @⌘i @✓j

���� ✓̂

(40)

p =

Z 1

�2 f(z;nd) dz (41)

5

^

^

20 39. Statistics

1 2 3 4 5 7 10 20 30 40 50 70 100 0.001

0.002

0.005

0.010

0.020

0.050

0.100

0.200

0.500

1.000

p -v

a lu

e fo

r te

st α

f or

c on

fi d

en ce

i n

te rv

a ls

3 42 6 8

10

15

20

25

30

40

50

n = 1

χ2

Figure 39.1: One minus the χ2 cumulative distribution, 1−F (χ2; n), for n degrees of freedom. This gives the p-value for the χ2 goodness-of-fit test as well as one minus the coverage probability for confidence regions (see Sec. 39.4.2.2).

0 10 20 30 40 50 0.0

0.5

1.0

1.5

2.0

2.5

Degrees of freedom n

50%

10%

90%

99%

95%

68%

32%

5%

1%

χ2/n

Figure 39.2: The ‘reduced’ χ2, equal to χ2/n, for n degrees of freedom. The curves show as a function of n the χ2/n that corresponds to a given p-value.

October 1, 2016 19:59

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Summary of LS ★ LS method can be applied without knowing the pdf of the

measurements y, only their covariance matrix Vij=cov[yi,yj] is needed

★ In LS the measurements are function 𝜆(x;𝜽) of a variable x known without errors and some free parameters 𝜽, to be estimated

★ When 𝜆(x;𝜽) is linear in the parameters 𝜽: ✦ the LS estimators are unbiased

✦ the LS estimator have the lowest possible variance (Rao-Crammer limit)

✦ the parameter errors are Gaussian (we can use the 𝛥𝜒2=𝜎2 rule)

✦ the value 𝜒2min (together with the number of degrees of freedom) provides a measure of the goodness of the fit

✦ (mind the difference of the latter two: estimators can have small errors and still the goodness of fit be very bad)

25

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Exercises 4. Using the error propagation formula, demonstrate that the error of y = f(x) = a + bx, is

given by σy2 = σa2 + x2σb2 + 2x cov[a,b], knowing the covariance matrix for a and b:

5. Use LS to:

a. Find the parameters and errors of the line that best fits the data below, with 𝜎=0.3

b. Plot the data and the fit for 0<x<5

c. Plot the result of the fit together with the error for x=0.5. Study the consequences of ignoring the correlation term in the covariance matrix

6. Verify the results of exercise 5c with the following numerical simulation

a. For each xi generate a random yi with a Gaussian distribution N(a+bx,𝜎2)

b. Fit the new data and predict the y for x=0.5

c. Repeat the steps above 1000 times and construct a histogram with the values of y(0.5). Obtain the standard deviation of these data and compare to the results of 5c

26

S2 = 1

n

nX

i=1

(xi � µ)2 = x2 � µ2 (29)

1

n

nX

i=1

xi + 1

n (30)

V̂xy = 1

n� 1

nX

i=1

(xi � x̄)(yi � ȳ) = n

n� 1 (xy � x̄ȳ) (31)

rxy = V̂xy

sxsy =

P n

i=1(xi � x̄)(yi � ȳ)qP n

j=1(xj � x̄)2 P

n

k=1(yk � ȳ)2 =

xy � x̄ȳq (x2 � x̄2)(y2 � ȳ2)

(32)

x̄ = nX

i=1

wixi (33)

E[x̄] = µ ) nX

i=1

wiE[xi] = µ ) nX

i=1

wiµ = µ ) µ nX

i=1

wi = µ ) nX

i=1

wi = 1 (34)

wi = 1/�2

iP n

j=1 1/� 2 j

(35)

V [s2] = 1

n

✓ µ4 �

n� 3

n� 1 µ2 2

◆ (36)

mk = 1

n� 1

nX

i=1

(xi � x̄)k (37)

Uij ⇡ nX

k,l=1

@⌘i @✓k

@⌘j @✓l

���� ✓̂

Vkl (38)

U ⇡ AV AT (39)

Aij = @⌘i @✓j

���� ✓̂

(40)

p =

Z 1

�2 f(z;nd) dz (41)

V =

✓ �2 a cov[a, b]

cov[a, b] �2 b

◆ (42)

5

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Maximum Likelihood

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Maximum Likelihood (ML) estimators ★ Consider a random variable x distributed according to a pdf f(x;𝜽), of known

functional form with some unknown parameters 𝜽 = (𝜃1,…,𝜃m) to be estimated

★ Suppose that we have n independent measurements of x: x1,…,xn.

★ If f(x;𝜽) is the true pdf, we have ✦ Probability of xi to be in interval [xi,xi+dxi] = f(xi;𝜽) dxi

✦ Since measurements are independent:

★ For the true pdf, one expects high probability for the data that were actually measured. Since the dxi do not depend on 𝜃, the same applies to the likelihood function:

i.e. the joint pdf regarded as a function of the parameters

★ In ML we find the parameters 𝜽 that maximize L(𝜽) (usually logL(𝜽), or minimize -2logL(𝜽), we’ll see why)

28

S2 = 1

n

nX

i=1

(xi � µ)2 = x2 � µ2 (29)

1

n

nX

i=1

xi + 1

n (30)

V̂xy = 1

n� 1

nX

i=1

(xi � x̄)(yi � ȳ) = n

n� 1 (xy � x̄ȳ) (31)

rxy = V̂xy

sxsy =

P n

i=1(xi � x̄)(yi � ȳ)qP n

j=1(xj � x̄)2 P

n

k=1(yk � ȳ)2 =

xy � x̄ȳq (x2 � x̄2)(y2 � ȳ2)

(32)

x̄ = nX

i=1

wixi (33)

E[x̄] = µ ) nX

i=1

wiE[xi] = µ ) nX

i=1

wiµ = µ ) µ nX

i=1

wi = µ ) nX

i=1

wi = 1 (34)

wi = 1/�2

iP n

j=1 1/� 2 j

(35)

V [s2] = 1

n

✓ µ4 �

n� 3

n� 1 µ2 2

◆ (36)

mk = 1

n� 1

nX

i=1

(xi � x̄)k (37)

Uij ⇡ nX

k,l=1

@⌘i @✓k

@⌘j @✓l

���� ✓̂

Vkl (38)

U ⇡ AV AT (39)

Aij = @⌘i @✓j

���� ✓̂

(40)

p =

Z 1

�2 f(z;nd) dz (41)

V =

✓ �2 a cov[a, b]

cov[a, b] �2 b

◆ (42)

probability that xi in [xi, xi + dxi] for all i = nY

i=1

f(xi;✓)dxi (43)

5

S2 = 1

n

nX

i=1

(xi � µ)2 = x2 � µ2 (29)

1

n

nX

i=1

xi + 1

n (30)

V̂xy = 1

n� 1

nX

i=1

(xi � x̄)(yi � ȳ) = n

n� 1 (xy � x̄ȳ) (31)

rxy = V̂xy

sxsy =

P n

i=1(xi � x̄)(yi � ȳ)qP n

j=1(xj � x̄)2 P

n

k=1(yk � ȳ)2 =

xy � x̄ȳq (x2 � x̄2)(y2 � ȳ2)

(32)

x̄ = nX

i=1

wixi (33)

E[x̄] = µ ) nX

i=1

wiE[xi] = µ ) nX

i=1

wiµ = µ ) µ nX

i=1

wi = µ ) nX

i=1

wi = 1 (34)

wi = 1/�2

iP n

j=1 1/� 2 j

(35)

V [s2] = 1

n

✓ µ4 �

n� 3

n� 1 µ2 2

◆ (36)

mk = 1

n� 1

nX

i=1

(xi � x̄)k (37)

Uij ⇡ nX

k,l=1

@⌘i @✓k

@⌘j @✓l

���� ✓̂

Vkl (38)

U ⇡ AV AT (39)

Aij = @⌘i @✓j

���� ✓̂

(40)

p =

Z 1

�2 f(z;nd) dz (41)

V =

✓ �2 a cov[a, b]

cov[a, b] �2 b

◆ (42)

probability that xi in [xi, xi + dxi] for all i = nY

i=1

f(xi;✓)dxi (43)

L(✓) = nY

i=1

f(xi;✓) (44)

5

x

f(x)

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

ML estimator for the Poisson distribution

★ We have N measurements n=(n1,…,ni) of Poisson process with parameter 𝜇 ✦ Remember both the mean and variance of Poisson distribution are 𝜇

✦ Poisson pdf:

★ Then, the joint likelihood for the parameter 𝜇 (given the measurements n) is:

★ Taking the derivative with respect to 𝜇, equaling to zero and solving for 𝜇 gives us the ML estimator 𝜇ML:

★ The ML estimator of the Poisson mean is the arithmetic mean ✦ Recall that for Poisson E[n] = 𝜇, so this results makes perfect sense

✦ Also 𝜇ML is unbiased, since arithmetic mean is always unbiased estimator of the mean

29

✓̂ = (ATV �1A)�1ATV �1y ⌘ By (44)

E[✓̂] = E[(ATV �1A)�1ATV �1y] = (ATV �1A)�1ATV �1E[y]

= (ATV �1A)�1ATV �1A✓ = ✓

U = B V BT = (ATV �1A)�1 (45)

(U�1)ij = 1

2

 @2�2

@✓i@✓j

�

✓=✓̂

(46)

�2(✓) = �2(✓̂) + 1

2

mX

i,j=1

 @2�2

@✓i@✓j

�

✓=✓̂

(✓i � ✓̂i)(✓j � ✓̂j)

= �2(✓̂) + mX

i,j=1

(✓i � ✓̂i)(U �1)ij(✓j � ✓̂j) (47)

��2 = U�1 ii

�̂i 2 =

1

�̂i 2 �̂i

2 = 1 (48)

P (n;µ) = µn

n! e�µ (49)

6

^

L(✓) = nY

i=1

f(xi;✓) (44)

�(x;✓) = ✓0 + ✓1x

(1 + ✓2x+ ✓3x2) (45)

L(µ;n) = NY

i=1

µni

ni! e�µ ) logL = �Nµ+

NX

i=1

(ni logµ� log ni!) (46)

@ logL

@µ

���� µ=µ̂ML

= 0 ) N � 1

µ̂ML

NX

i=1

ni = 0 ) µ̂ML = 1

N

NX

i=1

ni = n̄ (47)

6

^

U ⇡ AV AT (75)

Aij = @⌘i @✓j

���� ✓=µ

(76)

p =

Z 1

�2 f(z;nd) dz (77)

V =

✓ �2 a cov[a, b]

cov[a, b] �2 b

◆ (78)

probability that xi in [xi, xi + dxi] for all i = nY

i=1

f(xi;✓)dxi (79)

L(✓) = nY

i=1

f(xi;✓) (80)

�(x;✓) = ✓0 + ✓1x

(1 + ✓2x+ ✓3x2) (81)

L(µ;n) = NY

i=1

µni

ni! e�µ ) logL = �Nµ+

NX

i=1

(ni log µ� log ni!) (82)

@ logL

@µ

���� µ=µ̂ML

= 0 ) N � 1

µ̂ML

NX

i=1

ni = 0 ) µ̂ML = 1

N

NX

i=1

ni = n̄ (83)

f(t; ⌧) = 1

⌧ e�

t ⌧ (84)

logL(⌧) = NX

i=1

✓ log

1

⌧ � ti

⌧

◆ (85)

9

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

ML estimator for exponential distribution

★ Suppose we have data t={t1,…,tn} with exponential pdf:

★ The logL function is:

★ We maximize it and solve it for 𝜏 to obtain the ML estimator

which again is t, i.e. the arithmetic mean of the ti ✦ therefore 𝜏ML = t is an unbiased estimator of the mean of the

exponential distribution (𝜏, or equivalently E[t])

30

L(✓) = nY

i=1

f(xi;✓) (44)

�(x;✓) = ✓0 + ✓1x

(1 + ✓2x+ ✓3x2) (45)

L(µ;n) = NY

i=1

µni

ni! e�µ ) logL = �Nµ+

NX

i=1

(ni logµ� log ni!) (46)

@ logL

@µ

���� µ=µ̂ML

= 0 ) N � 1

µ̂ML

NX

i=1

ni = 0 ) µ̂ML = 1

N

NX

i=1

ni = n̄ (47)

f(t; ⌧) = 1

⌧ e�

t ⌧ (48)

6

@ logL

@⌧

���� ⌧=⌧̂ML

= 0 ) nX

i=1

✓ � 1

⌧̂ML +

ti ⌧̂2ML

◆ = 0 ) ⌧̂ML =

1

n

nX

i=1

ti = t̄ (1)

⌧̂ML = 1

n

nX

i=1

ti (2)

Finally, we compute the radius ⌅ of the MAGIC e↵ective field of view, defined such that observations of an isotropic gamma-ray flux with a hypothetical instrument with a flat-top acceptance R(⇠) = R(0) for ⇠ < ⌅, and R(⇠) = 0 for ⇠ > ⌅, would yield the same number of detected gamma rays as with MAGIC, when no cuts on the arrival direction are applied. We can therefore obtain ⌅ from the condition

R ⌅ 0 2⇡ ⇠R(0) d⇠ =

R1 0 2⇡ ⇠R(⇠) d⇠, where R(⇠) is shown

in Figure 19-bottom, yielding ⌅ = 1�. We note, however, that standard observations of sources with an extension larger than 0.4� are technically di�cult, as in that case the edge of the source would fall into the background estimation region. Nevertheless, the e↵ective field of view is a useful quantity for non-standard observations of di↵use signals like, e.g. the cosmic electron flux [Borla-Tridon ICRC, HESS electron spectrum paper].

The di↵erences between our results for VERITAS obtained using the conventional likelihood (Eq. 2.1, Figure 8a) and those published by the Collaboration [16] (shown in Figures 8b-d) are due to the di↵erent Aeff assumed, since all the remaining values relevant for the computation of the exclusion contours (TOBS , J and upper limit to the number of signal events) are taken from Ref. [16]. In this work, we have used the values of Aeff reported by Wood [ref] for Segue observation/analysis. In addition, we have checked that the Improvement Factors that we obtain do not di↵er significantly (< 5%) if we use the values of Aeff reported by McCutcheon [ref] (assumed for the analysis of the globular cluster W13). From this, we infer the validity of the obtained Improvement Factors also for the Aeff actually used in Ref. [16].

In addition, according to these results, the sensitivity gain of the CTA with respect to VERITAS would be marginal, or even nonexistent, for certain annihilation channels and mass ranges. We have traced this inconsistency down to a probable overestimation of the VERI- TAS performance assumed in [16]. For that, we have used the response functions assumed for VERITAS, MAGIC and CTA to compute the integral sensitivity (5� significance in 50 hours of observations) for a Crab-like spectrum1 at the analysis threshold, for the di↵erent in- struments. The results obtained for MAGIC (1.3% of Crab flux above 110 GeV) and CTA (0.30% of Crab flux above 75 GeV) are consistent with those published by the respective collaborations (Refs. [MAGIC] and [CTA] respectively). On the other hand, VERITAS re- sults imply a sensitivity of 0.32% of Crab flux above 165 GeV, more than a factor 2 bet- ter than what reported by the collaboration [ref: http://veritas.sao.arizona.edu/about-veritas- mainmenu-81/veritas-specifications-mainmenu-111]. It must also be stressed that this result is obtained assuming the e↵ective areas as reported by Wood [ref], which for most of the decay channels and mass ranges produce slightly less constraining results than what published by VERITAS [16]. Therefore we expect that repeating this calculation with the Aeff actually used in [16] would produce an even larger disagreement between the implied and reference sensitivities to the Crab Nebula spectrum.

d�

dE (⌦) =

d�PP

dE ⇥ J(⌦) (3)

1dN/dE = 5.8�13(E/300GeV)�2.32�0.13 log10(E/300GeV) GeV�1 cm�2 s�1 [ref MAGIC performance]

1

-

-^

U ⇡ AV AT (75)

Aij = @⌘i @✓j

���� ✓=µ

(76)

p =

Z 1

�2 f(z;nd) dz (77)

V =

✓ �2 a cov[a, b]

cov[a, b] �2 b

◆ (78)

probability that xi in [xi, xi + dxi] for all i = nY

i=1

f(xi;✓)dxi (79)

L(✓) = nY

i=1

f(xi;✓) (80)

�(x;✓) = ✓0 + ✓1x

(1 + ✓2x+ ✓3x2) (81)

L(µ;n) = NY

i=1

µni

ni! e�µ ) logL = �Nµ+

NX

i=1

(ni log µ� log ni!) (82)

@ logL

@µ

���� µ=µ̂ML

= 0 ) N � 1

µ̂ML

NX

i=1

ni = 0 ) µ̂ML = 1

N

NX

i=1

ni = n̄ (83)

f(t; ⌧) = 1

⌧ e�

t ⌧ (84)

logL(⌧) = nX

i=1

✓ log

1

⌧ � ti

⌧

◆ (85)

9

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

ML estimator of a function of a parameter

★ Following with the exponential distribution, now let us assume that we are interested in estimating 𝜆=1/𝜏

★ In general given a function a(𝜃) of a parameter 𝜃, one has

thus ∂L/∂𝜃=0 at 𝜃 implies ∂L/∂a=0 at a=a(𝜃) (unless ∂a/∂𝜃=0, i.e. if a does not really depend on 𝜃)

★ Therefore, the estimator for 𝜆 is:

★ It can be shown that the bias is given by [✪]:

i.e. even if 𝜏ML is unbiased 𝜆ML is only asymptotically unbiased. This is a general property and not only for the exponential distribution

31

@L

@✓ =

@L

@a

@a

@✓ = 0 (1)

@ logL

@⌧

���� ⌧=⌧̂ML

= 0 ) nX

i=1

✓ � 1

⌧̂ML +

ti ⌧̂2ML

◆ = 0 ) ⌧̂ML =

1

n

nX

i=1

ti = t̄ (2)

⌧̂ML = 1

n

nX

i=1

ti (3)

Finally, we compute the radius ⌅ of the MAGIC e↵ective field of view, defined such that observations of an isotropic gamma-ray flux with a hypothetical instrument with a flat-top acceptance R(⇠) = R(0) for ⇠ < ⌅, and R(⇠) = 0 for ⇠ > ⌅, would yield the same number of detected gamma rays as with MAGIC, when no cuts on the arrival direction are applied. We can therefore obtain ⌅ from the condition

R ⌅ 0 2⇡ ⇠R(0) d⇠ =

R1 0 2⇡ ⇠R(⇠) d⇠, where R(⇠) is shown

in Figure 19-bottom, yielding ⌅ = 1�. We note, however, that standard observations of sources with an extension larger than 0.4� are technically di�cult, as in that case the edge of the source would fall into the background estimation region. Nevertheless, the e↵ective field of view is a useful quantity for non-standard observations of di↵use signals like, e.g. the cosmic electron flux [Borla-Tridon ICRC, HESS electron spectrum paper].

The di↵erences between our results for VERITAS obtained using the conventional likelihood (Eq. 2.1, Figure 8a) and those published by the Collaboration [16] (shown in Figures 8b-d) are due to the di↵erent Aeff assumed, since all the remaining values relevant for the computation of the exclusion contours (TOBS , J and upper limit to the number of signal events) are taken from Ref. [16]. In this work, we have used the values of Aeff reported by Wood [ref] for Segue observation/analysis. In addition, we have checked that the Improvement Factors that we obtain do not di↵er significantly (< 5%) if we use the values of Aeff reported by McCutcheon [ref] (assumed for the analysis of the globular cluster W13). From this, we infer the validity of the obtained Improvement Factors also for the Aeff actually used in Ref. [16].

In addition, according to these results, the sensitivity gain of the CTA with respect to VERITAS would be marginal, or even nonexistent, for certain annihilation channels and mass ranges. We have traced this inconsistency down to a probable overestimation of the VERI- TAS performance assumed in [16]. For that, we have used the response functions assumed for VERITAS, MAGIC and CTA to compute the integral sensitivity (5� significance in 50 hours of observations) for a Crab-like spectrum1 at the analysis threshold, for the di↵erent in- struments. The results obtained for MAGIC (1.3% of Crab flux above 110 GeV) and CTA (0.30% of Crab flux above 75 GeV) are consistent with those published by the respective collaborations (Refs. [MAGIC] and [CTA] respectively). On the other hand, VERITAS re- sults imply a sensitivity of 0.32% of Crab flux above 165 GeV, more than a factor 2 bet- ter than what reported by the collaboration [ref: http://veritas.sao.arizona.edu/about-veritas- mainmenu-81/veritas-specifications-mainmenu-111]. It must also be stressed that this result is obtained assuming the e↵ective areas as reported by Wood [ref], which for most of the decay channels and mass ranges produce slightly less constraining results than what published by VERITAS [16]. Therefore we expect that repeating this calculation with the Aeff actually used in [16] would produce an even larger disagreement between the implied and reference sensitivities to the Crab Nebula spectrum.

1dN/dE = 5.8�13(E/300GeV)�2.32�0.13 log10(E/300GeV) GeV�1 cm�2 s�1 [ref MAGIC performance]

1

^ ^^

�̂ML = 1

⌧̂ML =

nP m

i=1 ti (1)

@L

@✓ =

@L

@a

@a

@✓ = 0 (2)

@ logL

@⌧

���� ⌧=⌧̂ML

= 0 ) nX

i=1

✓ � 1

⌧̂ML +

ti ⌧̂2ML

◆ = 0 ) ⌧̂ML =

1

n

nX

i=1

ti = t̄ (3)

⌧̂ML = 1

n

nX

i=1

ti (4)

Finally, we compute the radius ⌅ of the MAGIC e↵ective field of view, defined such that observations of an isotropic gamma-ray flux with a hypothetical instrument with a flat-top acceptance R(⇠) = R(0) for ⇠ < ⌅, and R(⇠) = 0 for ⇠ > ⌅, would yield the same number of detected gamma rays as with MAGIC, when no cuts on the arrival direction are applied. We can therefore obtain ⌅ from the condition

R ⌅ 0 2⇡ ⇠R(0) d⇠ =

R1 0 2⇡ ⇠R(⇠) d⇠, where R(⇠) is shown

in Figure 19-bottom, yielding ⌅ = 1�. We note, however, that standard observations of sources with an extension larger than 0.4� are technically di�cult, as in that case the edge of the source would fall into the background estimation region. Nevertheless, the e↵ective field of view is a useful quantity for non-standard observations of di↵use signals like, e.g. the cosmic electron flux [Borla-Tridon ICRC, HESS electron spectrum paper].

The di↵erences between our results for VERITAS obtained using the conventional likelihood (Eq. 2.1, Figure 8a) and those published by the Collaboration [16] (shown in Figures 8b-d) are due to the di↵erent Aeff assumed, since all the remaining values relevant for the computation of the exclusion contours (TOBS , J and upper limit to the number of signal events) are taken from Ref. [16]. In this work, we have used the values of Aeff reported by Wood [ref] for Segue observation/analysis. In addition, we have checked that the Improvement Factors that we obtain do not di↵er significantly (< 5%) if we use the values of Aeff reported by McCutcheon [ref] (assumed for the analysis of the globular cluster W13). From this, we infer the validity of the obtained Improvement Factors also for the Aeff actually used in Ref. [16].

In addition, according to these results, the sensitivity gain of the CTA with respect to VERITAS would be marginal, or even nonexistent, for certain annihilation channels and mass ranges. We have traced this inconsistency down to a probable overestimation of the VERI- TAS performance assumed in [16]. For that, we have used the response functions assumed for VERITAS, MAGIC and CTA to compute the integral sensitivity (5� significance in 50 hours of observations) for a Crab-like spectrum1 at the analysis threshold, for the di↵erent in- struments. The results obtained for MAGIC (1.3% of Crab flux above 110 GeV) and CTA (0.30% of Crab flux above 75 GeV) are consistent with those published by the respective collaborations (Refs. [MAGIC] and [CTA] respectively). On the other hand, VERITAS re- sults imply a sensitivity of 0.32% of Crab flux above 165 GeV, more than a factor 2 bet- ter than what reported by the collaboration [ref: http://veritas.sao.arizona.edu/about-veritas- mainmenu-81/veritas-specifications-mainmenu-111]. It must also be stressed that this result is obtained assuming the e↵ective areas as reported by Wood [ref], which for most of the decay channels and mass ranges produce slightly less constraining results than what published by VERITAS [16]. Therefore we expect that repeating this calculation with the Aeff actually used

1dN/dE = 5.8�13(E/300GeV)�2.32�0.13 log10(E/300GeV) GeV�1 cm�2 s�1 [ref MAGIC performance]

1

b[�̂ML] = E[�̂ML]� � = � n

n� 1 � � =

1

n� 1 � (1)

�̂ML = 1

⌧̂ML =

nP m

i=1 ti (2)

@L

@✓ =

@L

@a

@a

@✓ = 0 (3)

@ logL

@⌧

���� ⌧=⌧̂ML

= 0 ) nX

i=1

✓ � 1

⌧̂ML +

ti ⌧̂2ML

◆ = 0 ) ⌧̂ML =

1

n

nX

i=1

ti = t̄ (4)

⌧̂ML = 1

n

nX

i=1

ti (5)

Finally, we compute the radius ⌅ of the MAGIC e↵ective field of view, defined such that observations of an isotropic gamma-ray flux with a hypothetical instrument with a flat-top acceptance R(⇠) = R(0) for ⇠ < ⌅, and R(⇠) = 0 for ⇠ > ⌅, would yield the same number of detected gamma rays as with MAGIC, when no cuts on the arrival direction are applied. We can therefore obtain ⌅ from the condition

R ⌅ 0 2⇡ ⇠R(0) d⇠ =

R1 0 2⇡ ⇠R(⇠) d⇠, where R(⇠) is shown

in Figure 19-bottom, yielding ⌅ = 1�. We note, however, that standard observations of sources with an extension larger than 0.4� are technically di�cult, as in that case the edge of the source would fall into the background estimation region. Nevertheless, the e↵ective field of view is a useful quantity for non-standard observations of di↵use signals like, e.g. the cosmic electron flux [Borla-Tridon ICRC, HESS electron spectrum paper].

The di↵erences between our results for VERITAS obtained using the conventional likelihood (Eq. 2.1, Figure 8a) and those published by the Collaboration [16] (shown in Figures 8b-d) are due to the di↵erent Aeff assumed, since all the remaining values relevant for the computation of the exclusion contours (TOBS , J and upper limit to the number of signal events) are taken from Ref. [16]. In this work, we have used the values of Aeff reported by Wood [ref] for Segue observation/analysis. In addition, we have checked that the Improvement Factors that we obtain do not di↵er significantly (< 5%) if we use the values of Aeff reported by McCutcheon [ref] (assumed for the analysis of the globular cluster W13). From this, we infer the validity of the obtained Improvement Factors also for the Aeff actually used in Ref. [16].

In addition, according to these results, the sensitivity gain of the CTA with respect to VERITAS would be marginal, or even nonexistent, for certain annihilation channels and mass ranges. We have traced this inconsistency down to a probable overestimation of the VERI- TAS performance assumed in [16]. For that, we have used the response functions assumed for VERITAS, MAGIC and CTA to compute the integral sensitivity (5� significance in 50 hours of observations) for a Crab-like spectrum1 at the analysis threshold, for the di↵erent in- struments. The results obtained for MAGIC (1.3% of Crab flux above 110 GeV) and CTA (0.30% of Crab flux above 75 GeV) are consistent with those published by the respective collaborations (Refs. [MAGIC] and [CTA] respectively). On the other hand, VERITAS re- sults imply a sensitivity of 0.32% of Crab flux above 165 GeV, more than a factor 2 bet- ter than what reported by the collaboration [ref: http://veritas.sao.arizona.edu/about-veritas- mainmenu-81/veritas-specifications-mainmenu-111]. It must also be stressed that this result is obtained assuming the e↵ective areas as reported by Wood [ref], which for most of the decay

1dN/dE = 5.8�13(E/300GeV)�2.32�0.13 log10(E/300GeV) GeV�1 cm�2 s�1 [ref MAGIC performance]

1

^ ^

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

ML estimators of Gaussian mean ★ Suppose we have n measurements of a random variable x distributed according

to Gaussian pdf of unknown 𝜇 and 𝜎:

★ The log-likelihood function is given by

✦ The last term is -1/2 𝜒2 from LS estimators ⟹ -2logL is distributed as a 𝜒2 close to µML for Gaussian pdf

✦ We can use this to define confidence intervals for µML (e.g. for 1D: 𝛥-2logL = 1 for 1σ interval), see more later

★ Setting the derivative of logL with respect to 𝜇 to zero and solving for 𝜇 yields:

i.e. the arithmetic mean, therefore µML is unbiased estimator of Gaussian mean 𝜇

32

f(x;µ,�) = 1p 2⇡�

e �(x�µ)2

�2 (1)

logL(µ,�) = nX

i=1

log f(xi;µ,�) = n log 1p 2⇡

+ n

2 log

1

�2 � 1

2

nX

i=1

(xi � µ)2

�2 (2)

b[�̂ML] = E[�̂ML]� � = � n

n� 1 � � =

1

n� 1 � (3)

�̂ML = 1

⌧̂ML =

nP m

i=1 ti (4)

@L

@✓ =

@L

@a

@a

@✓ = 0 (5)

@ logL

@⌧

���� ⌧=⌧̂ML

= 0 ) nX

i=1

✓ � 1

⌧̂ML +

ti ⌧̂2ML

◆ = 0 ) ⌧̂ML =

1

n

nX

i=1

ti = t̄ (6)

⌧̂ML = 1

n

nX

i=1

ti (7)

Finally, we compute the radius ⌅ of the MAGIC e↵ective field of view, defined such that observations of an isotropic gamma-ray flux with a hypothetical instrument with a flat-top acceptance R(⇠) = R(0) for ⇠ < ⌅, and R(⇠) = 0 for ⇠ > ⌅, would yield the same number of detected gamma rays as with MAGIC, when no cuts on the arrival direction are applied. We can therefore obtain ⌅ from the condition

R ⌅ 0 2⇡ ⇠R(0) d⇠ =

R1 0 2⇡ ⇠R(⇠) d⇠, where R(⇠) is shown

in Figure 19-bottom, yielding ⌅ = 1�. We note, however, that standard observations of sources with an extension larger than 0.4� are technically di�cult, as in that case the edge of the source would fall into the background estimation region. Nevertheless, the e↵ective field of view is a useful quantity for non-standard observations of di↵use signals like, e.g. the cosmic electron flux [Borla-Tridon ICRC, HESS electron spectrum paper].

The di↵erences between our results for VERITAS obtained using the conventional likelihood (Eq. 2.1, Figure 8a) and those published by the Collaboration [16] (shown in Figures 8b-d) are due to the di↵erent Aeff assumed, since all the remaining values relevant for the computation of the exclusion contours (TOBS , J and upper limit to the number of signal events) are taken from Ref. [16]. In this work, we have used the values of Aeff reported by Wood [ref] for Segue observation/analysis. In addition, we have checked that the Improvement Factors that we obtain do not di↵er significantly (< 5%) if we use the values of Aeff reported by McCutcheon [ref] (assumed for the analysis of the globular cluster W13). From this, we infer the validity of the obtained Improvement Factors also for the Aeff actually used in Ref. [16].

In addition, according to these results, the sensitivity gain of the CTA with respect to VERITAS would be marginal, or even nonexistent, for certain annihilation channels and mass ranges. We have traced this inconsistency down to a probable overestimation of the VERI- TAS performance assumed in [16]. For that, we have used the response functions assumed for VERITAS, MAGIC and CTA to compute the integral sensitivity (5� significance in 50 hours of observations) for a Crab-like spectrum1 at the analysis threshold, for the di↵erent in- struments. The results obtained for MAGIC (1.3% of Crab flux above 110 GeV) and CTA

1dN/dE = 5.8�13(E/300GeV)�2.32�0.13 log10(E/300GeV) GeV�1 cm�2 s�1 [ref MAGIC performance]

1

^

^

^

@ logL

@µ

���� µ=µ̂ML

= 0 ) 1

2

nX

i=1

2(xi � µ̂ML)

�2 = 0 )

nX

i=1

xi = nµ̂ML ) µ̂ML = 1

n

nX

i=1

xi (1)

µ̂ML = 1

n

nX

i=1

xi (2)

f(x;µ,�) = 1p 2⇡�

e �(x�µ)2

�2 (3)

logL(µ,�) = nX

i=1

log f(xi;µ,�) = n log 1p 2⇡

+ n

2 log

1

�2 � 1

2

nX

i=1

(xi � µ)2

�2 (4)

b[�̂ML] = E[�̂ML]� � = � n

n� 1 � � =

1

n� 1 � (5)

�̂ML = 1

⌧̂ML =

nP m

i=1 ti (6)

@L

@✓ =

@L

@a

@a

@✓ = 0 (7)

@ logL

@⌧

���� ⌧=⌧̂ML

= 0 ) nX

i=1

✓ � 1

⌧̂ML +

ti ⌧̂2ML

◆ = 0 ) ⌧̂ML =

1

n

nX

i=1

ti = t̄ (8)

⌧̂ML = 1

n

nX

i=1

ti (9)

Finally, we compute the radius ⌅ of the MAGIC e↵ective field of view, defined such that observations of an isotropic gamma-ray flux with a hypothetical instrument with a flat-top acceptance R(⇠) = R(0) for ⇠ < ⌅, and R(⇠) = 0 for ⇠ > ⌅, would yield the same number of detected gamma rays as with MAGIC, when no cuts on the arrival direction are applied. We can therefore obtain ⌅ from the condition

R ⌅ 0 2⇡ ⇠R(0) d⇠ =

R1 0 2⇡ ⇠R(⇠) d⇠, where R(⇠) is shown

in Figure 19-bottom, yielding ⌅ = 1�. We note, however, that standard observations of sources with an extension larger than 0.4� are technically di�cult, as in that case the edge of the source would fall into the background estimation region. Nevertheless, the e↵ective field of view is a useful quantity for non-standard observations of di↵use signals like, e.g. the cosmic electron flux [Borla-Tridon ICRC, HESS electron spectrum paper].

The di↵erences between our results for VERITAS obtained using the conventional likelihood (Eq. 2.1, Figure 8a) and those published by the Collaboration [16] (shown in Figures 8b-d) are due to the di↵erent Aeff assumed, since all the remaining values relevant for the computation of the exclusion contours (TOBS , J and upper limit to the number of signal events) are taken from Ref. [16]. In this work, we have used the values of Aeff reported by Wood [ref] for Segue observation/analysis. In addition, we have checked that the Improvement Factors that we obtain do not di↵er significantly (< 5%) if we use the values of Aeff reported by McCutcheon [ref] (assumed for the analysis of the globular cluster W13). From this, we infer the validity of the obtained Improvement Factors also for the Aeff actually used in Ref. [16].

1

U ⇡ AV AT (75)

Aij = @⌘i @✓j

���� ✓=µ

(76)

p =

Z 1

�2 f(z;nd) dz (77)

V =

✓ �2 a cov[a, b]

cov[a, b] �2 b

◆ (78)

probability that xi in [xi, xi + dxi] for all i = nY

i=1

f(xi;✓)dxi (79)

L(✓) = nY

i=1

f(xi;✓) (80)

�(x;✓) = ✓0 + ✓1x

(1 + ✓2x+ ✓3x2) (81)

L(µ;n) = NY

i=1

µni

ni! e�µ ) logL = �Nµ+

NX

i=1

(ni log µ� log ni!) (82)

@ logL

@µ

���� µ=µ̂ML

= 0 ) N � 1

µ̂ML

NX

i=1

ni = 0 ) µ̂ML = 1

N

NX

i=1

ni = n̄ (83)

f(t; ⌧) = 1

⌧ e�

t ⌧ (84)

logL(⌧) = nX

i=1

✓ log

1

⌧ � ti

⌧

◆ (85)

f(x;µ,�) = 1p 2⇡�

e �(x�µ)2

2�2 (86)

9

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

ML estimator of Gaussian variance ★ Setting the derivative of logL with respect to 𝜎2 to zero and solving

for 𝜎2 yields:

★ We already know this is a biased estimator with bias:

(i.e. it is asymptotically unbiased)

33

@ logL

@�2

���� �2=c�2

ML

= 0 ) nX

i=0

�1

2

1 c�2

ML

+ (xi � µ)2

2(c�2 ML)2

! = 0 ) c�2

ML = 1

n

nX

i=1

(xi � µ)2 (1)

@ logL

@µ

���� µ=µ̂ML

= 0 ) 1

2

nX

i=1

2(xi � µ̂ML)

�2 = 0 )

nX

i=1

xi = nµ̂ML ) µ̂ML = 1

n

nX

i=1

xi (2)

µ̂ML = 1

n

nX

i=1

xi (3)

f(x;µ,�) = 1p 2⇡�

e �(x�µ)2

�2 (4)

logL(µ,�) = nX

i=1

log f(xi;µ,�) = n log 1p 2⇡

+ n

2 log

1

�2 � 1

2

nX

i=1

(xi � µ)2

�2 (5)

b[�̂ML] = E[�̂ML]� � = � n

n� 1 � � =

1

n� 1 � (6)

�̂ML = 1

⌧̂ML =

nP m

i=1 ti (7)

@L

@✓ =

@L

@a

@a

@✓ = 0 (8)

@ logL

@⌧

���� ⌧=⌧̂ML

= 0 ) nX

i=1

✓ � 1

⌧̂ML +

ti ⌧̂2ML

◆ = 0 ) ⌧̂ML =

1

n

nX

i=1

ti = t̄ (9)

⌧̂ML = 1

n

nX

i=1

ti (10)

Finally, we compute the radius ⌅ of the MAGIC e↵ective field of view, defined such that observations of an isotropic gamma-ray flux with a hypothetical instrument with a flat-top acceptance R(⇠) = R(0) for ⇠ < ⌅, and R(⇠) = 0 for ⇠ > ⌅, would yield the same number of detected gamma rays as with MAGIC, when no cuts on the arrival direction are applied. We can therefore obtain ⌅ from the condition

R ⌅ 0 2⇡ ⇠R(0) d⇠ =

R1 0 2⇡ ⇠R(⇠) d⇠, where R(⇠) is shown

in Figure 19-bottom, yielding ⌅ = 1�. We note, however, that standard observations of sources with an extension larger than 0.4� are technically di�cult, as in that case the edge of the source would fall into the background estimation region. Nevertheless, the e↵ective field of view is a useful quantity for non-standard observations of di↵use signals like, e.g. the cosmic electron flux [Borla-Tridon ICRC, HESS electron spectrum paper].

The di↵erences between our results for VERITAS obtained using the conventional likelihood (Eq. 2.1, Figure 8a) and those published by the Collaboration [16] (shown in Figures 8b-d) are due to the di↵erent Aeff assumed, since all the remaining values relevant for the computation of the exclusion contours (TOBS , J and upper limit to the number of signal events) are taken from Ref. [16]. In this work, we have used the values of Aeff reported by Wood [ref] for Segue

1

b[c�2 ML] = E[c�2

ML]� �2 = n� 1

n �2 � �2 =

�1

n �2 (1)

@ logL

@�2

���� �2=c�2

ML

= 0 ) nX

i=0

�1

2

1 c�2

ML

+ (xi � µ)2

2(c�2 ML)2

! = 0 ) c�2

ML = 1

n

nX

i=1

(xi � µ)2 (2)

@ logL

@µ

���� µ=µ̂ML

= 0 ) 1

2

nX

i=1

2(xi � µ̂ML)

�2 = 0 )

nX

i=1

xi = nµ̂ML ) µ̂ML = 1

n

nX

i=1

xi (3)

µ̂ML = 1

n

nX

i=1

xi (4)

f(x;µ,�) = 1p 2⇡�

e �(x�µ)2

�2 (5)

logL(µ,�) = nX

i=1

log f(xi;µ,�) = n log 1p 2⇡

+ n

2 log

1

�2 � 1

2

nX

i=1

(xi � µ)2

�2 (6)

b[�̂ML] = E[�̂ML]� � = � n

n� 1 � � =

1

n� 1 � (7)

�̂ML = 1

⌧̂ML =

nP m

i=1 ti (8)

@L

@✓ =

@L

@a

@a

@✓ = 0 (9)

@ logL

@⌧

���� ⌧=⌧̂ML

= 0 ) nX

i=1

✓ � 1

⌧̂ML +

ti ⌧̂2ML

◆ = 0 ) ⌧̂ML =

1

n

nX

i=1

ti = t̄ (10)

⌧̂ML = 1

n

nX

i=1

ti (11)

Finally, we compute the radius ⌅ of the MAGIC e↵ective field of view, defined such that observations of an isotropic gamma-ray flux with a hypothetical instrument with a flat-top acceptance R(⇠) = R(0) for ⇠ < ⌅, and R(⇠) = 0 for ⇠ > ⌅, would yield the same number of detected gamma rays as with MAGIC, when no cuts on the arrival direction are applied. We can therefore obtain ⌅ from the condition

R ⌅ 0 2⇡ ⇠R(0) d⇠ =

R1 0 2⇡ ⇠R(⇠) d⇠, where R(⇠) is shown

in Figure 19-bottom, yielding ⌅ = 1�. We note, however, that standard observations of sources with an extension larger than 0.4� are technically di�cult, as in that case the edge of the source would fall into the background estimation region. Nevertheless, the e↵ective field of view is a useful quantity for non-standard observations of di↵use signals like, e.g. the cosmic electron flux [Borla-Tridon ICRC, HESS electron spectrum paper].

The di↵erences between our results for VERITAS obtained using the conventional likelihood (Eq. 2.1, Figure 8a) and those published by the Collaboration [16] (shown in Figures 8b-d) are

1

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

ML estimator of Gaussians with different variances

★ Let us assume the measurements x=(x1,…xn) are distributed as Gaussian all with same unknown mean 𝜇 and different known variances 𝜎i2 ✦ this is the same problem we solved with the weighted mean, i.e: we have several experiments

measuring the same parameter with different error, and want to get the global result

★ The logL function is:

and the estimator of the mean

i.e, the weighted mean

★ So far: ✦ ML seems to generalize previous results: estimators, LS, weighted mean…

✦ It provides an unified approach

✦ Estimators may be biased, but asymptotically unbiased

34

logL(µ, {�i}i=1,...n) = n log 1p 2⇡

+ n

2 log

1

�2 i

� 1

2

nX

i=1

(xi � µ)2

�2 i

(1)

b[c�2 ML] = E[c�2

ML]� �2 = n� 1

n �2 � �2 =

�1

n �2 (2)

@ logL

@�2

���� �2=c�2

ML

= 0 ) nX

i=0

�1

2

1 c�2

ML

+ (xi � µ)2

2(c�2 ML)2

! = 0 ) c�2

ML = 1

n

nX

i=1

(xi � µ)2 (3)

@ logL

@µ

���� µ=µ̂ML

= 0 ) 1

2

nX

i=1

2(xi � µ̂ML)

�2 = 0 )

nX

i=1

xi = nµ̂ML ) µ̂ML = 1

n

nX

i=1

xi (4)

µ̂ML = 1

n

nX

i=1

xi (5)

f(x;µ,�) = 1p 2⇡�

e �(x�µ)2

�2 (6)

logL(µ,�) = nX

i=1

log f(xi;µ,�) = n log 1p 2⇡

+ n

2 log

1

�2 � 1

2

nX

i=1

(xi � µ)2

�2 (7)

b[�̂ML] = E[�̂ML]� � = � n

n� 1 � � =

1

n� 1 � (8)

�̂ML = 1

⌧̂ML =

nP m

i=1 ti (9)

@L

@✓ =

@L

@a

@a

@✓ = 0 (10)

@ logL

@⌧

���� ⌧=⌧̂ML

= 0 ) nX

i=1

✓ � 1

⌧̂ML +

ti ⌧̂2ML

◆ = 0 ) ⌧̂ML =

1

n

nX

i=1

ti = t̄ (11)

⌧̂ML = 1

n

nX

i=1

ti (12)

Finally, we compute the radius ⌅ of the MAGIC e↵ective field of view, defined such that observations of an isotropic gamma-ray flux with a hypothetical instrument with a flat-top acceptance R(⇠) = R(0) for ⇠ < ⌅, and R(⇠) = 0 for ⇠ > ⌅, would yield the same number of detected gamma rays as with MAGIC, when no cuts on the arrival direction are applied. We can therefore obtain ⌅ from the condition

R ⌅ 0 2⇡ ⇠R(0) d⇠ =

R1 0 2⇡ ⇠R(⇠) d⇠, where R(⇠) is shown

in Figure 19-bottom, yielding ⌅ = 1�. We note, however, that standard observations of sources with an extension larger than 0.4� are technically di�cult, as in that case the edge of the source would fall into the background estimation region. Nevertheless, the e↵ective field of view is a

1

@ logL

@µ

���� µ=µ̂ML

= 0 ) � nX

i=1

(xi � µ̂ML)

�2 i

= 0 ) µ̂ML =

P n

i=1 xi/� 2 iP

n

i=1 1/� 2 i

(1)

logL(µ, {�i}i=1,...n) = n log 1p 2⇡

+ n

2 log

1

�2 i

� 1

2

nX

i=1

(xi � µ)2

�2 i

(2)

b[c�2 ML] = E[c�2

ML]� �2 = n� 1

n �2 � �2 =

�1

n �2 (3)

@ logL

@�2

���� �2=c�2

ML

= 0 ) nX

i=0

�1

2

1 c�2

ML

+ (xi � µ)2

2(c�2 ML)2

! = 0 ) c�2

ML = 1

n

nX

i=1

(xi � µ)2 (4)

@ logL

@µ

���� µ=µ̂ML

= 0 ) 1

2

nX

i=1

2(xi � µ̂ML)

�2 = 0 )

nX

i=1

xi = nµ̂ML ) µ̂ML = 1

n

nX

i=1

xi (5)

µ̂ML = 1

n

nX

i=1

xi (6)

f(x;µ,�) = 1p 2⇡�

e �(x�µ)2

�2 (7)

logL(µ,�) = nX

i=1

log f(xi;µ,�) = n log 1p 2⇡

+ n

2 log

1

�2 � 1

2

nX

i=1

(xi � µ)2

�2 (8)

b[�̂ML] = E[�̂ML]� � = � n

n� 1 � � =

1

n� 1 � (9)

�̂ML = 1

⌧̂ML =

nP m

i=1 ti (10)

@L

@✓ =

@L

@a

@a

@✓ = 0 (11)

@ logL

@⌧

���� ⌧=⌧̂ML

= 0 ) nX

i=1

✓ � 1

⌧̂ML +

ti ⌧̂2ML

◆ = 0 ) ⌧̂ML =

1

n

nX

i=1

ti = t̄ (12)

⌧̂ML = 1

n

nX

i=1

ti (13)

Finally, we compute the radius ⌅ of the MAGIC e↵ective field of view, defined such that observations of an isotropic gamma-ray flux with a hypothetical instrument with a flat-top acceptance R(⇠) = R(0) for ⇠ < ⌅, and R(⇠) = 0 for ⇠ > ⌅, would yield the same number of detected gamma rays as with MAGIC, when no cuts on the arrival direction are applied. We can therefore obtain ⌅ from the condition

R ⌅ 0 2⇡ ⇠R(0) d⇠ =

R1 0 2⇡ ⇠R(⇠) d⇠, where R(⇠) is shown

in Figure 19-bottom, yielding ⌅ = 1�. We note, however, that standard observations of sources

1

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Variance of ML estimator: analytical method ★ Suppose we have n measurements of a random variable x with pdf f(x;𝜽)

★ So far we have computed several estimators of different 𝜃’s using ML. How to compute the variance (V[𝜃] or 𝜎2𝜃) of those estimators? ✦ Remember: 𝜎𝜃 quantifies the spread of the values of 𝜃 we would obtain after

performing many times the same experiment (always with n measurements)

★ In principle we could apply the definition of variance and compute it analytically. ✦ E.g. for the exponential distribution with mean 𝜏 estimated by 𝜏 = 1/n ∑i=1

nti

✦ In this case we already knew because that’s the estimator of the variance of the arithmetic mean

35

^ ^

^ ^

^

�(t) =

Z dt(t� E[t])2g(t)

=

Z dxn...

Z dx1 (t(x)� t(E[x]))2 f(x1)...f(xn) (21)

V [⌧̂ ] = E[⌧̂2]� (E[⌧̂ ])2

=

Z dtn...

Z dt1

1

n

nX

i=1

ti

!2 1

⌧ e�t1/⌧ ...

1

⌧ e�tn/⌧

� Z

dtn...

Z dt1

1

n

nX

i=1

ti

! 1

⌧ e�t1/⌧ ...

1

⌧ e�tn/⌧

!2

= ⌧2

n (22)

@ logL

@µ

���� µ=µ̂ML

= 0 ) � nX

i=1

(xi � µ̂ML)

�2 i

= 0 ) µ̂ML =

P n

i=1 xi/� 2 iP

n

i=1 1/� 2 i

(23)

logL(µ, {�i}i=1,...n) = n log 1p 2⇡

+ n

2 log

1

�2 i

� 1

2

nX

i=1

(xi � µ)2

�2 i

(24)

b[c�2 ML] = E[c�2

ML]� �2 = n� 1

n �2 � �2 =

�1

n �2 (25)

@ logL

@�2

���� �2=c�2

ML

= 0 ) nX

i=0

�1

2

1 c�2

ML

+ (xi � µ)2

2(c�2 ML)2

! = 0 ) c�2

ML = 1

n

nX

i=1

(xi � µ)2 (26)

@ logL

@µ

���� µ=µ̂ML

= 0 ) 1

2

nX

i=1

2(xi � µ̂ML)

�2 = 0 )

nX

i=1

xi = nµ̂ML ) µ̂ML = 1

n

nX

i=1

xi (27)

µ̂ML = 1

n

nX

i=1

xi (28)

f(x;µ,�) = 1p 2⇡�

e �(x�µ)2

�2 (29)

logL(µ,�) = nX

i=1

log f(xi;µ,�) = n log 1p 2⇡

+ n

2 log

1

�2 � 1

2

nX

i=1

(xi � µ)2

�2 (30)

b[�̂ML] = E[�̂ML]� � = � n

n� 1 � � =

1

n� 1 � (31)

�̂ML = 1

⌧̂ML =

nP m

i=1 ti (32)

3

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Comments to ML estimator variance analytical computation

★ Previous result V[𝜏] = 𝜏2/n we could have guessed since we already know that the variance of the arithmetic mean estimator x is 1/n times the variance of x (i.e. Var[x]=Var[x]/n), which for the exponential distribution is 𝜏2

★ 𝜏 is normally unknown, so in practice we will use an estimator of the

variance: V[𝜏] = 𝜏2/n = t2/n

✦ Remember that 𝜎𝜏 = 𝜏/√n, but that does not mean that the interval 𝜏 ± 𝜎𝜏 has 68.3% CL, which is true for a Gaussian distribution only

✦ When we report the result of an experiment as x = x ± 𝛥x, the “error” 𝛥x is normally interpreted as the 68.3% CL, and NOT as the standard deviation (they coincide for Gaussian estimators though, many according to the CLT)

★ The analytical method is (of course) always valid but normally too complicated or impracticable ✦ In the following slides we revise practical methods to estimate the variance

36

^

⋀ ^^ -

^̂ ^ ^^^

-

^ ^

-

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Variance of ML estimator: Monte Carlo method

★ We can estimate the variance of our estimator using Monte Carlo simulations ✦ Simulate a large number of experiments like the one we have performed

✦ As “true” value of the parameter we use our estimation

✦ For each simulated experiment we can compute the value of the estimator

✦ Using s2 (see before) we can estimate the variance (unbiased)

37

★ E.g. ✦ Expononential pdf, 𝜏=1, n=50 measurements lead to,

e.g. 𝜏 = t = 1.062

✦ Then we generate m=1000 MC “realizations”, each with 𝜏=1.062 (we don’t know the true 𝜏) and n=50

✦ For each realization we compute 𝜏i = ti

✦ Mean is 𝜏 = 1.059 (close to input)

✦ s = √s2 = √(1/(m-1)∑i=1 m(𝜏i-1.059)2) = 0.151 (close to

analytical estimator value 𝜎𝜏 = 𝜏/√n = 0.150

,-.... • to'

150

100

50

o o 0.5 1.5

Variance of ML estimators: the ReF bound 77

2

Fig. 6.3 A histogram of the ML es- timate T from 1000 Monte Carlo ex- periments with 50 observations per ex- periment. For the Monte Carlo 'true' parameter T, the result of Fig. 6.2 was used. The sample standard deviation is s = 0.151.

here holds are almost always met in practical situations (cf. [Ead71] Section 7.4.5). In the case of equality (i.e. minimum variance) the estimator is said to be efficient. It can be shown that if efficient estimators exist for a given problem, the maximum likelihood method will find them. Furthermore it can be shown that ML estimators are always efficient in the large sample limit, except when the extent of the sample space depends on the estimated parameter. In practice, one often assumes efficiency and zero bias. In cases of doubt one should check the results with a Monte Carlo study. The -general conditions for efficiency are discussed in, for example, [Ead71] Section 7.4.5, [Stu91] Chapter 18.

For the example of the exponential distribution with mean r one has from equation (6.5)

o2logL =!!..- =!!..- (1- 2f) (6.17) or2 r2 r n r2 r

i=1

and objor = 0 since b = 0 (see equation (6.7)). Thus the RCF bound for the variance (also called the minimum variance bound, or MVB) of T is

A 1 V[r]:2: E[-;2(1- 2:)]

1 n

(6.18)

where we have used equation (6.7) for E[ f]. Since r2 j n is also the variance obtained from the exact calculation (equation (6.15)) we see that equality holds and f = 2::7=1 ti is an efficient estimator for the parameter T.

For the case of more than one parameter, () = (()1, ... , Om), the correspond- ing formula for the inverse of the covariance matrix of their estimators Vij = cov [Oi , OJ] is (assuming efficiency and zero bias)

(6.19)

^ -

^ -

^

-̂

^^

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

The bootstrap method ★ Suppose we have our data sample x=(x1,…,xn) and we do not know

anything about their pdf f(x), not even a parametric expression f(x;𝜽), and Monte Carlo simulations are very costly

★ We compute some statistic of the data t(x), e.g. the mean, and we want to compute the statistical properties of the statistic t (e.g. its variance 𝜎t), but we do not know the pdf for the statistics g(t) either ✦ We cannot use the analytical method (

)

✦ We cannot use the Monte Carlo method

★ Bootstrap: ✦ We fabricate copies of the original distribution by drawing n measurements out of

x but allowing repetition, e.g. if n=10 (x7, x3, x2, x10, x1, x3, x5, x7, x9, x7)

✦ We compute t(x) for each of the copies and study the properties of t (e.g. variance, correlation, etc) from the resulting distribution like in the MC method

38

�(t) =

Z dt(t� E[t])2g(t)

=

Z dxn...

Z dx1 (t(x)� t(E[x]))2 f(x1)...f(xn) (1)

V [⌧̂ ] = E[⌧̂2]� (E[⌧ ])2

=

Z dtn...

Z dt1

1

n

nX

i=1

ti

!2 1

⌧ e�t1/⌧ ...

1

⌧ e�tn/⌧

� Z

dtn...

Z dt1

1

n

nX

i=1

ti

! 1

⌧ e�t1/⌧ ...

1

⌧ e�tn/⌧

!2

= ⌧2

n (2)

@ logL

@µ

���� µ=µ̂ML

= 0 ) � nX

i=1

(xi � µ̂ML)

�2 i

= 0 ) µ̂ML =

P n

i=1 xi/� 2 iP

n

i=1 1/� 2 i

(3)

logL(µ, {�i}i=1,...n) = n log 1p 2⇡

+ n

2 log

1

�2 i

� 1

2

nX

i=1

(xi � µ)2

�2 i

(4)

b[c�2 ML] = E[c�2

ML]� �2 = n� 1

n �2 � �2 =

�1

n �2 (5)

@ logL

@�2

���� �2=c�2

ML

= 0 ) nX

i=0

�1

2

1 c�2

ML

+ (xi � µ)2

2(c�2 ML)2

! = 0 ) c�2

ML = 1

n

nX

i=1

(xi � µ)2 (6)

@ logL

@µ

���� µ=µ̂ML

= 0 ) 1

2

nX

i=1

2(xi � µ̂ML)

�2 = 0 )

nX

i=1

xi = nµ̂ML ) µ̂ML = 1

n

nX

i=1

xi (7)

µ̂ML = 1

n

nX

i=1

xi (8)

f(x;µ,�) = 1p 2⇡�

e �(x�µ)2

�2 (9)

logL(µ,�) = nX

i=1

log f(xi;µ,�) = n log 1p 2⇡

+ n

2 log

1

�2 � 1

2

nX

i=1

(xi � µ)2

�2 (10)

b[�̂ML] = E[�̂ML]� � = � n

n� 1 � � =

1

n� 1 � (11)

�̂ML = 1

⌧̂ML =

nP m

i=1 ti (12)

1

�(t) =

Z dt(t� E[t])2g(t)

=

Z dxn...

Z dx1 (t(x)� t(E[x]))2 f(x1)...f(xn) (1)

V [⌧̂ ] = E[⌧̂2]� (E[⌧ ])2

=

Z dtn...

Z dt1

1

n

nX

i=1

ti

!2 1

⌧ e�t1/⌧ ...

1

⌧ e�tn/⌧

� Z

dtn...

Z dt1

1

n

nX

i=1

ti

! 1

⌧ e�t1/⌧ ...

1

⌧ e�tn/⌧

!2

= ⌧2

n (2)

@ logL

@µ

���� µ=µ̂ML

= 0 ) � nX

i=1

(xi � µ̂ML)

�2 i

= 0 ) µ̂ML =

P n

i=1 xi/� 2 iP

n

i=1 1/� 2 i

(3)

logL(µ, {�i}i=1,...n) = n log 1p 2⇡

+ n

2 log

1

�2 i

� 1

2

nX

i=1

(xi � µ)2

�2 i

(4)

b[c�2 ML] = E[c�2

ML]� �2 = n� 1

n �2 � �2 =

�1

n �2 (5)

@ logL

@�2

���� �2=c�2

ML

= 0 ) nX

i=0

�1

2

1 c�2

ML

+ (xi � µ)2

2(c�2 ML)2

! = 0 ) c�2

ML = 1

n

nX

i=1

(xi � µ)2 (6)

@ logL

@µ

���� µ=µ̂ML

= 0 ) 1

2

nX

i=1

2(xi � µ̂ML)

�2 = 0 )

nX

i=1

xi = nµ̂ML ) µ̂ML = 1

n

nX

i=1

xi (7)

µ̂ML = 1

n

nX

i=1

xi (8)

f(x;µ,�) = 1p 2⇡�

e �(x�µ)2

�2 (9)

logL(µ,�) = nX

i=1

log f(xi;µ,�) = n log 1p 2⇡

+ n

2 log

1

�2 � 1

2

nX

i=1

(xi � µ)2

�2 (10)

b[�̂ML] = E[�̂ML]� � = � n

n� 1 � � =

1

n� 1 � (11)

�̂ML = 1

⌧̂ML =

nP m

i=1 ti (12)

1

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

The jacknife method ★ Similar in spirit to bootstrap, we have a data sample x=(x1,…,xn)

and we want to study statistical properties of a given statistic t(x), e.g. its variance 𝜎2t, without knowing the pdf’s of x or t

★ Jacknife: ✦ Recompute t n times, each using the sample (x1,…xi-1,xi+1,…,xn) ∀i=1, …,n; i.e.: n samples of (n-1) measurements each, removing a different xi each time

✦ Thus we get n estimates of t and we can compute statistics of t, e.g. its variance

✦ For instance the variance of t can be estimated as

with and ti the value of t computed without xi

39

V [t] = n� 1

n

nX

i=1

(ti � t̄)2 (1)

�(t) =

Z dt(t� E[t])2g(t)

=

Z dxn...

Z dx1 (t(x)� t(E[x]))2 f(x1)...f(xn) (2)

V [⌧̂ ] = E[⌧̂2]� (E[⌧ ])2

=

Z dtn...

Z dt1

1

n

nX

i=1

ti

!2 1

⌧ e�t1/⌧ ...

1

⌧ e�tn/⌧

� Z

dtn...

Z dt1

1

n

nX

i=1

ti

! 1

⌧ e�t1/⌧ ...

1

⌧ e�tn/⌧

!2

= ⌧2

n (3)

@ logL

@µ

���� µ=µ̂ML

= 0 ) � nX

i=1

(xi � µ̂ML)

�2 i

= 0 ) µ̂ML =

P n

i=1 xi/� 2 iP

n

i=1 1/� 2 i

(4)

logL(µ, {�i}i=1,...n) = n log 1p 2⇡

+ n

2 log

1

�2 i

� 1

2

nX

i=1

(xi � µ)2

�2 i

(5)

b[c�2 ML] = E[c�2

ML]� �2 = n� 1

n �2 � �2 =

�1

n �2 (6)

@ logL

@�2

���� �2=c�2

ML

= 0 ) nX

i=0

�1

2

1 c�2

ML

+ (xi � µ)2

2(c�2 ML)2

! = 0 ) c�2

ML = 1

n

nX

i=1

(xi � µ)2 (7)

@ logL

@µ

���� µ=µ̂ML

= 0 ) 1

2

nX

i=1

2(xi � µ̂ML)

�2 = 0 )

nX

i=1

xi = nµ̂ML ) µ̂ML = 1

n

nX

i=1

xi (8)

µ̂ML = 1

n

nX

i=1

xi (9)

f(x;µ,�) = 1p 2⇡�

e �(x�µ)2

�2 (10)

logL(µ,�) = nX

i=1

log f(xi;µ,�) = n log 1p 2⇡

+ n

2 log

1

�2 � 1

2

nX

i=1

(xi � µ)2

�2 (11)

b[�̂ML] = E[�̂ML]� � = � n

n� 1 � � =

1

n� 1 � (12)

1

t̄ = 1

n

nX

i=1

ti (1)

V [t] = n� 1

n

nX

i=1

(ti � t̄)2 (2)

�(t) =

Z dt(t� E[t])2g(t)

=

Z dxn...

Z dx1 (t(x)� t(E[x]))2 f(x1)...f(xn) (3)

V [⌧̂ ] = E[⌧̂2]� (E[⌧ ])2

=

Z dtn...

Z dt1

1

n

nX

i=1

ti

!2 1

⌧ e�t1/⌧ ...

1

⌧ e�tn/⌧

� Z

dtn...

Z dt1

1

n

nX

i=1

ti

! 1

⌧ e�t1/⌧ ...

1

⌧ e�tn/⌧

!2

= ⌧2

n (4)

@ logL

@µ

���� µ=µ̂ML

= 0 ) � nX

i=1

(xi � µ̂ML)

�2 i

= 0 ) µ̂ML =

P n

i=1 xi/� 2 iP

n

i=1 1/� 2 i

(5)

logL(µ, {�i}i=1,...n) = n log 1p 2⇡

+ n

2 log

1

�2 i

� 1

2

nX

i=1

(xi � µ)2

�2 i

(6)

b[c�2 ML] = E[c�2

ML]� �2 = n� 1

n �2 � �2 =

�1

n �2 (7)

@ logL

@�2

���� �2=c�2

ML

= 0 ) nX

i=0

�1

2

1 c�2

ML

+ (xi � µ)2

2(c�2 ML)2

! = 0 ) c�2

ML = 1

n

nX

i=1

(xi � µ)2 (8)

@ logL

@µ

���� µ=µ̂ML

= 0 ) 1

2

nX

i=1

2(xi � µ̂ML)

�2 = 0 )

nX

i=1

xi = nµ̂ML ) µ̂ML = 1

n

nX

i=1

xi (9)

µ̂ML = 1

n

nX

i=1

xi (10)

f(x;µ,�) = 1p 2⇡�

e �(x�µ)2

�2 (11)

1

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

The Rao-Cramér-Frechet bound ★ Sometimes it is difficult to compute the variances analytically and

Monte Carlo simulations can be computationally expensive

★ The Rao-Cramér-Frechet (RCF) inequality provides a lower bound to an estimator’s variance ✦ “The variance of an unbiased estimator is at least as high as the inverse

of the Fisher information”, i.e.:

★ If an estimator has the lowest variance is said to be efficient ★ It can be shown that if efficient estimators exists for a given

problem, the ML method will find them

★ It can also be shown that ML estimators are always efficient for n→∞

40

V [✓̂] � E

 �@2 logL

@✓2

��1

(1)

t̄ = 1

n

nX

i=1

ti (2)

V [t] = n� 1

n

nX

i=1

(ti � t̄)2 (3)

�(t) =

Z dt(t� E[t])2g(t)

=

Z dxn...

Z dx1 (t(x)� t(E[x]))2 f(x1)...f(xn) (4)

V [⌧̂ ] = E[⌧̂2]� (E[⌧ ])2

=

Z dtn...

Z dt1

1

n

nX

i=1

ti

!2 1

⌧ e�t1/⌧ ...

1

⌧ e�tn/⌧

� Z

dtn...

Z dt1

1

n

nX

i=1

ti

! 1

⌧ e�t1/⌧ ...

1

⌧ e�tn/⌧

!2

= ⌧2

n (5)

@ logL

@µ

���� µ=µ̂ML

= 0 ) � nX

i=1

(xi � µ̂ML)

�2 i

= 0 ) µ̂ML =

P n

i=1 xi/� 2 iP

n

i=1 1/� 2 i

(6)

logL(µ, {�i}i=1,...n) = n log 1p 2⇡

+ n

2 log

1

�2 i

� 1

2

nX

i=1

(xi � µ)2

�2 i

(7)

b[c�2 ML] = E[c�2

ML]� �2 = n� 1

n �2 � �2 =

�1

n �2 (8)

@ logL

@�2

���� �2=c�2

ML

= 0 ) nX

i=0

�1

2

1 c�2

ML

+ (xi � µ)2

2(c�2 ML)2

! = 0 ) c�2

ML = 1

n

nX

i=1

(xi � µ)2 (9)

@ logL

@µ

���� µ=µ̂ML

= 0 ) 1

2

nX

i=1

2(xi � µ̂ML)

�2 = 0 )

nX

i=1

xi = nµ̂ML ) µ̂ML = 1

n

nX

i=1

xi (10)

µ̂ML = 1

n

nX

i=1

xi (11)

1

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

RCF bound for the exponential distribution

★ E.g., for the exponential distribution we have

so

and therefore the equality holds, which means that 𝜏 = 1/n ∑i=1 nti is an efficient

estimator of the parameter 𝜏, as expected because it’s the ML estimator

41

L(⌧ ; {ti}i=1,...n) = nY

i=1

✓ 1

⌧ e�ti/⌧

◆ )

logL = �n log ⌧ � 1

⌧

nX

i=1

ti )

@2 logL

@⌧2 =

n

⌧2

1� 2

⌧

1

n

nX

i=1

ti

! =

n

⌧2

✓ 1� 2⌧̂

⌧

◆ (1)

V [✓̂] � E

 �@2 logL

@✓2

��1

(2)

t̄ = 1

n

nX

i=1

ti (3)

V [t] = n� 1

n

nX

i=1

(ti � t̄)2 (4)

�(t) =

Z dt(t� E[t])2g(t)

=

Z dxn...

Z dx1 (t(x)� t(E[x]))2 f(x1)...f(xn) (5)

V [⌧̂ ] = E[⌧̂2]� (E[⌧ ])2

=

Z dtn...

Z dt1

1

n

nX

i=1

ti

!2 1

⌧ e�t1/⌧ ...

1

⌧ e�tn/⌧

� Z

dtn...

Z dt1

1

n

nX

i=1

ti

! 1

⌧ e�t1/⌧ ...

1

⌧ e�tn/⌧

!2

= ⌧2

n (6)

@ logL

@µ

���� µ=µ̂ML

= 0 ) � nX

i=1

(xi � µ̂ML)

�2 i

= 0 ) µ̂ML =

P n

i=1 xi/� 2 iP

n

i=1 1/� 2 i

(7)

logL(µ, {�i}i=1,...n) = n log 1p 2⇡

+ n

2 log

1

�2 i

� 1

2

nX

i=1

(xi � µ)2

�2 i

(8)

b[c�2 ML] = E[c�2

ML]� �2 = n� 1

n �2 � �2 =

�1

n �2 (9)

@ logL

@�2

���� �2=c�2

ML

= 0 ) nX

i=0

�1

2

1 c�2

ML

+ (xi � µ)2

2(c�2 ML)2

! = 0 ) c�2

ML = 1

n

nX

i=1

(xi � µ)2 (10)

1

V [⌧̂ ] � 1

E ⇥ � n

⌧2 (1� 2⌧̂

⌧ ) ⇤ = 1

� n

⌧2 (1� 2E[⌧̂ ]

⌧ ) =

⌧2

n (1)

L(⌧ ; {ti}i=1,...n) = nY

i=1

✓ 1

⌧ e�ti/⌧

◆ )

logL = �n log ⌧ � 1

⌧

nX

i=1

ti )

@2 logL

@⌧2 =

n

⌧2

1� 2

⌧

1

n

nX

i=1

ti

! =

n

⌧2

✓ 1� 2⌧̂

⌧

◆ (2)

V [✓̂] � E

 �@2 logL

@✓2

��1

(3)

t̄ = 1

n

nX

i=1

ti (4)

V [t] = n� 1

n

nX

i=1

(ti � t̄)2 (5)

�(t) =

Z dt(t� E[t])2g(t)

=

Z dxn...

Z dx1 (t(x)� t(E[x]))2 f(x1)...f(xn) (6)

V [⌧̂ ] = E[⌧̂2]� (E[⌧ ])2

=

Z dtn...

Z dt1

1

n

nX

i=1

ti

!2 1

⌧ e�t1/⌧ ...

1

⌧ e�tn/⌧

� Z

dtn...

Z dt1

1

n

nX

i=1

ti

! 1

⌧ e�t1/⌧ ...

1

⌧ e�tn/⌧

!2

= ⌧2

n (7)

@ logL

@µ

���� µ=µ̂ML

= 0 ) � nX

i=1

(xi � µ̂ML)

�2 i

= 0 ) µ̂ML =

P n

i=1 xi/� 2 iP

n

i=1 1/� 2 i

(8)

logL(µ, {�i}i=1,...n) = n log 1p 2⇡

+ n

2 log

1

�2 i

� 1

2

nX

i=1

(xi � µ)2

�2 i

(9)

1

^

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

RCF with multiple parameters ★ For multiple parameters 𝜽=(𝜃1,…,𝜃m), the RCF formula for the inverse of the

covariance matrix of efficient estimators is (assuming zero bias):

called the Fisher information matrix ★ Computing the expectation values is not trivial. Typically for large n, one

replaces this formula with its value at 𝜽=𝜽:

which, for one parameter 𝜃 reduces to:

(These formulas are normally in numerical methods for estimating the covariance matrix of the ML estimators)

42

(V �1)ij = E

 �@2 logL

@✓i@✓j

� (1)

V [⌧̂ ] � 1

E ⇥ � n

⌧2 (1� 2⌧̂

⌧ ) ⇤ = 1

� n

⌧2 (1� 2E[⌧̂ ]

⌧ ) =

⌧2

n (2)

L(⌧ ; {ti}i=1,...n) = nY

i=1

✓ 1

⌧ e�ti/⌧

◆ )

logL = �n log ⌧ � 1

⌧

nX

i=1

ti )

@2 logL

@⌧2 =

n

⌧2

1� 2

⌧

1

n

nX

i=1

ti

! =

n

⌧2

✓ 1� 2⌧̂

⌧

◆ (3)

V [✓̂] � E

 �@2 logL

@✓2

��1

(4)

t̄ = 1

n

nX

i=1

ti (5)

V [t] = n� 1

n

nX

i=1

(ti � t̄)2 (6)

�(t) =

Z dt(t� E[t])2g(t)

=

Z dxn...

Z dx1 (t(x)� t(E[x]))2 f(x1)...f(xn) (7)

V [⌧̂ ] = E[⌧̂2]� (E[⌧ ])2

=

Z dtn...

Z dt1

1

n

nX

i=1

ti

!2 1

⌧ e�t1/⌧ ...

1

⌧ e�tn/⌧

� Z

dtn...

Z dt1

1

n

nX

i=1

ti

! 1

⌧ e�t1/⌧ ...

1

⌧ e�tn/⌧

!2

= ⌧2

n (8)

@ logL

@µ

���� µ=µ̂ML

= 0 ) � nX

i=1

(xi � µ̂ML)

�2 i

= 0 ) µ̂ML =

P n

i=1 xi/� 2 iP

n

i=1 1/� 2 i

(9)

1

^

(dV �1)ij = �@2 logL

@✓i@✓j

���� ✓=✓̂

(1)

(V �1)ij = E

 �@2 logL

@✓i@✓j

� (2)

V [⌧̂ ] � 1

E ⇥ � n

⌧2 (1� 2⌧̂

⌧ ) ⇤ = 1

� n

⌧2 (1� 2E[⌧̂ ]

⌧ ) =

⌧2

n (3)

L(⌧ ; {ti}i=1,...n) = nY

i=1

✓ 1

⌧ e�ti/⌧

◆ )

logL = �n log ⌧ � 1

⌧

nX

i=1

ti )

@2 logL

@⌧2 =

n

⌧2

1� 2

⌧

1

n

nX

i=1

ti

! =

n

⌧2

✓ 1� 2⌧̂

⌧

◆ (4)

V [✓̂] � E

 �@2 logL

@✓2

��1

(5)

t̄ = 1

n

nX

i=1

ti (6)

V [t] = n� 1

n

nX

i=1

(ti � t̄)2 (7)

�(t) =

Z dt(t� E[t])2g(t)

=

Z dxn...

Z dx1 (t(x)� t(E[x]))2 f(x1)...f(xn) (8)

V [⌧̂ ] = E[⌧̂2]� (E[⌧ ])2

=

Z dtn...

Z dt1

1

n

nX

i=1

ti

!2 1

⌧ e�t1/⌧ ...

1

⌧ e�tn/⌧

� Z

dtn...

Z dt1

1

n

nX

i=1

ti

! 1

⌧ e�t1/⌧ ...

1

⌧ e�tn/⌧

!2

= ⌧2

n (9)

1

logL(✓) = logLmax � (✓ � ✓̂)2

2c�2 ✓̂

(1)

logL(✓) = logL(✓̂) +

 @ logL

@✓

�

✓=✓̂

(✓ � ✓̂) + 1

2!

 @2 logL

@✓2

�

✓=✓̂

(✓ � ✓̂)2 + ... (2)

c�2 ✓̂ =

1

@2 logL @✓2

!����� ✓=✓̂

(3)

(V �1)ij = E

 �@2 logL

@✓i@✓j

� (4)

V [⌧̂ ] � 1

E ⇥ � n

⌧2 (1� 2⌧̂

⌧ ) ⇤ = 1

� n

⌧2 (1� 2E[⌧̂ ]

⌧ ) =

⌧2

n (5)

L(⌧ ; {ti}i=1,...n) = nY

i=1

✓ 1

⌧ e�ti/⌧

◆ )

logL = �n log ⌧ � 1

⌧

nX

i=1

ti )

@2 logL

@⌧2 =

n

⌧2

1� 2

⌧

1

n

nX

i=1

ti

! =

n

⌧2

✓ 1� 2⌧̂

⌧

◆ (6)

V [✓̂] � E

 �@2 logL

@✓2

��1

(7)

t̄ = 1

n

nX

i=1

ti (8)

V [t] = n� 1

n

nX

i=1

(ti � t̄)2 (9)

�(t) =

Z dt(t� E[t])2g(t)

=

Z dxn...

Z dx1 (t(x)� t(E[x]))2 f(x1)...f(xn) (10)

1

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Variance of ML estimators: the logLmax-1/2 rule

★ Let us consider the case of one parameter 𝜃 and Taylor-expand logL(𝜃) around its ML estimator 𝜃:

✦ By definition logL(𝜃)= logLmax(𝜃)

✦ The first derivative (second term) is 0 at the maximum

★ Ignoring higher order terms (0 if L is Gaussian) and using the expression for 𝜎2𝜃 found in previous slide we have

which, evaluated at 𝜃±𝜎2𝜃 yields:

That is, a change in the parameter 𝜃 by 1 standard deviation from its ML estimate leads to a decrease in the log-likelihood of 1/2 from the maximum

43

logL(✓) = logL(✓̂) +

 @ logL

@✓

�

✓=✓̂

(✓ � ✓̂) + 1

2!

 @2 logL

@✓2

�

✓=✓̂

(✓ � ✓̂)2 + ... (1)

c�2 ✓̂ =

1

@2 logL @✓2

!����� ✓=✓̂

(2)

(V �1)ij = E

 �@2 logL

@✓i@✓j

� (3)

V [⌧̂ ] � 1

E ⇥ � n

⌧2 (1� 2⌧̂

⌧ ) ⇤ = 1

� n

⌧2 (1� 2E[⌧̂ ]

⌧ ) =

⌧2

n (4)

L(⌧ ; {ti}i=1,...n) = nY

i=1

✓ 1

⌧ e�ti/⌧

◆ )

logL = �n log ⌧ � 1

⌧

nX

i=1

ti )

@2 logL

@⌧2 =

n

⌧2

1� 2

⌧

1

n

nX

i=1

ti

! =

n

⌧2

✓ 1� 2⌧̂

⌧

◆ (5)

V [✓̂] � E

 �@2 logL

@✓2

��1

(6)

t̄ = 1

n

nX

i=1

ti (7)

V [t] = n� 1

n

nX

i=1

(ti � t̄)2 (8)

�(t) =

Z dt(t� E[t])2g(t)

=

Z dxn...

Z dx1 (t(x)� t(E[x]))2 f(x1)...f(xn) (9)

V [⌧̂ ] = E[⌧̂2]� (E[⌧ ])2

=

Z dtn...

Z dt1

1

n

nX

i=1

ti

!2 1

⌧ e�t1/⌧ ...

1

⌧ e�tn/⌧

� Z

dtn...

Z dt1

1

n

nX

i=1

ti

! 1

⌧ e�t1/⌧ ...

1

⌧ e�tn/⌧

!2

= ⌧2

n (10)

1

^

^

logL(✓) = logLmax � (✓ � ✓̂)2

2c�2 ✓̂

(1)

logL(✓) = logL(✓̂) +

 @ logL

@✓

�

✓=✓̂

(✓ � ✓̂) + 1

2!

 @2 logL

@✓2

�

✓=✓̂

(✓ � ✓̂)2 + ... (2)

c�2 ✓̂ =

1

@2 logL @✓2

!����� ✓=✓̂

(3)

(V �1)ij = E

 �@2 logL

@✓i@✓j

� (4)

V [⌧̂ ] � 1

E ⇥ � n

⌧2 (1� 2⌧̂

⌧ ) ⇤ = 1

� n

⌧2 (1� 2E[⌧̂ ]

⌧ ) =

⌧2

n (5)

L(⌧ ; {ti}i=1,...n) = nY

i=1

✓ 1

⌧ e�ti/⌧

◆ )

logL = �n log ⌧ � 1

⌧

nX

i=1

ti )

@2 logL

@⌧2 =

n

⌧2

1� 2

⌧

1

n

nX

i=1

ti

! =

n

⌧2

✓ 1� 2⌧̂

⌧

◆ (6)

V [✓̂] � E

 �@2 logL

@✓2

��1

(7)

t̄ = 1

n

nX

i=1

ti (8)

V [t] = n� 1

n

nX

i=1

(ti � t̄)2 (9)

�(t) =

Z dt(t� E[t])2g(t)

=

Z dxn...

Z dx1 (t(x)� t(E[x]))2 f(x1)...f(xn) (10)

1

^

^ ^

^

^

⌫i(✓) = ⌫tot

Z x max i

xmin i

f(x;✓)dx (1)

logL(⌫tot,✓) = �⌫tot + NX

i=1

ni log ⌫i(⌫tot,✓) (2)

f(n;⌫) = ⌫ntot tot e�⌫tot

ntot!

ntot!

n1!...nN !

✓ ⌫1 ⌫tot

◆n1

...

✓ ⌫N ⌫tot

◆nN

= NY

i=1

⌫ni i

ni! e�⌫ (3)

logL(✓) = NX

i=1

ni log ⌫i(✓) (4)

f(n;⌫) = ntot!

n1!...nN !

✓ ⌫1 ntot

◆n1

...

✓ ⌫N ntot

◆nN

(5)

⌫i(✓) = ntot

Z x max i

xmin i

f(x;✓)dx (6)

logL(✓̂ ±c�2 ✓̂ ) = logLmax �

1

2 (7)

logL(✓) = logLmax � (✓ � ✓̂)2

2c�2 ✓̂

(8)

logL(✓) = logL(✓̂) +

 @ logL

@✓

�

✓=✓̂

(✓ � ✓̂) + 1

2!

 @2 logL

@✓2

�

✓=✓̂

(✓ � ✓̂)2 + ... (9)

c�2 ✓̂ =

1

@2 logL @✓2

!����� ✓=✓̂

(10)

(V �1)ij = E

 �@2 logL

@✓i@✓j

� (11)

V [⌧̂ ] � 1

E ⇥ � n

⌧2 (1� 2⌧̂

⌧ ) ⇤ = 1

� n

⌧2 (1� 2E[⌧̂ ]

⌧ ) =

⌧2

n (12)

1

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Variance of ML estimators: the logLmax-1/2 rule

★ Similar use as for the LS 𝜒2+1 rule: ✦ Provides estimator of standard deviation of the

estimator (i.e. 𝜎2𝜃)

✦ If the pdf of 𝜃 is Gaussian: provides contours with tabulated CL

★ Remember the assumptions we have made to arrive here: ✦ The estimator 𝜃 is unbiased

✦ The estimator 𝜃 is efficient

✦ The expectation value of logL derivatives are approximated by their values at 𝜃

✦ The likelihood function is Gaussian (or equivalently the logL is a parabola)

★ All assumptions are fulfilled for n→∞

44

Variance of ML estimators: graphical method 79

By definition of B we know that log L (B) = log Lmax and that the second term in the expansion is zero. Using equation (6.22) and ignoring higher order terms gIves

(8 - 0)2 log L(8) = log Lmax - __ ,

20-2 e (6.24)

or

A 1 log L(8 ± O"e) = log Lmax - "2. (6.25)

That is, a change in the parameter 8 of one standard deviation from its ML estimate leads to a decrease in the log-likelihood of 1/2 from its maximum value.

It can be shown that the log-likelihood function becomes a parabola (i.e. the likelihood function becomes a Gaussian curve) in the large sample limit. Even if log L is not parabolic, one can nevertheless adopt equation (6.25) as the definition of the statistical error. The interpretation of such errors is discussed further in Chapter 9.

As an example of the graphical method for determining the variance of an es- timator, consider again the examples of Sections 6.2 and 6.5 with the exponential distribution. Figure 6.4 shows the log-likelihood function log L( r) as a function of the parameter r for a Monte Carlo experiment consisting of 50 measurements. The standard deviation of f is estimated by changing r until log L( r) decreases by 1/2, giving Llf_ = 0.137, Llf+ = 0.165. In this case logL(r) is reasonably close to a parabola and one can approximate 0" f Ll f _ Ll f + 0.15. This leads to approximately the same answer as from the exact standard deviation r /...;n evaluated with r = f. In Chapter 9 the interval [f - Ll f _ , f + Ll f +] will be reinterpreted as an approximation for the 68.3% central confidence interval (cf. Section 9.6).

-53

-53.5

-54 0.8 1.2 1.4 1.6

Fig. 6.4 The log-likelihood function logL(T). In the large sample limit, the widths of the intervals [i-Lli-,il and [i,f + Llf+l correspond to one stan- dard deviation at.

82 The method of maximum likelihood

The sample means, standard deviations, covariance and correlation coefficient (see Section 5.2) from the Monte Carlo experiments are:

0.499 0.051 0.0024

(3

r

0.498 0.111 0.42.

(6.29)

Note that & and /3 are in good agreement with the 'true' values put into the Monte Carlo (0' = 0.5 and (3 = 0.5) and the sample (co)variances are close to the values estimated numerically from the ReF bound.

The fact that a and /3 are correlated is seen from the fact that the band of points in the scatter plot is tilted. That is, if one required a > 0', this would lead to an enhanced probability to also find /3 > (3. In other words, the conditional p.d.f. for a given /3 > (3 is centered at a higher mean value and has a smaller variance than the marginal p.d.f. for a.

Figure 6.7 shows the positions of the ML estimates in the parameter space along with a contour corresponding to log L = log Lmax - 1/2.

0.7

0.6 ····t··································· ----·---------r----

0.5 true value Fig. 6.7 The contour of constant likelihood logL = logLmax - 1/2 shown with the true values for the par- ameters (a,.6) and the ML estimates (a,t1). In the large sample limit the tangents to the curve correspond to a ± u& and t1 ± uiJ.

0.4

i MLf", •• ,"

0.3

·· .. ··········l· .. ·········· .. · .. ······ .. ··· .. · ..

0.3 0.4 0.5 0.6 0.7

<X

In the large sample limit, the log-likelihood function takes on the form

log L( 0', (3) = log Lmax

[C';'&)' + -2P(";'&) ,(6.30)

where p = cov[a,/3]/(O"&O"{§) is the correlation coefficient for a and /3. The contour of log L (0', (3) = log Lmax - 1/2 is thus given by

^ 𝛽±σ𝛽

α±σα

^ ^

^

^ ^

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Summary of ML estimators ★ ML estimators are consistent (𝜃→𝜃0 for n→∞)

★ ML estimators are invariant under transformation: if a(𝜃) is monotonic function, then a = a(𝜃)

★ Often the ML estimator is biased: E[𝜃]≠𝜃0

★ ML estimators are asymptotically unbiased: E[𝜃]=𝜃0 for n→∞

★ ML estimators are asymptotically efficient: V[𝜃] is minimal (according to RCF inequality) for n→∞

★ ML estimators are asymptotically normal (Gaussian distributed): g(𝜃) is Gaussian about 𝜃 for n→∞ (from CLT)

45

^

^ ^

^

^

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

ML with binned data ★ Sometimes, the number of

measurements n is so high, that ML is impractical, since one must sum log f(xi;𝜃) for each value xi.

★ In such cases, instead of recording individual xi values, we make an histogram with the number of times a measurement lies within each of its N bins, i.e. n = (n1,…,nN).

★ The expectation values 𝝂 = (𝜈1,…,𝜈N) of the n are given by

where ximin and ximax are the bin limits and ntot = ∑i=1

Nni

46

⌫i(✓) = ntot

Z x max i

xmin i

f(x;✓)dx (1)

logL(✓ ±c�2 ✓̂ ) = logLmax �

1

2 (2)

logL(✓) = logLmax � (✓ � ✓̂)2

2c�2 ✓̂

(3)

logL(✓) = logL(✓̂) +

 @ logL

@✓

�

✓=✓̂

(✓ � ✓̂) + 1

2!

 @2 logL

@✓2

�

✓=✓̂

(✓ � ✓̂)2 + ... (4)

c�2 ✓̂ =

1

@2 logL @✓2

!����� ✓=✓̂

(5)

(V �1)ij = E

 �@2 logL

@✓i@✓j

� (6)

V [⌧̂ ] � 1

E ⇥ � n

⌧2 (1� 2⌧̂

⌧ ) ⇤ = 1

� n

⌧2 (1� 2E[⌧̂ ]

⌧ ) =

⌧2

n (7)

L(⌧ ; {ti}i=1,...n) = nY

i=1

✓ 1

⌧ e�ti/⌧

◆ )

logL = �n log ⌧ � 1

⌧

nX

i=1

ti )

@2 logL

@⌧2 =

n

⌧2

1� 2

⌧

1

n

nX

i=1

ti

! =

n

⌧2

✓ 1� 2⌧̂

⌧

◆ (8)

V [✓̂] � E

 �@2 logL

@✓2

��1

(9)

t̄ = 1

n

nX

i=1

ti (10)

V [t] = n� 1

n

nX

i=1

(ti � t̄)2 (11)

1

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

ML with binned data ★ We can regard the histogram as a single measurement of an N-dimensional random vector for

which the joint pdf is the multinomial distribution, i.e.

where pi has been expressed as pi =𝜈i/ntot

★ Taking the log and removing terms that do not depend on the parameters 𝝂 we get

which we maximize by any means at our hands (generally numerically)

★ We may want to consider ntot as a random variable from a Poisson distribution with mean 𝜈tot, in which case

i.e. equivalent to treating the number of entries in each bin as independent Poisson variable with mean 𝝂i. The logL is:

47

f(n;⌫) = ntot!

n1!...nN !

✓ ⌫1 ntot

◆n1

...

✓ ⌫N ntot

◆nN

(1)

⌫i(✓) = ntot

Z x max i

xmin i

f(x;✓)dx (2)

logL(✓ ±c�2 ✓̂ ) = logLmax �

1

2 (3)

logL(✓) = logLmax � (✓ � ✓̂)2

2c�2 ✓̂

(4)

logL(✓) = logL(✓̂) +

 @ logL

@✓

�

✓=✓̂

(✓ � ✓̂) + 1

2!

 @2 logL

@✓2

�

✓=✓̂

(✓ � ✓̂)2 + ... (5)

c�2 ✓̂ =

1

@2 logL @✓2

!����� ✓=✓̂

(6)

(V �1)ij = E

 �@2 logL

@✓i@✓j

� (7)

V [⌧̂ ] � 1

E ⇥ � n

⌧2 (1� 2⌧̂

⌧ ) ⇤ = 1

� n

⌧2 (1� 2E[⌧̂ ]

⌧ ) =

⌧2

n (8)

L(⌧ ; {ti}i=1,...n) = nY

i=1

✓ 1

⌧ e�ti/⌧

◆ )

logL = �n log ⌧ � 1

⌧

nX

i=1

ti )

@2 logL

@⌧2 =

n

⌧2

1� 2

⌧

1

n

nX

i=1

ti

! =

n

⌧2

✓ 1� 2⌧̂

⌧

◆ (9)

V [✓̂] � E

 �@2 logL

@✓2

��1

(10)

t̄ = 1

n

nX

i=1

ti (11)

1

logL(⌫tot,✓) = �⌫tot + NX

i=1

ni log ⌫i(⌫tot,✓) (1)

f(n;⌫) = ⌫ntot tot e�⌫tot

ntot!

ntot!

n1!...nN !

✓ ⌫1 ⌫tot

◆n1

...

✓ ⌫N ⌫tot

◆nN

= NY

i=1

⌫ni i

ni! e�⌫ (2)

logL(✓) = NX

i=1

ni log ⌫i(✓) (3)

f(n;⌫) = ntot!

n1!...nN !

✓ ⌫1 ntot

◆n1

...

✓ ⌫N ntot

◆nN

(4)

⌫i(✓) = ntot

Z x max i

xmin i

f(x;✓)dx (5)

logL(✓ ±c�2 ✓̂ ) = logLmax �

1

2 (6)

logL(✓) = logLmax � (✓ � ✓̂)2

2c�2 ✓̂

(7)

logL(✓) = logL(✓̂) +

 @ logL

@✓

�

✓=✓̂

(✓ � ✓̂) + 1

2!

 @2 logL

@✓2

�

✓=✓̂

(✓ � ✓̂)2 + ... (8)

c�2 ✓̂ =

1

@2 logL @✓2

!����� ✓=✓̂

(9)

(V �1)ij = E

 �@2 logL

@✓i@✓j

� (10)

V [⌧̂ ] � 1

E ⇥ � n

⌧2 (1� 2⌧̂

⌧ ) ⇤ = 1

� n

⌧2 (1� 2E[⌧̂ ]

⌧ ) =

⌧2

n (11)

1

logL(⌫tot,✓) = �⌫tot + NX

i=1

ni log ⌫i(⌫tot,✓) (1)

f(n;⌫) = ⌫ntot tot e�⌫tot

ntot!

ntot!

n1!...nN !

✓ ⌫1 ⌫tot

◆n1

...

✓ ⌫N ⌫tot

◆nN

= NY

i=1

⌫ni i

ni! e�⌫ (2)

logL(✓) = NX

i=1

ni log ⌫i(✓) (3)

f(n;⌫) = ntot!

n1!...nN !

✓ ⌫1 ntot

◆n1

...

✓ ⌫N ntot

◆nN

(4)

⌫i(✓) = ntot

Z x max i

xmin i

f(x;✓)dx (5)

logL(✓ ±c�2 ✓̂ ) = logLmax �

1

2 (6)

logL(✓) = logLmax � (✓ � ✓̂)2

2c�2 ✓̂

(7)

logL(✓) = logL(✓̂) +

 @ logL

@✓

�

✓=✓̂

(✓ � ✓̂) + 1

2!

 @2 logL

@✓2

�

✓=✓̂

(✓ � ✓̂)2 + ... (8)

c�2 ✓̂ =

1

@2 logL @✓2

!����� ✓=✓̂

(9)

(V �1)ij = E

 �@2 logL

@✓i@✓j

� (10)

V [⌧̂ ] � 1

E ⇥ � n

⌧2 (1� 2⌧̂

⌧ ) ⇤ = 1

� n

⌧2 (1� 2E[⌧̂ ]

⌧ ) =

⌧2

n (11)

1

⌫i(✓) = ⌫tot

Z x max i

xmin i

f(x;✓)dx (1)

logL(⌫tot,✓) = �⌫tot + NX

i=1

ni log ⌫i(⌫tot,✓) (2)

f(n;⌫) = ⌫ntot tot e�⌫tot

ntot!

ntot!

n1!...nN !

✓ ⌫1 ⌫tot

◆n1

...

✓ ⌫N ⌫tot

◆nN

= NY

i=1

⌫ni i

ni! e�⌫ (3)

logL(✓) = NX

i=1

ni log ⌫i(✓) (4)

f(n;⌫) = ntot!

n1!...nN !

✓ ⌫1 ntot

◆n1

...

✓ ⌫N ntot

◆nN

(5)

⌫i(✓) = ntot

Z x max i

xmin i

f(x;✓)dx (6)

logL(✓ ±c�2 ✓̂ ) = logLmax �

1

2 (7)

logL(✓) = logLmax � (✓ � ✓̂)2

2c�2 ✓̂

(8)

logL(✓) = logL(✓̂) +

 @ logL

@✓

�

✓=✓̂

(✓ � ✓̂) + 1

2!

 @2 logL

@✓2

�

✓=✓̂

(✓ � ✓̂)2 + ... (9)

c�2 ✓̂ =

1

@2 logL @✓2

!����� ✓=✓̂

(10)

(V �1)ij = E

 �@2 logL

@✓i@✓j

� (11)

V [⌧̂ ] � 1

E ⇥ � n

⌧2 (1� 2⌧̂

⌧ ) ⇤ = 1

� n

⌧2 (1� 2E[⌧̂ ]

⌧ ) =

⌧2

n (12)

1

(with )

[a, b] = [✓̂obs � � ✓̂ , ✓̂obs + �

✓̂ ] (1)

a = ✓̂obs � � ✓̂ ��1(1� ↵)

b = ✓̂obs + � ✓̂ ��1(1� �) (2)

↵ =

Z 1

✓̂obs

1p 2⇡�

✓̂

e � 1

2 ( ✓̂�a � ✓̂ )2

d✓̂ =

( x = ✓̂�a

�✓̂

dx = 1 �✓̂ d✓̂

) =

Z 1

✓̂obs�a � ✓̂

1p 2⇡

e�x 2 /2 dx = 1� �(

✓̂obs � a

� ✓̂

)

� =

Z ✓̂obs

�1

1p 2⇡�

✓̂

e � 1

2 ( ✓̂�b � ✓̂ )2

d✓̂ =

( x = ✓̂�b

�✓̂

dx = 1 �✓̂ d✓̂

) =

Z ✓̂obs�b � ✓̂

�1

1p 2⇡

e�x 2 /2 dx = �(

✓̂obs � b

� ✓̂

)

✓ = ✓̂obs +(a� ✓̂obs) �(✓̂obs � b)

(3)

✓ = ✓̂obs ±�✓̂ (4)

⌫i(✓) = ⌫tot

Z x max i

xmin i

f(x;✓)dx (5)

logL(⌫tot,✓) = �⌫tot + NX

i=1

ni log ⌫i(⌫tot,✓) (6)

f(n;⌫) = ⌫ntot tot e�⌫tot

ntot!

ntot!

n1!...nN !

✓ ⌫1 ⌫tot

◆n1

...

✓ ⌫N ⌫tot

◆nN

= NY

i=1

⌫ni i

ni! e�⌫i (7)

logL(✓) = NX

i=1

ni log ⌫i(✓) (8)

f(n;⌫) = ntot!

n1!...nN !

✓ ⌫1 ntot

◆n1

...

✓ ⌫N ntot

◆nN

(9)

⌫i(✓) = ntot

Z x max i

xmin i

f(x;✓)dx (10)

1

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Exercises 7. Consider the height data from Exercise 1 (first column in

ParEst_input_1.txt), asume a Gaussian pdf and estimate, maximizing numerically the likelihood:

a. the population mean and standard deviation

b. the standard deviation of both estimators using i. Monte Carlo simulations

ii. The bootstrap method

iii. The jacknife method

iv. The logLmin + 1/2 rule

c. bin the data in a histogram of 7 bins between 140 and 210 cm and estimate the population mean and standard deviation

d. compare results with those from Exercise 1

48

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Confidence intervals

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Standard deviation as statistical error ★ So far, when reporting the error of an estimator 𝜃, we have quoted its standard deviation 𝜎𝜃 or an estimate of it 𝜎𝜃

★ This actually might work as definition of error: it provides a measure of how much the estimators would be spread should we repeat the same experiments many times,

★ But knowing the mean and standard deviation of a pdf is not all there is to know about the pdf, unless the functional form (e.g. Gaussian, Poisson, exponential, …) is known as well

★ In practice many estimators have a Gaussian pdf in the large sample limit (n→∞), so knowing mean and standard deviation of the estimator IS in may cases all there is to know about the sampling distribution

★ The standard deviation 𝜎𝜃 could in principle be what we mean by statistical error, but in the general case (in particular, if g(𝜃;𝜃) is not Gaussian) it is not

★ By convention, instead, what we quote as statistical errors are values 𝛥𝜃1 and 𝛥𝜃2 for which the confidence interval [𝜃-𝛥𝜃1,𝜃-𝛥𝜃2] has a well defined confidence level CL or coverage, that is, should the experiment be repeated a high number of times, the confidence intervals computed using that same rule, would contain the true value a fraction CL of the experiments.

50

^ ^ ^^

^ ^

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Classical confidence intervals (Neyman construction)

★ Suppose we have n observations of a random variable x, used to compute an estimator 𝜃(x1,…,xn) for a parameter 𝜃, obtaining the value 𝜃obs

★ In addition, we know (by calibration process or Monte Carlo study) the pdf of 𝜃, g(𝜃;𝜃) as a function of 𝜃 (we don’t know the true value of 𝜃)

51

^ ^

^

120 Statistical errors, confidence intervals and limits

0.5

o o 2

5

4

3

2

o o 2

9

3 4

b

3 4

5

5

Fig. 9.1 A p.d.f. g(8; 8) for an esti- mator e for a given value of the true parameter 8. The two shaded regions indicate the values of 8 :::; v{3, which has a probability {J, and e u a , which has a probabilit:y a.

Fig. 9.2 Construction of the confi- dence interval [a, b] given an observed value 80bs of the estimator 8 for the parameter 8 (see text).

Figure 9.2 shows an example of how the functions ua(O) and vf3(O) might appear as a function of the true value of O. The region between the two curves is called the confidence belt. The probability for the estimator to be inside the belt, regardless of the value of 0, is given by

P(Vf3(O) 0 ua(O)) = 1 - ex - (3. (9.3)

As long as ua( 0) and vf3 (0) are monotonically increasing functions of 0, which in general should be the case if (; is to be a good estimator for 0, one can determine the inverse functions

The ineqllalities

a(O) == b(O) == v;! (0).

(9.4)

★ For a given 𝜃, we can find the value u𝛼 such that integrating g above u𝛼 provides a given fixed probability 𝛼 of our choice, i.e.:

★ Similarly, we can find the value v𝛽 such that integrating g until v𝛽 provides a given fixed probability 𝛽 of our choice, i.e.:

^

P (✓̂ � u↵(✓)) =

Z 1

u↵(✓) g(✓̂; ✓)d✓̂ = ↵ (1)

Finally, we compute the radius ⌅ of the MAGIC e↵ective field of view, defined such that observations of an isotropic gamma-ray flux with a hypothetical instrument with a flat-top acceptance R(⇠) = R(0) for ⇠ < ⌅, and R(⇠) = 0 for ⇠ > ⌅, would yield the same number of detected gamma rays as with MAGIC, when no cuts on the arrival direction are applied. We can therefore obtain ⌅ from the condition

R ⌅ 0 2⇡ ⇠R(0) d⇠ =

R1 0 2⇡ ⇠R(⇠) d⇠, where R(⇠) is shown

in Figure 19-bottom, yielding ⌅ = 1�. We note, however, that standard observations of sources with an extension larger than 0.4� are technically di�cult, as in that case the edge of the source would fall into the background estimation region. Nevertheless, the e↵ective field of view is a useful quantity for non-standard observations of di↵use signals like, e.g. the cosmic electron flux [Borla-Tridon ICRC, HESS electron spectrum paper].

The di↵erences between our results for VERITAS obtained using the conventional likelihood (Eq. 2.1, Figure 8a) and those published by the Collaboration [16] (shown in Figures 8b-d) are due to the di↵erent Aeff assumed, since all the remaining values relevant for the computation of the exclusion contours (TOBS , J and upper limit to the number of signal events) are taken from Ref. [16]. In this work, we have used the values of Aeff reported by Wood [ref] for Segue observation/analysis. In addition, we have checked that the Improvement Factors that we obtain do not di↵er significantly (< 5%) if we use the values of Aeff reported by McCutcheon [ref] (assumed for the analysis of the globular cluster W13). From this, we infer the validity of the obtained Improvement Factors also for the Aeff actually used in Ref. [16].

In addition, according to these results, the sensitivity gain of the CTA with respect to VERITAS would be marginal, or even nonexistent, for certain annihilation channels and mass ranges. We have traced this inconsistency down to a probable overestimation of the VERI- TAS performance assumed in [16]. For that, we have used the response functions assumed for VERITAS, MAGIC and CTA to compute the integral sensitivity (5� significance in 50 hours of observations) for a Crab-like spectrum1 at the analysis threshold, for the di↵erent in- struments. The results obtained for MAGIC (1.3% of Crab flux above 110 GeV) and CTA (0.30% of Crab flux above 75 GeV) are consistent with those published by the respective collaborations (Refs. [MAGIC] and [CTA] respectively). On the other hand, VERITAS re- sults imply a sensitivity of 0.32% of Crab flux above 165 GeV, more than a factor 2 bet- ter than what reported by the collaboration [ref: http://veritas.sao.arizona.edu/about-veritas- mainmenu-81/veritas-specifications-mainmenu-111]. It must also be stressed that this result is obtained assuming the e↵ective areas as reported by Wood [ref], which for most of the decay channels and mass ranges produce slightly less constraining results than what published by VERITAS [16]. Therefore we expect that repeating this calculation with the Aeff actually used in [16] would produce an even larger disagreement between the implied and reference sensitivities to the Crab Nebula spectrum.

d�

dE (⌦) =

d�PP

dE ⇥ J(⌦) (2)

d�PP

dE = A

dN�

dE (3)

1dN/dE = 5.8�13(E/300GeV)�2.32�0.13 log10(E/300GeV) GeV�1 cm�2 s�1 [ref MAGIC performance]

1

Z 1

✓̂obs

g(✓̂; a) d✓̂ = ↵

Z ✓̂obs

�1 g(✓̂; a) d✓̂ = � (1)

✓ = ✓̂obs+(a�✓̂obs)

�(✓̂obs�b) (2)

P (a(✓̂)  ✓  b(✓̂)) = 1� ↵� � (3)

P (a(✓̂) � ✓) = ↵

P (b(✓̂)  ✓) = �

a(✓̂) � ✓

b(✓̂)  ✓

✓̂ � u↵(✓)

✓̂  v�(✓)

a(✓̂) ⌘ u�1 ↵ (✓̂)

b(✓̂) ⌘ v�1 �

(✓̂)

P (v�(✓) � ✓̂ � u↵(✓)) = 1� ↵� � (4)

P (✓̂  v�(✓)) =

Z v�(✓)

�1 g(✓̂; ✓)d✓̂ = � (5)

P (✓̂ � u↵(✓)) =

Z 1

u↵(✓) g(✓̂; ✓)d✓̂ = ↵ (6)

Finally, we compute the radius ⌅ of the MAGIC e↵ective field of view, defined such that observations of an isotropic gamma-ray flux with a hypothetical instrument with a flat-top acceptance R(⇠) = R(0) for ⇠ < ⌅, and R(⇠) = 0 for ⇠ > ⌅, would yield the same number of

1

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Classical confidence intervals (Neyman construction)

★ If we know g(𝜃;𝜃) we can compute the values of u𝛼(𝜃) and v𝛽(𝜃) for any value of 𝜃 (see the figure)

★ The region between the two curves is called the confidence belt.

★ Independently on the value of 𝜃, the probability for the estimator 𝜃 to be inside the belt is given by

★ As long as u𝛼(𝜃) and v𝛽(𝜃) are monotonically increasing functions of 𝜃 (true in general for good estimators), one can determine the inverse functions

52

120 Statistical errors, confidence intervals and limits

0.5

o o 2

5

4

3

2

o o 2

9

3 4

b

3 4

5

5

Fig. 9.1 A p.d.f. g(8; 8) for an esti- mator e for a given value of the true parameter 8. The two shaded regions indicate the values of 8 :::; v{3, which has a probability {J, and e u a , which has a probabilit:y a.

Fig. 9.2 Construction of the confi- dence interval [a, b] given an observed value 80bs of the estimator 8 for the parameter 8 (see text).

Figure 9.2 shows an example of how the functions ua(O) and vf3(O) might appear as a function of the true value of O. The region between the two curves is called the confidence belt. The probability for the estimator to be inside the belt, regardless of the value of 0, is given by

P(Vf3(O) 0 ua(O)) = 1 - ex - (3. (9.3)

As long as ua( 0) and vf3 (0) are monotonically increasing functions of 0, which in general should be the case if (; is to be a good estimator for 0, one can determine the inverse functions

The ineqllalities

a(O) == b(O) == v;! (0).

(9.4)

^

^

P (v�(✓) � ✓̂ � u↵(✓)) = 1� ↵� � (1)

P (✓̂  v↵(✓)) =

Z v↵(✓)

1 g(✓̂; ✓)d✓̂ = � (2)

P (✓̂ � u↵(✓)) =

Z 1

u↵(✓) g(✓̂; ✓)d✓̂ = ↵ (3)

Finally, we compute the radius ⌅ of the MAGIC e↵ective field of view, defined such that observations of an isotropic gamma-ray flux with a hypothetical instrument with a flat-top acceptance R(⇠) = R(0) for ⇠ < ⌅, and R(⇠) = 0 for ⇠ > ⌅, would yield the same number of detected gamma rays as with MAGIC, when no cuts on the arrival direction are applied. We can therefore obtain ⌅ from the condition

R ⌅ 0 2⇡ ⇠R(0) d⇠ =

R1 0 2⇡ ⇠R(⇠) d⇠, where R(⇠) is shown

in Figure 19-bottom, yielding ⌅ = 1�. We note, however, that standard observations of sources with an extension larger than 0.4� are technically di�cult, as in that case the edge of the source would fall into the background estimation region. Nevertheless, the e↵ective field of view is a useful quantity for non-standard observations of di↵use signals like, e.g. the cosmic electron flux [Borla-Tridon ICRC, HESS electron spectrum paper].

The di↵erences between our results for VERITAS obtained using the conventional likelihood (Eq. 2.1, Figure 8a) and those published by the Collaboration [16] (shown in Figures 8b-d) are due to the di↵erent Aeff assumed, since all the remaining values relevant for the computation of the exclusion contours (TOBS , J and upper limit to the number of signal events) are taken from Ref. [16]. In this work, we have used the values of Aeff reported by Wood [ref] for Segue observation/analysis. In addition, we have checked that the Improvement Factors that we obtain do not di↵er significantly (< 5%) if we use the values of Aeff reported by McCutcheon [ref] (assumed for the analysis of the globular cluster W13). From this, we infer the validity of the obtained Improvement Factors also for the Aeff actually used in Ref. [16].

In addition, according to these results, the sensitivity gain of the CTA with respect to VERITAS would be marginal, or even nonexistent, for certain annihilation channels and mass ranges. We have traced this inconsistency down to a probable overestimation of the VERI- TAS performance assumed in [16]. For that, we have used the response functions assumed for VERITAS, MAGIC and CTA to compute the integral sensitivity (5� significance in 50 hours of observations) for a Crab-like spectrum1 at the analysis threshold, for the di↵erent in- struments. The results obtained for MAGIC (1.3% of Crab flux above 110 GeV) and CTA (0.30% of Crab flux above 75 GeV) are consistent with those published by the respective collaborations (Refs. [MAGIC] and [CTA] respectively). On the other hand, VERITAS re- sults imply a sensitivity of 0.32% of Crab flux above 165 GeV, more than a factor 2 bet- ter than what reported by the collaboration [ref: http://veritas.sao.arizona.edu/about-veritas- mainmenu-81/veritas-specifications-mainmenu-111]. It must also be stressed that this result is obtained assuming the e↵ective areas as reported by Wood [ref], which for most of the decay channels and mass ranges produce slightly less constraining results than what published by VERITAS [16]. Therefore we expect that repeating this calculation with the Aeff actually used in [16] would produce an even larger disagreement between the implied and reference sensitivities to the Crab Nebula spectrum.

1dN/dE = 5.8�13(E/300GeV)�2.32�0.13 log10(E/300GeV) GeV�1 cm�2 s�1 [ref MAGIC performance]

1

a(✓̂) ⌘ u�1 ↵ (✓̂)

b(✓̂) ⌘ v�1 �

(✓̂)

P (v�(✓) � ✓̂ � u↵(✓)) = 1� ↵� � (1)

P (✓̂  v↵(✓)) =

Z v↵(✓)

1 g(✓̂; ✓)d✓̂ = � (2)

P (✓̂ � u↵(✓)) =

Z 1

u↵(✓) g(✓̂; ✓)d✓̂ = ↵ (3)

Finally, we compute the radius ⌅ of the MAGIC e↵ective field of view, defined such that observations of an isotropic gamma-ray flux with a hypothetical instrument with a flat-top acceptance R(⇠) = R(0) for ⇠ < ⌅, and R(⇠) = 0 for ⇠ > ⌅, would yield the same number of detected gamma rays as with MAGIC, when no cuts on the arrival direction are applied. We can therefore obtain ⌅ from the condition

R ⌅ 0 2⇡ ⇠R(0) d⇠ =

R1 0 2⇡ ⇠R(⇠) d⇠, where R(⇠) is shown

in Figure 19-bottom, yielding ⌅ = 1�. We note, however, that standard observations of sources with an extension larger than 0.4� are technically di�cult, as in that case the edge of the source would fall into the background estimation region. Nevertheless, the e↵ective field of view is a useful quantity for non-standard observations of di↵use signals like, e.g. the cosmic electron flux [Borla-Tridon ICRC, HESS electron spectrum paper].

The di↵erences between our results for VERITAS obtained using the conventional likelihood (Eq. 2.1, Figure 8a) and those published by the Collaboration [16] (shown in Figures 8b-d) are due to the di↵erent Aeff assumed, since all the remaining values relevant for the computation of the exclusion contours (TOBS , J and upper limit to the number of signal events) are taken from Ref. [16]. In this work, we have used the values of Aeff reported by Wood [ref] for Segue observation/analysis. In addition, we have checked that the Improvement Factors that we obtain do not di↵er significantly (< 5%) if we use the values of Aeff reported by McCutcheon [ref] (assumed for the analysis of the globular cluster W13). From this, we infer the validity of the obtained Improvement Factors also for the Aeff actually used in Ref. [16].

In addition, according to these results, the sensitivity gain of the CTA with respect to VERITAS would be marginal, or even nonexistent, for certain annihilation channels and mass ranges. We have traced this inconsistency down to a probable overestimation of the VERI- TAS performance assumed in [16]. For that, we have used the response functions assumed for VERITAS, MAGIC and CTA to compute the integral sensitivity (5� significance in 50 hours of observations) for a Crab-like spectrum1 at the analysis threshold, for the di↵erent in- struments. The results obtained for MAGIC (1.3% of Crab flux above 110 GeV) and CTA (0.30% of Crab flux above 75 GeV) are consistent with those published by the respective collaborations (Refs. [MAGIC] and [CTA] respectively). On the other hand, VERITAS re- sults imply a sensitivity of 0.32% of Crab flux above 165 GeV, more than a factor 2 bet- ter than what reported by the collaboration [ref: http://veritas.sao.arizona.edu/about-veritas- mainmenu-81/veritas-specifications-mainmenu-111]. It must also be stressed that this result is

1dN/dE = 5.8�13(E/300GeV)�2.32�0.13 log10(E/300GeV) GeV�1 cm�2 s�1 [ref MAGIC performance]

1

★ That is, e.g: given a value of 𝜃, a(𝜃) will give back the 𝜃 for which that 𝜃 is the value from which we should integrate g(𝜃;𝜃) in order to get probability 𝛼

^ ^ ^

^

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Classical confidence intervals (Neyman construction)

★ Then

or, taken together:

★ If a(𝜃) and b(𝜃) are evaluated for the observed value of the estimator 𝜃obs, that determines a pair of values a and b. The interval [a,b] is called confidence interval with confidence level or coverage CL=1-𝛼-𝛽

★ Should the experiment be repeated many times, we would obtain many values of [a,b] which would contain the real value for a fraction 1-𝛼-𝛽 of them

★ In practice a (b) is the value of the true parameter for which a fraction 𝛼 (𝛽) of repeated experiments would be higher (lower) than 𝜃obs, i.e. are determined by solving (probably numerically):

53

✓̂ � u↵(✓)

✓̂  v�(✓)

(1)

a(✓̂) ⌘ u�1 ↵ (✓̂)

b(✓̂) ⌘ v�1 �

(✓̂)

P (v�(✓) � ✓̂ � u↵(✓)) = 1� ↵� � (2)

P (✓̂  v↵(✓)) =

Z v↵(✓)

1 g(✓̂; ✓)d✓̂ = � (3)

P (✓̂ � u↵(✓)) =

Z 1

u↵(✓) g(✓̂; ✓)d✓̂ = ↵ (4)

Finally, we compute the radius ⌅ of the MAGIC e↵ective field of view, defined such that observations of an isotropic gamma-ray flux with a hypothetical instrument with a flat-top acceptance R(⇠) = R(0) for ⇠ < ⌅, and R(⇠) = 0 for ⇠ > ⌅, would yield the same number of detected gamma rays as with MAGIC, when no cuts on the arrival direction are applied. We can therefore obtain ⌅ from the condition

R ⌅ 0 2⇡ ⇠R(0) d⇠ =

R1 0 2⇡ ⇠R(⇠) d⇠, where R(⇠) is shown

in Figure 19-bottom, yielding ⌅ = 1�. We note, however, that standard observations of sources with an extension larger than 0.4� are technically di�cult, as in that case the edge of the source would fall into the background estimation region. Nevertheless, the e↵ective field of view is a useful quantity for non-standard observations of di↵use signals like, e.g. the cosmic electron flux [Borla-Tridon ICRC, HESS electron spectrum paper].

The di↵erences between our results for VERITAS obtained using the conventional likelihood (Eq. 2.1, Figure 8a) and those published by the Collaboration [16] (shown in Figures 8b-d) are due to the di↵erent Aeff assumed, since all the remaining values relevant for the computation of the exclusion contours (TOBS , J and upper limit to the number of signal events) are taken from Ref. [16]. In this work, we have used the values of Aeff reported by Wood [ref] for Segue observation/analysis. In addition, we have checked that the Improvement Factors that we obtain do not di↵er significantly (< 5%) if we use the values of Aeff reported by McCutcheon [ref] (assumed for the analysis of the globular cluster W13). From this, we infer the validity of the obtained Improvement Factors also for the Aeff actually used in Ref. [16].

In addition, according to these results, the sensitivity gain of the CTA with respect to VERITAS would be marginal, or even nonexistent, for certain annihilation channels and mass ranges. We have traced this inconsistency down to a probable overestimation of the VERI- TAS performance assumed in [16]. For that, we have used the response functions assumed for VERITAS, MAGIC and CTA to compute the integral sensitivity (5� significance in 50

1

P (a(✓̂) � ✓) = ↵

P (b(✓̂)  ✓) = �

a(✓̂) � ✓

b(✓̂)  ✓

✓̂ � u↵(✓)

✓̂  v�(✓)

a(✓̂) ⌘ u�1 ↵ (✓̂)

b(✓̂) ⌘ v�1 �

(✓̂)

P (v�(✓) � ✓̂ � u↵(✓)) = 1� ↵� � (1)

P (✓̂  v↵(✓)) =

Z v↵(✓)

1 g(✓̂; ✓)d✓̂ = � (2)

P (✓̂ � u↵(✓)) =

Z 1

u↵(✓) g(✓̂; ✓)d✓̂ = ↵ (3)

Finally, we compute the radius ⌅ of the MAGIC e↵ective field of view, defined such that observations of an isotropic gamma-ray flux with a hypothetical instrument with a flat-top acceptance R(⇠) = R(0) for ⇠ < ⌅, and R(⇠) = 0 for ⇠ > ⌅, would yield the same number of detected gamma rays as with MAGIC, when no cuts on the arrival direction are applied. We can therefore obtain ⌅ from the condition

R ⌅ 0 2⇡ ⇠R(0) d⇠ =

R1 0 2⇡ ⇠R(⇠) d⇠, where R(⇠) is shown

in Figure 19-bottom, yielding ⌅ = 1�. We note, however, that standard observations of sources with an extension larger than 0.4� are technically di�cult, as in that case the edge of the source would fall into the background estimation region. Nevertheless, the e↵ective field of view is a useful quantity for non-standard observations of di↵use signals like, e.g. the cosmic electron flux [Borla-Tridon ICRC, HESS electron spectrum paper].

The di↵erences between our results for VERITAS obtained using the conventional likelihood (Eq. 2.1, Figure 8a) and those published by the Collaboration [16] (shown in Figures 8b-d) are due to the di↵erent Aeff assumed, since all the remaining values relevant for the computation of the exclusion contours (TOBS , J and upper limit to the number of signal events) are taken from Ref. [16]. In this work, we have used the values of Aeff reported by Wood [ref] for Segue observation/analysis. In addition, we have checked that the Improvement Factors that we obtain do not di↵er significantly (< 5%) if we use the values of Aeff reported by McCutcheon [ref]

1

P (a(✓̂) � ✓) = ↵

P (b(✓̂)  ✓) = �

a(✓̂) � ✓

b(✓̂)  ✓

✓̂ � u↵(✓)

✓̂  v�(✓)

a(✓̂) ⌘ u�1 ↵ (✓̂)

b(✓̂) ⌘ v�1 �

(✓̂)

P (v�(✓) � ✓̂ � u↵(✓)) = 1� ↵� � (1)

P (✓̂  v↵(✓)) =

Z v↵(✓)

1 g(✓̂; ✓)d✓̂ = � (2)

P (✓̂ � u↵(✓)) =

Z 1

u↵(✓) g(✓̂; ✓)d✓̂ = ↵ (3)

Finally, we compute the radius ⌅ of the MAGIC e↵ective field of view, defined such that observations of an isotropic gamma-ray flux with a hypothetical instrument with a flat-top acceptance R(⇠) = R(0) for ⇠ < ⌅, and R(⇠) = 0 for ⇠ > ⌅, would yield the same number of detected gamma rays as with MAGIC, when no cuts on the arrival direction are applied. We can therefore obtain ⌅ from the condition

R ⌅ 0 2⇡ ⇠R(0) d⇠ =

R1 0 2⇡ ⇠R(⇠) d⇠, where R(⇠) is shown

in Figure 19-bottom, yielding ⌅ = 1�. We note, however, that standard observations of sources with an extension larger than 0.4� are technically di�cult, as in that case the edge of the source would fall into the background estimation region. Nevertheless, the e↵ective field of view is a useful quantity for non-standard observations of di↵use signals like, e.g. the cosmic electron flux [Borla-Tridon ICRC, HESS electron spectrum paper].

The di↵erences between our results for VERITAS obtained using the conventional likelihood (Eq. 2.1, Figure 8a) and those published by the Collaboration [16] (shown in Figures 8b-d) are due to the di↵erent Aeff assumed, since all the remaining values relevant for the computation of the exclusion contours (TOBS , J and upper limit to the number of signal events) are taken from Ref. [16]. In this work, we have used the values of Aeff reported by Wood [ref] for Segue observation/analysis. In addition, we have checked that the Improvement Factors that we obtain do not di↵er significantly (< 5%) if we use the values of Aeff reported by McCutcheon [ref]

1

P (a(✓̂)  ✓  b(✓̂)) = 1� ↵� � (1)

P (a(✓̂) � ✓) = ↵

P (b(✓̂)  ✓) = �

a(✓̂) � ✓

b(✓̂)  ✓

✓̂ � u↵(✓)

✓̂  v�(✓)

a(✓̂) ⌘ u�1 ↵ (✓̂)

b(✓̂) ⌘ v�1 �

(✓̂)

P (v�(✓) � ✓̂ � u↵(✓)) = 1� ↵� � (2)

P (✓̂  v↵(✓)) =

Z v↵(✓)

1 g(✓̂; ✓)d✓̂ = � (3)

P (✓̂ � u↵(✓)) =

Z 1

u↵(✓) g(✓̂; ✓)d✓̂ = ↵ (4)

Finally, we compute the radius ⌅ of the MAGIC e↵ective field of view, defined such that observations of an isotropic gamma-ray flux with a hypothetical instrument with a flat-top acceptance R(⇠) = R(0) for ⇠ < ⌅, and R(⇠) = 0 for ⇠ > ⌅, would yield the same number of detected gamma rays as with MAGIC, when no cuts on the arrival direction are applied. We can therefore obtain ⌅ from the condition

R ⌅ 0 2⇡ ⇠R(0) d⇠ =

R1 0 2⇡ ⇠R(⇠) d⇠, where R(⇠) is shown

in Figure 19-bottom, yielding ⌅ = 1�. We note, however, that standard observations of sources with an extension larger than 0.4� are technically di�cult, as in that case the edge of the source would fall into the background estimation region. Nevertheless, the e↵ective field of view is a useful quantity for non-standard observations of di↵use signals like, e.g. the cosmic electron flux [Borla-Tridon ICRC, HESS electron spectrum paper].

The di↵erences between our results for VERITAS obtained using the conventional likelihood (Eq. 2.1, Figure 8a) and those published by the Collaboration [16] (shown in Figures 8b-d) are due to the di↵erent Aeff assumed, since all the remaining values relevant for the computation

1

^ ^ ^

Z 1

✓̂obs

g(✓̂; a) d✓̂ = ↵

Z ✓̂obs

�1 g(✓̂; b) d✓̂ = � (1)

✓ = ✓̂obs+(a�✓̂obs)

�(✓̂obs�b) (2)

P (a(✓̂)  ✓  b(✓̂)) = 1� ↵� � (3)

P (a(✓̂) � ✓) = ↵

P (b(✓̂)  ✓) = �

a(✓̂) � ✓

b(✓̂)  ✓

✓̂ � u↵(✓)

✓̂  v�(✓)

a(✓̂) ⌘ u�1 ↵ (✓̂)

b(✓̂) ⌘ v�1 �

(✓̂)

P (v�(✓) � ✓̂ � u↵(✓)) = 1� ↵� � (4)

P (✓̂  v�(✓)) =

Z v�(✓)

�1 g(✓̂; ✓)d✓̂ = � (5)

P (✓̂ � u↵(✓)) =

Z 1

u↵(✓) g(✓̂; ✓)d✓̂ = ↵ (6)

Finally, we compute the radius ⌅ of the MAGIC e↵ective field of view, defined such that observations of an isotropic gamma-ray flux with a hypothetical instrument with a flat-top acceptance R(⇠) = R(0) for ⇠ < ⌅, and R(⇠) = 0 for ⇠ > ⌅, would yield the same number of

1

122 Statistical errors, confidence intervals and limits

one actually obtained, Bobs, as is illustrated in Fig. 9.3. Similarly, b is the value of () for which a fraction {3 of the estimates would be lower than Bobs. That is, taking eobs = ua(a) = v,6(b), equations (9.1) and (9.2) become

(9.9)

{3

The previously described procedure to determine the confidence interval is thus equivalent to solving (9.9) for a and b, e.g. numerically.

(a)

0.5

o o 2 3 4 5 9 Fig. 9.3 (a) The p.d.f. g(B; a), where

a is the lower limit of the confidence

b (b) interval. If the true parameter B were

equal to a, the estimates 0 would be &reater than the one actually observed Bobs with a probability Q. (b) The

0.5 p.d.f. g(O; b), where b is the upper limit of the confidence interval. If B were equal to b, B would be observed less

2 3 4 5 than Bobs with probability {3. o o 9

Figure 9.3 also illustrates the relationship between a confidence interval and a test of goodness-of-fit, cf. Section 4.5. For example, we could test the hypothesis () == a using B as a test statistic. If we define the region e Bobs as having equal or less agreement with the hypothesis than the result obtained (a one-sided test), then the resulting P-value of the test is a. For the confidence interval, however, the probability a is specified first, and the value a is a random quantity depending on the data. For a goodness-of-fit test, the hypothesis, here () = a, is specified and the P-value is treated as a random variable.

Note that one sometimes calls the P-value, here equal to a, the 'confidence level' of the test, whereas the one-sided confidence interval () a has a confidence level of 1 - a. That is, for a test, small a indicates a low level of confidence in the hypothesis () = a. For a confidence interval, small a indicates a high level of

122 Statistical errors, confidence intervals and limits

one actually obtained, Bobs, as is illustrated in Fig. 9.3. Similarly, b is the value of () for which a fraction {3 of the estimates would be lower than Bobs. That is, taking eobs = ua(a) = v,6(b), equations (9.1) and (9.2) become

(9.9)

{3

The previously described procedure to determine the confidence interval is thus equivalent to solving (9.9) for a and b, e.g. numerically.

(a)

0.5

o o 2 3 4 5 9 Fig. 9.3 (a) The p.d.f. g(B; a), where

a is the lower limit of the confidence

b (b) interval. If the true parameter B were

equal to a, the estimates 0 would be &reater than the one actually observed Bobs with a probability Q. (b) The

0.5 p.d.f. g(O; b), where b is the upper limit of the confidence interval. If B were equal to b, B would be observed less

2 3 4 5 than Bobs with probability {3. o o 9

Figure 9.3 also illustrates the relationship between a confidence interval and a test of goodness-of-fit, cf. Section 4.5. For example, we could test the hypothesis () == a using B as a test statistic. If we define the region e Bobs as having equal or less agreement with the hypothesis than the result obtained (a one-sided test), then the resulting P-value of the test is a. For the confidence interval, however, the probability a is specified first, and the value a is a random quantity depending on the data. For a goodness-of-fit test, the hypothesis, here () = a, is specified and the P-value is treated as a random variable.

Note that one sometimes calls the P-value, here equal to a, the 'confidence level' of the test, whereas the one-sided confidence interval () a has a confidence level of 1 - a. That is, for a test, small a indicates a low level of confidence in the hypothesis () = a. For a confidence interval, small a indicates a high level of

^

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Classical confidence intervals (Neyman construction)

★ Given a CL, there is a freedom on how we “share” the probability 1- CL between 𝛼 and 𝛽: ✦ 𝛼 = 𝛽 ≡ 𝛾/2 = (1-CL)/2 ⇒ defines central confidence intervals which

does not imply that a and b are equidistant from 𝜃obs, i.e. the confidence interval can be asymmetric wrt 𝜃obs, normally quoted as

or, in case of equidistance (a-𝜃obs=𝜃obs-b≡𝛥𝜃), we quote it as

✦ 𝛼 = 1-CL; 𝛽 = 0 ⇒ defines an upper limit

✦ 𝛼 = 0; 𝛽 = 1-CL ⇒ defines a lower limit

★ When represented graphically the quantities a-𝜃obs and 𝜃obs-b are shown as error bars around 𝜃obs

54

✓ = ✓̂obs ±�✓̂ (1)

⌫i(✓) = ⌫tot

Z x max i

xmin i

f(x;✓)dx (2)

logL(⌫tot,✓) = �⌫tot + NX

i=1

ni log ⌫i(⌫tot,✓) (3)

f(n;⌫) = ⌫ntot tot e�⌫tot

ntot!

ntot!

n1!...nN !

✓ ⌫1 ⌫tot

◆n1

...

✓ ⌫N ⌫tot

◆nN

= NY

i=1

⌫ni i

ni! e�⌫ (4)

logL(✓) = NX

i=1

ni log ⌫i(✓) (5)

f(n;⌫) = ntot!

n1!...nN !

✓ ⌫1 ntot

◆n1

...

✓ ⌫N ntot

◆nN

(6)

⌫i(✓) = ntot

Z x max i

xmin i

f(x;✓)dx (7)

logL(✓̂ ±c�2 ✓̂ ) = logLmax �

1

2 (8)

logL(✓) = logLmax � (✓ � ✓̂)2

2c�2 ✓̂

(9)

logL(✓) = logL(✓̂) +

 @ logL

@✓

�

✓=✓̂

(✓ � ✓̂) + 1

2!

 @2 logL

@✓2

�

✓=✓̂

(✓ � ✓̂)2 + ... (10)

c�2 ✓̂ =

1

@2 logL @✓2

!����� ✓=✓̂

(11)

(V �1)ij = E

 �@2 logL

@✓i@✓j

� (12)

1

^ ^ ^

^ ^

^ ^

↵ =

Z 1

✓̂obs

1p 2⇡�

✓̂

e � 1

2 ( ✓̂�a � ✓̂ )2

d✓̂ =

( x = ✓̂�a

�✓̂

dx = 1 �✓̂ d✓̂

) =

Z 1

✓̂obs�a � ✓̂

1p 2⇡

e�x 2 /2 dx = 1� �(

✓̂obs � a

� ✓̂

)

� =

Z ✓̂obs

�1

1p 2⇡�

✓̂

e � 1

2 ( ✓̂�b � ✓̂ )2

d✓̂ =

( x = ✓̂�b

�✓̂

dx = 1 �✓̂ d✓̂

) =

Z ✓̂obs�b � ✓̂

�1

1p 2⇡

e�x 2 /2 dx = �(

✓̂obs � b

� ✓̂

)

✓ = ✓̂obs +(a� ✓̂obs) �(✓̂obs � b)

(1)

✓ = ✓̂obs ±�✓̂ (2)

⌫i(✓) = ⌫tot

Z x max i

xmin i

f(x;✓)dx (3)

logL(⌫tot,✓) = �⌫tot + NX

i=1

ni log ⌫i(⌫tot,✓) (4)

f(n;⌫) = ⌫ntot tot e�⌫tot

ntot!

ntot!

n1!...nN !

✓ ⌫1 ⌫tot

◆n1

...

✓ ⌫N ⌫tot

◆nN

= NY

i=1

⌫ni i

ni! e�⌫ (5)

logL(✓) = NX

i=1

ni log ⌫i(✓) (6)

f(n;⌫) = ntot!

n1!...nN !

✓ ⌫1 ntot

◆n1

...

✓ ⌫N ntot

◆nN

(7)

⌫i(✓) = ntot

Z x max i

xmin i

f(x;✓)dx (8)

logL(✓̂ ±c�2 ✓̂ ) = logLmax �

1

2 (9)

logL(✓) = logLmax � (✓ � ✓̂)2

2c�2 ✓̂

(10)

1

^

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Confidence interval for Gaussian pdf ★ Consider pdf of the estimator 𝜃 is Gaussian (common situation, because of the CLT)

★ a and b are found by solving:

with ϕ the cumulative distribution of the standard Gaussian (𝜇=0, 𝜎=1), with the argument being the distance between 𝜃obs and a (or b) in units of 𝜎𝜃.

★ Then, we can use the inverse of ϕ (ϕ-1) to solve for a and b:

★ Values of ϕ-1(CL) or equivalently ϕ-1(1-𝛾) are tabulated, both for rounded values of ϕ-1 or the CL (see next)

★ As mentioned, in many practical cases we do not know 𝜎𝜃, and we use 𝜎𝜃 (valid if 𝜎𝜃 unbiased for n→∞)

55

^

↵ =

Z 1

✓̂obs

1p 2⇡�

✓̂

e � 1

2 ( ✓̂�a � ✓̂ )2

d✓̂ =

( x = ✓̂�a

�✓̂

dx = 1 �✓̂ d✓̂

) =

Z 1

✓̂obs�a � ✓̂

1p 2⇡

e�x 2 /2 dx = 1� �(

✓̂obs � a

� ✓̂

)

� =

Z ✓̂obs

�1

1p 2⇡�

✓̂

e � 1

2 ( ✓̂�b � ✓̂ )2

d✓̂ =

( x = ✓̂�b

�✓̂

dx = 1 �✓̂ d✓̂

) =

Z ✓̂obs�b � ✓̂

�1

1p 2⇡

e�x 2 /2 dx = �(

✓̂obs � b

� ✓̂

)

✓ = ✓̂obs +(a� ✓̂obs) �(✓̂obs � b)

(1)

✓ = ✓̂obs ±�✓̂ (2)

⌫i(✓) = ⌫tot

Z x max i

xmin i

f(x;✓)dx (3)

logL(⌫tot,✓) = �⌫tot + NX

i=1

ni log ⌫i(⌫tot,✓) (4)

f(n;⌫) = ⌫ntot tot e�⌫tot

ntot!

ntot!

n1!...nN !

✓ ⌫1 ⌫tot

◆n1

...

✓ ⌫N ⌫tot

◆nN

= NY

i=1

⌫ni i

ni! e�⌫ (5)

logL(✓) = NX

i=1

ni log ⌫i(✓) (6)

f(n;⌫) = ntot!

n1!...nN !

✓ ⌫1 ntot

◆n1

...

✓ ⌫N ntot

◆nN

(7)

⌫i(✓) = ntot

Z x max i

xmin i

f(x;✓)dx (8)

logL(✓̂ ±c�2 ✓̂ ) = logLmax �

1

2 (9)

logL(✓) = logLmax � (✓ � ✓̂)2

2c�2 ✓̂

(10)

1

a = ✓̂obs � � ✓̂ ��1(1� ↵)

b = ✓̂obs + � ✓̂ ��1(1� �) (1)

↵ =

Z 1

✓̂obs

1p 2⇡�

✓̂

e � 1

2 ( ✓̂�a � ✓̂ )2

d✓̂ =

( x = ✓̂�a

�✓̂

dx = 1 �✓̂ d✓̂

) =

Z 1

✓̂obs�a � ✓̂

1p 2⇡

e�x 2 /2 dx = 1� �(

✓̂obs � a

� ✓̂

)

� =

Z ✓̂obs

�1

1p 2⇡�

✓̂

e � 1

2 ( ✓̂�b � ✓̂ )2

d✓̂ =

( x = ✓̂�b

�✓̂

dx = 1 �✓̂ d✓̂

) =

Z ✓̂obs�b � ✓̂

�1

1p 2⇡

e�x 2 /2 dx = �(

✓̂obs � b

� ✓̂

)

✓ = ✓̂obs +(a� ✓̂obs) �(✓̂obs � b)

(2)

✓ = ✓̂obs ±�✓̂ (3)

⌫i(✓) = ⌫tot

Z x max i

xmin i

f(x;✓)dx (4)

logL(⌫tot,✓) = �⌫tot + NX

i=1

ni log ⌫i(⌫tot,✓) (5)

f(n;⌫) = ⌫ntot tot e�⌫tot

ntot!

ntot!

n1!...nN !

✓ ⌫1 ⌫tot

◆n1

...

✓ ⌫N ⌫tot

◆nN

= NY

i=1

⌫ni i

ni! e�⌫ (6)

logL(✓) = NX

i=1

ni log ⌫i(✓) (7)

f(n;⌫) = ntot!

n1!...nN !

✓ ⌫1 ntot

◆n1

...

✓ ⌫N ntot

◆nN

(8)

⌫i(✓) = ntot

Z x max i

xmin i

f(x;✓)dx (9)

logL(✓̂ ±c�2 ✓̂ ) = logLmax �

1

2 (10)

1

^

^

^^ ^ ^^

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Confidence interval for Gaussian pdf

★ By convention we report the 68.3% central confidence interval, i.e. 𝛼 = 𝛽 = 𝛾/ 2, with ϕ-1(1-𝛾/2) =1, i.e. a 1𝜎 error bar, which results in the simple prescription

, i.e. for Gaussian the error is the standard deviation

56

124 Statistical errors, confidence intervals and limits

3 0.6 3 0.6

S. (a) s. (b)

<1>-1 (Y/2) <1>-1 (1-y/2)

0.4 0.4

0.2 0.2 a

o o -4 -2 o 2 4 -4 -2 o 2 4

x x

Fig. 9.4 The standard Gaussian p.d.f. <p(x) showing the relationship between the quantiles and the confidence level for (a) a central confidence interval and (b) a one-sided confidence

interval.

somewhat complicated procedure explained in the previous section results in a simple prescription for determining the confidence interval.

Suppose that the standard deviation (J'§ is known, and that the experiment has resulted in an estimate Bobs. According to equations (9.9), the confidence interval [a, b] is determined by solving the equations

(9.11)

(3

for a and b, where G has been expressed using the cumulative distribution of the standard Gaussian 4> (2.26) (see also (2.27)). This gives

a = Bobs - (J'§4>-1(1- a), A -1 b=()obs+(J'§4> (1-(3).

(9.12)

Here 4>-1 is the inverse function of 4>, i.e. the quantile of the standard Gaussian, and in order to make the two equations symmetric we have used 4>-1 ((3) = _4>-1(1 - (3).

The quantiles 4>-1(1_ a) and 4>-1 (1- (3) represent how far away the interval limits a and b are located with respect to the estimate Bobs in units of the standard deviation (J'§. The relationship between the quantiles of the standard Gaussian distribution and the confidence level is illustrated in Fig. 9.4( a) for central and Fig. 9.4(b) for one-sided confidence intervals.

Confidence interval for a Gaussian distributed estimator 125

Consider a central confidence interval with a =.(3 = ,/2. The confidence level 1-, is often chosen such that the quantile is a small integer, e.g. cI>-1(1-,/2) = 1,2,3, .... Similarly, for one-sided intervals (limits) one often chooses a small integer for cI>-1 (1 - a). Commonly used values for both central and one-sided intervals are shown in Table 9.1. Alternatively one can choose a round number for the confidence level instead of for the quantile. Commonly used values are shown in Table 9.2. Other possible values can be obtained from [Bra92, Fr079 , Dud88] or from computer routines (e.g. the routine GAUSIN in [CER97]).

Table 9.1 The values of the confidence level for different values of the quantile of the standard Gaussian for central intervals (left) the quantile (1-,/2) and confidence level 1-,; for one-sided intervals (right) the quantile - Q) and confidence level 1- Q.

cI> (1 - ,/2) 1 2 3

0.6827 0.9544 0.9973

cI> (1-0') 1 2 3

1 - a 0.8413 0.9772 0.9987

Table 9.2 The values of the quantile of the standard Gaussian for different values of the confidence level: for central intervals (left) the confidence level 1 - , and the quan- tile (1 - ,/2); for one-sided intervals (right) the confidence level I - Q and the quantile

(I - Q).

0.90 0.95 0.99

-1.645 1.960 2.576

1 - a 0.90 0.95 0.99

1.282 1.645 2.326

For the conventional 68.3% central confidence interval one has a = {3 = ,/2, with cI>-1 (1-, /2) = 1, i.e. a' 1 (J' error bar'. This results in the simple prescription

(9.13)

Thus for the case of a Gaussian distributed estimator, the 68.3% central confi- dence interval is given by the estimated value plus or minus one standard de- viation. The final result of the measurement of () is then simply reported as Oobs±(J'o·

If the standard deviation (J'o is not known a priori but rather is estimated from the data, then the situation is in principle somewhat more complicated. If, for example, the estimated standard deviation (;-0 had been used instead of (J'o' then it would not have been so simple to relate the cumulative distribution G(e; (), (;-g) to cI>, the cumulative distribution of the standard Gaussian, since (;-{} depends in general on O. In practice, however, the recipe given above can still

Confidence interval for a Gaussian distributed estimator 125

Consider a central confidence interval with a =.(3 = ,/2. The confidence level 1-, is often chosen such that the quantile is a small integer, e.g. cI>-1(1-,/2) = 1,2,3, .... Similarly, for one-sided intervals (limits) one often chooses a small integer for cI>-1 (1 - a). Commonly used values for both central and one-sided intervals are shown in Table 9.1. Alternatively one can choose a round number for the confidence level instead of for the quantile. Commonly used values are shown in Table 9.2. Other possible values can be obtained from [Bra92, Fr079 , Dud88] or from computer routines (e.g. the routine GAUSIN in [CER97]).

Table 9.1 The values of the confidence level for different values of the quantile of the standard Gaussian for central intervals (left) the quantile (1-,/2) and confidence level 1-,; for one-sided intervals (right) the quantile - Q) and confidence level 1- Q.

cI> (1 - ,/2) 1 2 3

0.6827 0.9544 0.9973

cI> (1-0') 1 2 3

1 - a 0.8413 0.9772 0.9987

Table 9.2 The values of the quantile of the standard Gaussian for different values of the confidence level: for central intervals (left) the confidence level 1 - , and the quan- tile (1 - ,/2); for one-sided intervals (right) the confidence level I - Q and the quantile

(I - Q).

0.90 0.95 0.99

-1.645 1.960 2.576

1 - a 0.90 0.95 0.99

1.282 1.645 2.326

For the conventional 68.3% central confidence interval one has a = {3 = ,/2, with cI>-1 (1-, /2) = 1, i.e. a' 1 (J' error bar'. This results in the simple prescription

(9.13)

Thus for the case of a Gaussian distributed estimator, the 68.3% central confi- dence interval is given by the estimated value plus or minus one standard de- viation. The final result of the measurement of () is then simply reported as Oobs±(J'o·

If the standard deviation (J'o is not known a priori but rather is estimated from the data, then the situation is in principle somewhat more complicated. If, for example, the estimated standard deviation (;-0 had been used instead of (J'o' then it would not have been so simple to relate the cumulative distribution G(e; (), (;-g) to cI>, the cumulative distribution of the standard Gaussian, since (;-{} depends in general on O. In practice, however, the recipe given above can still

Central CI 1-sided UL

Rounded quantile (ϕ-1) Rounded CL

[a, b] = [✓̂obs � � ✓̂ , ✓̂obs + �

✓̂ ] (1)

a = ✓̂obs � � ✓̂ ��1(1� ↵)

b = ✓̂obs + � ✓̂ ��1(1� �) (2)

↵ =

Z 1

✓̂obs

1p 2⇡�

✓̂

e � 1

2 ( ✓̂�a � ✓̂ )2

d✓̂ =

( x = ✓̂�a

�✓̂

dx = 1 �✓̂ d✓̂

) =

Z 1

✓̂obs�a � ✓̂

1p 2⇡

e�x 2 /2 dx = 1� �(

✓̂obs � a

� ✓̂

)

� =

Z ✓̂obs

�1

1p 2⇡�

✓̂

e � 1

2 ( ✓̂�b � ✓̂ )2

d✓̂ =

( x = ✓̂�b

�✓̂

dx = 1 �✓̂ d✓̂

) =

Z ✓̂obs�b � ✓̂

�1

1p 2⇡

e�x 2 /2 dx = �(

✓̂obs � b

� ✓̂

)

✓ = ✓̂obs +(a� ✓̂obs) �(✓̂obs � b)

(3)

✓ = ✓̂obs ±�✓̂ (4)

⌫i(✓) = ⌫tot

Z x max i

xmin i

f(x;✓)dx (5)

logL(⌫tot,✓) = �⌫tot + NX

i=1

ni log ⌫i(⌫tot,✓) (6)

f(n;⌫) = ⌫ntot tot e�⌫tot

ntot!

ntot!

n1!...nN !

✓ ⌫1 ⌫tot

◆n1

...

✓ ⌫N ⌫tot

◆nN

= NY

i=1

⌫ni i

ni! e�⌫ (7)

logL(✓) = NX

i=1

ni log ⌫i(✓) (8)

f(n;⌫) = ntot!

n1!...nN !

✓ ⌫1 ntot

◆n1

...

✓ ⌫N ntot

◆nN

(9)

⌫i(✓) = ntot

Z x max i

xmin i

f(x;✓)dx (10)

1

Master of Interdisciplinary Research in Experimental Sciences J. Rico - Parameter Estimation UPF, Barcelona, September 2018

Exercises ★ Six independent observations from a Gaussian

distribution N(𝜇,𝜎) are given by x=(12.11, 9.83, 11.44, 7.52, 6.69, 10.73). If 𝜎=2.1 is known, find the symmetric confidence intervals for 𝜇 with confidence levels 1-𝛼=0.6827, 0.90, and 0.95, respectively. Discuss the case in which 𝜎 is unknown

57

ParEst_input_3.txt

# You can read this file with e.g. np.loadtxt # First column are x values # Second column are standard deviations 5.12 1.00 2.89 1.00 2.97 1.00 1.35 2.00 5.23 2.00 -1.10 2.00 4.37 0.50 2.97 0.70 3.98 1.50 3.13 1.50 5.69 1.50 1.44 1.00 3.18 1.00 3.12 1.00 4.41 0.80