Estimating parameters in spatio- temporal Quermass- in - - PowerPoint PPT Presentation

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

Estimating parameters in spatio-temporal Quermass-interaction process

Kateˇ rina Staˇ nkov´ a Helisov´ a

Czech Technical University in Prague helisova@math.feld.cvut.cz joint work with Viktor Beneˇ s and Mark´ eta Zikmundov´ a Charles University in Prague 9th May 2012

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

Outline

• 1. Quermass-interaction process and its extension
• 2. Simulation
• 3. Spatio-temporal Quermass-interaction process
• 4. Maximum likelihood method using MCMC
• 5. Particle filtering
• 6. Particle MCMC

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

Outline

• 1. Quermass-interaction process and its extension
• 2. Simulation
• 3. Spatio-temporal Quermass-interaction process
• 4. Maximum likelihood method using MCMC
• 5. Particle filtering
• 6. Particle MCMC

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

Notation

• x = b(u, r) ... a disc with centre in u ∈ R2 and radius r ∈ (0, ∞)
• x = {x1, . . . , xn} ... finite configuration of n discs
• Ux ... the union of discs from the configuration x.
• Y ... random disc Boolean model (i.e. union of discs without any

interactions) with an intensity function of discs centers ρ(u) and probability distribution of the discs radii Q

• X ... random disc process which is absolutely continuous with re-

spect to the process Y

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

Assumptions

• The intensity function ρ(u) = 1 on a bounded set S and ρ(u) = 0
• therwise, i.e. the centers of the reference Boolean model form unit

Poisson process on S.

• For any finite configuration of discs x = {x1, . . . , xn}, the proba-

bility measure of X with respect to the probability measure of Y is given by a density fθ(x) =exp{θ · T(Ux)} cθ , where – cθ is the normalizing constant, – θ is d-dimensional vector of parameters, – T = T(Ux) is a d-dimensional vector of geometrical characteris- tics of the union Ux of the discs from the configuration x.

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

Quermass-interaction process

The density is of the form fθ(x) = 1 cθ exp{θ1A(Ux) + θ2L(Ux) + θ3χ(Ux)}, where

• A = A(Ux) is the area,
• L = L(Ux) is the perimeter,
• χ = χ(Ux) is the Euler-Poincar´

e characteristic (the number of con- nected components minus the number of holes, i.e. χ(Ux) = Ncc(Ux) − Nh(Ux))

• f the union Ux.

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

Extended Quermass-interaction process

• Møller, Helisov´

a (2008): – In the density fθ(x) =exp{θ · T(Ux)} cθ , we have T = (A, L, χ, Nh, Nbv, Nid), where Nbv = Nbv(Ux) is the number of boundary vertices, Nid = Nid(Ux) is the number of isolated discs,

• f the union Ux.

– Theory and simulations studied.

• Møller, Helisov´

a (2010): – T = (A, L, Ncc, Nh). – Statistical analysis.

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

Outline

• 1. Quermass-interaction process and its extension
• 2. Simulation
• 3. Spatio-temporal Quermass-interaction process
• 4. Maximum likelihood method using MCMC
• 5. Particle filtering
• 6. Particle MCMC

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

Papangelou conditional intensity

Definition For a finite x ⊂ S × (0, ∞) and y ∈ S × (0, ∞) \ x, Papangelou conditional intensity is defined as λθ(x, y) = fθ(x ∪ {y})/fθ(x). Denoting A(x, y) =A(Ux∪y) − A(Ux), L(x, y) =L(Ux∪y) − L(Ux), χ(x, y) =χ(Ux∪y) − χ(Ux), we get λθ(x, y) = exp (θ1A(x, y) + θ2L(x, y) + θ3χ(x, y)) .

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

MCMC algorithm

• 1. Suppose that in iteration t, we have a configuration xt = {x1, . . . , xn}
• 2. Proposal in iteration t + 1:

(a) with probability 1/2, the proposal is xt ∪ {xn+1}

• i. we accept the proposal with probability min{1; H(xt, xn+1)}

and set xt+1 = xt ∪ {xn+1}

• ii. else we set xt+1 = xt

(b) else, the proposal is xt\{xi}

• i. we accept the proposal with probability min{1; 1/H(xt\{xi}, xi)}

and set xt+1 = xt\{xi}

• ii. else xt+1 = xt

where H(xt, xn+1) = λθ(xt, xn+1) |S|

n+1

and H(xt\{xi}, xi) = λθ(xt\{xi}, xi)|S|

n

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

Outline

• 1. Quermass-interaction process and its extension
• 2. Simulation
• 3. Spatio-temporal Quermass-interaction process
• 4. Maximum likelihood method using MCMC
• 5. Particle filtering
• 6. Particle MCMC

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

Spatio-temporal Quermass-interaction process

Zikmundov´ a, Staˇ nkov´ a Helisov´ a, Beneˇ s (2012): fθ(k)(x) = exp{θ(k) · T(Ux)} cθ(k) , where θ(k) = θ(k−1) + η(k), k = 1, 2 . . . , T, where θ(0) fixed is given and η(k) are iid random vectors with Gaussian distribution N(a, σ2I), where a ∈ Rd (d = 3 for basic Quermass- interactio process and d = 4, 5, 6 for extended versions), σ2 > 0 and I is the unit matrix.

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

Temporal dependence

is given within its simulation algorithm:

• 1. Choose a fixed θ(0)
• 2. Simulate parameter vectors θ(k), k = 1, 2, . . . , T.
• 3. Simulate a realization x0 (using M-H algorithm described above).
• 4. Simulate realizations xk, k = 1, 2 . . . , T (M-H alg.) with the pro-

posal distribution Propk of newly added disc at time k given by Propk = (1 − β) · Prop(RP) + β · Prop(emp)

k−1 ,

β ∈ (0, 1), where Prop(RP) is a distribution of the reference process, Prop(emp)

k−1

is the empirical distribution obtained from the configuration xk−1 and β is a chosen constant.

Remark: Idea is that β describs the power of time dependence so that (β × 100)% of the added discs are taken from the previous configuration and the remaining discs are simulated randomly, so the dependence is stronger when β is bigger.

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

Example

A realization of the (A, L, Ncc, Nh)-interaction process in S = [0, 10] × [0, 10] with Q the uniform distribution on the interval [0.2, 0.7], θ(0) = (0.5, −0.25, −0.5, 0.5), a = (−0.1, 0.05, 0.1, −0.1), σ2 = 0.001 and β = 0.5 in times k = 0, 1, 2 (upper row) and k = 3, ..., 10 (lower row).

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

Outline

• 1. Quermass-interaction process and its extension
• 2. Simulation
• 3. Spatio-temporal Quermass-interaction process
• 4. Maximum likelihood method using MCMC
• 5. Particle filtering
• 6. Particle MCMC

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

Maximum likelihood method using MCMC simulations (MCMC MLE)

• Denote fθ(k)(x) = hθ(k)(x)/cθ(k) (i.e. hθ(k)(x) = exp{θ(k) · T(Ux)} is

the unnormalized density).

• For an observation x, the log likelihood function is given by

l(θ(k)) = log hθ(k)(x) − log cθ(k) = θ(k) · T(Ux) − log cθ(k). Problem: cθ(k) has no explicit expression. Solution: We maximize the likelihood ratio fθ(k)/fθ(k)

0 for a fixed vector

θ(k) instead.

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

MCMC MLE

• For the fixed θ(k)

0 , the log likelihood ratio

l(θ(k)) − l(θ(k)

0 ) = log(hθ(k)(x)/hθ(k)

0 (x)) − log(cθ(k)/cθ(k) 0 )

can be approximated by l(θ(k))−l(θ(k)

0 ) = log(hθ(x)/hθ(k)

0 (x))−log 1

n

n

• i=1

hθ(k)(zj)/hθ(k)

0 (zj),

where zj are realizations from fθ(k)

0 (x) obtained by MCMC simula-

tions.

• This is applied to the observations in each time k.

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

Outline

• 1. Quermass-interaction process and its extension
• 2. Simulation
• 3. Spatio-temporal Quermass-interaction process
• 4. Maximum likelihood method using MCMC
• 5. Particle filtering
• 6. Particle MCMC

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

Particle filter estimate (PFE)

• 1. In time k = 0, sample particles θ(0,i), i = 1, . . . , m, independently

from a proposal density p(θ(0)).

• 2. For times k = 1, ..., T

(a) for i = 1, . . . , m, sample ˜ θ(k,i) from q(θ(k)|θ(k−1,i)) and denote ˜ θ(0:k,i) = (θ(0:k−1,i), ˜ θ(k,i)), (b) for i = 1, . . . , m, set wi

k = f˜ θ(k,i)(xk) and normalize them,

(c) for i = 1, . . . , m, sample with replacement θ(0:k,i) from ˜ θ(0:k,i) with normalized weights from (b).

• 3. Filtered estimate is ˆ

θ(0:T) = 1

m

m

i=1 θ(0:T,i).

Remark: Here, denoting ˆ θmle the MLE estimate of θ, we set p(θ(0)) ∼ U(0, 2ˆ θ(0)

mle) for ˆ

θ(0)

mle positive

• r p(θ(0)) ∼ U(2ˆ

θ(0)

mle, 0) for ˆ

θ(0)

mle negative, and q(θ(k)|θ(k−1,i)) ∼ N(ˆ

θ(k−1,i)

mle

+ ˆ a, ˆ σ2), where ˆ a, ˆ σ2 are obtained by standard linear regression methods from MLE estimates.

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

Comparing MCMC MLE and PFE

2 4 6 8 10 6 4 2 2 4 6 time theta1 2 4 6 8 10 6 4 2 2 4 6 time theta1 2 4 6 8 10 6 4 2 2 4 6 time theta2 2 4 6 8 10 6 4 2 2 4 6 time theta2

Comparing the real parameters (solid line) with envelopes (dotted lines) and averages (dashed lines)

• f estimates obtained from 39 realizations of the process by MCMC MLE (left) and PFE (right).

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

Comparing MCMC MLE and PFE

2 4 6 8 10 6 4 2 2 4 6 time theta3 2 4 6 8 10 6 4 2 2 4 6 time theta3 2 4 6 8 10 6 4 2 2 4 6 time theta4 2 4 6 8 10 6 4 2 2 4 6 time theta4

Comparing the real parameters (solid line) with envelopes (dotted lines) and averages (dashed lines)

• f estimates obtained from 39 realizations of the process by MCMC MLE (left) and PFE (right).

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

Comparing MCMC MLE and PFE

• MCMC MLE seems to be better in average in later times
• PFE has smaller variance

↓ Possible reason: small number of components and number of holes in earlier and later times, respectively. ↓ Next steps:

• Application to Quermass-interaction process
• Consider longer time sequence
• Try also particle MCMC

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

Outline

• 1. Quermass-interaction process and its extension
• 2. Simulation
• 3. Spatio-temporal Quermass-interaction process
• 4. Maximum likelihood method using MCMC
• 5. Particle filtering
• 6. Particle MCMC

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

Particle MCMC

• Idea: Particle filter running in more iterations.
• Algorithm:
• 1. Iteration 0:

(a) Set (θ(0)(0), a(0), σ2(0)) arbitrarily (e.g. ML estimates). (b) Using steps 2 and 3 from PFE algorithm, obtain θ(0:T)(0).

• 2. Iteration t + 1:

(a) Given (θ(0)(t), a(t), σ2(t)), propose (θ(0)∗, a∗, σ2∗) (e.g. ran- dom walk). (b) Using steps 2 and 3 from PFE algorithm, obtain θ(0:T)∗. (c) Accept this proposal (i.e. set (θ(0)(t+1), a(t+1), σ2(t+1)) = (θ(0)∗, a∗, σ2∗)) with probability min(1; MH), where MH is Metropolis-Hastings ratio.

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

Metropolis-Hastings ratio in particle MCMC

MH = p1(θ(0)∗, a∗, σ2∗) p1(θ(0)(t − 1), a(t − 1), σ2(t − 1)) · p2((θ(0)(t), a(t), σ2(t))|(θ(0)∗, a∗, σ2∗)) p2((θ(0)∗, a∗, σ2∗)|(θ(0)(t), a(t), σ2(t))) ·

• k ¯

w∗

k

• k ¯

wk , where ¯ wk = 1

m

m

i=1 fθ(k,i)(xk) and analogously ¯

w∗

k = 1 m

m

i=1 fθ(k,i)∗(xk)

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

Application

A realization of the Quermass-interaction process in S = [0, 10] × [0, 10] with Q the uniform distribution on the interval [0.2, 0.7], θ(0) = (1, −0.5, −1), a = (−0.1, 0.05, 0.1), σ2 = 0.001 and β = 0.5 in times k = 0, 5, 10 (upper row) and k = 15, 20, 25 (lower row).

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

Comparing all three methods

5 10 15 20 25 2 1 1 2 time theta1

Comparing the real parameters (black line) with envelopes obtained from 19 realizations of the process by MCMC MLE (red lines) and PFE (blue lines) and PMCMC (green lines).

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

Comparing all three methods

5 10 15 20 25 1.0 0.5 0.0 0.5 1.0 1.5 2.0 time theta2

Comparing the real parameters (black line) with envelopes obtained from 19 realizations of the process by MCMC MLE (red lines) and PFE (blue lines) and PMCMC (green lines).

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

Comparing all three methods

5 10 15 20 25 1.5 1.0 0.5 0.0 0.5 1.0 1.5 time theta3

Comparing the real parameters (black line) with envelopes obtained from 19 realizations of the process by MCMC MLE (red lines) and PFE (blue lines) and PMCMC (green lines).

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

References

• 1. Andrieu C., Doucet A. and Holenstein R. (2010): Particle Markov chain Monte Carlo methods.

Journal of the Royal Statistical Society 72(3), 269-342.

• 2. Doucet A., de Freitas N., Gordon N. (2001): Sequential Monte Carlo methods in practice.

Springer, New York.

• 3. Moeller J., Helisov´

a K. (2008): Power diagrams and interaction processes for unions of discs. Advances in Applied Probability 40(2), 321-347.

• 4. Moeller J., Helisov´

a K. (2010): Likelihood inference for unions of interacting discs. Scandina- vian Journal of Statistics 37(3), 365-381.

• 5. Zikmundov´

a M., Staˇ nkov´ a Helisov´ a K., Beneˇ s V. (2012): Spatio-temporal model for a random set given by a union of interacting discs. Methodology and Computing in Applied Probability. DOI: 10.1007/s11009-012-9287-6.

• 6. Zikmundov´

a M., Staˇ nkov´ a Helisov´ a K., Beneˇ s V. (2012): On the use of particle Markov chain Monte Carlo in stochastic geometry. In preparation.

May 9, 2012 Kateˇ rina Staˇ nkov´ a Helisov´ a Estimating parameters in spatio- temporal Quermass- interaction process

◭◭ ◮◮ ◭ ◮ Back Close

Thank you for your attention!