Lecture 14 Covariance Functions 3/08/2018 1 More on Covariance - - PowerPoint PPT Presentation

Lecture 14

Covariance Functions

3/08/2018

1

More on Covariance Functions

2

Nugget Covariance

𝐷𝑝𝑤(𝑧𝑢𝑗, 𝑧𝑢𝑘) = 𝜏21{ℎ=0} where ℎ = |𝑢𝑗 − 𝑢𝑘|

0.00 0.25 0.50 0.75 1.00 −2 −1 1 2 5 10 15 20 5 10 15 20

h x C y draw

Draw 1 Draw 2

3

(- / Power / Square) Exponential Covariance

𝐷𝑝𝑤(𝑧𝑢𝑗, 𝑧𝑢𝑘) = 𝜏2 exp (−(ℎ 𝑚)𝑞) where ℎ = |𝑢𝑗 − 𝑢𝑘|

0.00 0.25 0.50 0.75 1.00 −4 −2 2 −3 −2 −1 1 −2 2 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00

h x x x C y y y Cov

Exp Pow Exp (p=1.5) Sq Exp

Exponential Powered Exponential (p=1.5) Square Exponential Covariance − l=12, sigma2=1

4

Matern Covariance

𝐷𝑝𝑤(𝑧𝑢𝑗, 𝑧𝑢𝑘) = 𝜏2 21−𝜉 Γ(𝜉) ( √ 2𝜉 ℎ ⋅ 𝑚)

𝜉 𝐿𝜉 (

√ 2𝜉 ℎ ⋅ 𝑚) where ℎ = |𝑢𝑗−𝑢𝑘|

0.00 0.25 0.50 0.75 1.00 −2 −1 1 2 −3 −2 −1 1 −1 1 2 2 4 6 2 4 6 2 4 6 2 4 6

h x x x C y y y

v=1/2 v=3/2 v=5/2

Matern − v=1/2 Matern − v=3/2 Matern − v=5/2 Covariance − l=2, sigma2=1

5

Matern Covariance

• 𝐿𝜉 is the modified Bessel function of the second kind.
• A Gaussian process with Matérn covariance has sample functions that

are ⌈𝜉 − 1⌉ times differentiable.

• When 𝜉 = 1/2 + 𝑞 for 𝑞 ∈ N+ then the Matern has a simplified form

(product of an exponential and a polynomial of order 𝑞).

• When 𝜉 = 1/2 the Matern is equivalent to the exponential covariance.
• As 𝜉 → ∞ the Matern converges to the square exponential

covariance.

• A Gaussian process with Matérn covariance has paths that are ⌈𝜉⌉ − 1

times differentiable.

6

Rational Quadratic Covariance

𝐷𝑝𝑤(𝑧𝑢𝑗, 𝑧𝑢𝑘) = 𝜏2 (1 + ℎ2 𝑚2 𝛽 )

−𝛽

where ℎ = |𝑢𝑗 − 𝑢𝑘|

0.00 0.25 0.50 0.75 1.00 −2 −1 1 −1 1 2 −1 1 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00

h x x x y y y

alpha=1 alpha=3 alpha=10

Rational Quadratic − alpha=1 Rational Quadratic − alpha=10 Rational Quadratic − alpha=100 Covariance − l=12, sigma2=1

7

Rational Quadratic Covariance

• is a scaled mixture of squared exponential covariance functions with

different characteristic length-scales (𝑚).

• As 𝛽 → ∞ the rational quadratic converges to the square

exponential covariance.

• Has sample functions that are infinitely differentiable for any value of

𝛽

8

Spherical Covariance

𝐷𝑝𝑤(𝑧𝑢𝑗, 𝑧𝑢𝑘) = {𝜏2 (1 − 3

2ℎ ⋅ 𝑚 + 1 2(ℎ ⋅ 𝑚)3))

if 0 < ℎ < 1/𝑚

• therwise

where ℎ = |𝑢𝑗−𝑢𝑘|

0.00 0.25 0.50 0.75 1.00 −3 −2 −1 1 −2 2 −2 −1 1 2 0.0 0.3 0.6 0.9 0.0 0.3 0.6 0.9 0.0 0.3 0.6 0.9 0.0 0.3 0.6 0.9

h x x x y y y

l=1 l=3 l=10

Spherical − l=1 Spherical − l=3 Spherical − l=10 Covariance − sigma2=1

9

Periodic Covariance

𝐷𝑝𝑤(𝑧𝑢𝑗, 𝑧𝑢𝑘) = 𝜏2 exp (−2 𝑚2 sin2 (𝜌ℎ 𝑞 )) where ℎ = |𝑢𝑗 − 𝑢𝑘|

0.00 0.25 0.50 0.75 1.00 −1 1 −2 −1 1 −2 −1 1 2 1 2 3 4 2 4 6 2 4 6 2 4 6

h x x x y y y forcats::as_factor(Cov)

p=1 p=2 p=3

Periodic − p=1 Periodic − p=2 Periodic − p=3 Covariance − l=2, sigma2=1

10

Linear Covariance

𝐷𝑝𝑤(𝑧𝑢𝑗, 𝑧𝑢𝑘) = 𝜏2

𝑐 + 𝜏2 𝑤 (𝑢𝑗 − 𝑑)(𝑢𝑘 − 𝑑)

−1.0 −0.5 0.0 0.5 1.0 0.00 0.25 0.50 0.75 1.00

x y

11

Combining Covariances

If we definite two valid covariance functions, 𝐷𝑝𝑤𝑏(𝑧𝑢𝑗, 𝑧𝑢𝑘) and

𝐷𝑝𝑤𝑐(𝑧𝑢𝑗, 𝑧𝑢𝑘) then the following are also valid covariance functions, 𝐷𝑝𝑤𝑏(𝑧𝑢𝑗, 𝑧𝑢𝑘) + 𝐷𝑝𝑤𝑐(𝑧𝑢𝑗, 𝑧𝑢𝑘) 𝐷𝑝𝑤𝑏(𝑧𝑢𝑗, 𝑧𝑢𝑘) × 𝐷𝑝𝑤𝑐(𝑧𝑢𝑗, 𝑧𝑢𝑘)

12

Linear × Linear → Quadratic

𝐷𝑝𝑤𝑏(𝑧𝑢𝑗, 𝑧𝑢𝑘) = 1 + 2 (𝑢𝑗 × 𝑢𝑘) 𝐷𝑝𝑤𝑐(𝑧𝑢𝑗, 𝑧𝑢𝑘) = 2 + 1 (𝑢𝑗 × 𝑢𝑘)

−10 −5 5 −2 −1 1 2

x y

Cov_a * Cov_b

13

Linear × Periodic

𝐷𝑝𝑤𝑏(𝑧𝑢𝑗, 𝑧𝑢𝑘) = 1 + 1 (𝑢𝑗 × 𝑢𝑘) 𝐷𝑝𝑤𝑐(𝑧𝑢𝑗, 𝑧𝑢𝑘) = exp (−2 sin2 (2𝜌 ℎ))

−4 −2 2 1 2 3

x y

Cov_a * Cov_b

14

Linear + Periodic

𝐷𝑝𝑤𝑏(𝑧𝑢𝑗, 𝑧𝑢𝑘) = 1 + 1 (𝑢𝑗 × 𝑢𝑘) 𝐷𝑝𝑤𝑐(ℎ = |𝑢𝑗 − 𝑢𝑘|) = exp (−2 sin2 (2𝜌 ℎ))

−5 −4 −3 −2 −1 1 2 3

x y draw

Draw 1 Draw 2

Cov_a + Cov_b

15

Sq Exp × Periodic → Locally Periodic

𝐷𝑝𝑤𝑏(ℎ = |𝑢𝑗 − 𝑢𝑘|) = exp(−(1/3)ℎ2) 𝐷𝑝𝑤𝑐(ℎ = |𝑢𝑗 − 𝑢𝑘|) = exp (−2 sin2 (𝜌 ℎ))

−2 −1 1 2 2 4 6

x y

Cov_a * Cov_b

16

Sq Exp (short) + Sq Exp (long)

𝐷𝑝𝑤𝑏(ℎ = |𝑢𝑗 − 𝑢𝑘|) = (1/4) exp(−4 √ 3ℎ2) 𝐷𝑝𝑤𝑐(ℎ = |𝑢𝑗 − 𝑢𝑘|) = exp(−( √ 3/2)ℎ2)

−2 −1 1 0.0 2.5 5.0 7.5 10.0

x y

Cov_a + Cov_b

17

Sq Exp (short) + Sq Exp (long) (Seen another way)

Cov_A (short) Cov_B (long) Cov_A + Cov_B 0.0 2.5 5.0 7.5 10.0 0.0 2.5 5.0 7.5 10.0 0.0 2.5 5.0 7.5 10.0 −3 −2 −1 1 2

x y 18

BDA3 example

19

BDA3

http://research.cs.aalto.fi/pml/software/gpstuff/demo_births.shtml

20

Births (one year)

• 1. Smooth long term trend

(sq exp cov)

• 2. Seven day periodic trend with

decay (periodic × sq exp cov)

• 3. Constant mean

21

Component Contributions

We can view our GP in the following ways,

𝐳 ∼ 𝒪(𝝂, 𝚻1 + 𝚻2 + 𝜏2𝐉 )

but with appropriate conditioning we can also think of 𝐳 as being the sum

• f multipe independent GPs

𝐳 = 𝝂 + 𝑥1(𝐮) + 𝑥2(𝐮) + 𝑥3(𝐮)

where

𝑥1(𝐮) ∼ 𝒪(0, 𝚻1) 𝑥2(𝐮) ∼ 𝒪(0, 𝚻2) 𝑥3(𝐮) ∼ 𝒪(0, 𝜏2𝐉 )

22

Decomposition of Covariance Components

⎡ ⎢ ⎣ 𝑧 𝑥1 𝑥2 ⎤ ⎥ ⎦ ∼ 𝒪 ⎛ ⎜ ⎝ ⎡ ⎢ ⎣ 𝝂 ⎤ ⎥ ⎦ , ⎡ ⎢ ⎣ Σ1 + Σ2 + 𝜏2𝐉 Σ1 Σ2 Σ1 Σ1 Σ2 Σ2 ⎤ ⎥ ⎦ ⎞ ⎟ ⎠ therefore

𝑥1 | 𝐳, 𝝂, 𝜾 ∼ 𝒪(𝝂𝑑𝑝𝑜𝑒, 𝚻𝑑𝑝𝑜𝑒) 𝝂𝑑𝑝𝑜𝑒 = 0 + Σ1 (Σ1 + Σ2 + 𝜏2𝐽)−1(𝐳 − 𝝂) 𝚻𝑑𝑝𝑜𝑒 = Σ1 − Σ1(Σ1 + Σ2 + 𝜏2𝐉)−1Σ1

𝑢

23

Births (multiple years)

• 1. slowly changing trend (sq exp cov)
• 2. small time scale correlating noise (sq exp cov)
• 3. 7 day periodical component capturing day of week effect (periodic × sq exp cov)
• 4. 365.25 day periodical component capturing day of year effect (periodic × sq exp cov)
• 5. component to take into account the special days and interaction with weekends (linear

cov)

• 6. independent Gaussian noise (nugget cov)
• 7. constant mean

24

Mauna Loa Exampel

25

Atmospheric CO2

330 360 390 1960 1980 2000

x y Source

NOAA Scripps (co2 in R)

26

GP Model

Based on Rasmussen 5.4.3 (we are using slightly different data and parameterization)

𝐳 ∼ 𝒪(𝝂, 𝚻1 + 𝚻2 + 𝚻3 + 𝚻4 + 𝜏2I ) {𝝂}𝑗 = ̄ 𝑧 {𝚻1}𝑗𝑘 = 𝜏2

1 exp (−(𝑚1 ⋅ 𝑒𝑗𝑘)2)

{𝚻2}𝑗𝑘 = 𝜏2

2 exp (−(𝑚2 ⋅ 𝑒𝑗𝑘)2) exp (−2 (𝑚3)2 sin2(𝜌 𝑒𝑗𝑘/𝑞))

{𝚻3}𝑗𝑘 = 𝜏2

3 (1 + (𝑚4 ⋅ 𝑒𝑗𝑘)2

𝛽 )

−𝛽

{𝚻4}𝑗𝑘 = 𝜏2

4 exp (−(𝑚5 ⋅ 𝑒𝑗𝑘)2)

27

JAGS Model

ml_model = ”model{ y ~ dmnorm(mu, inverse(Sigma)) for (i in 1:(length(y)-1)) { for (j in (i+1):length(y)) { k1[i,j] <- sigma2[1] * exp(- pow(l[1] * d[i,j],2)) k2[i,j] <- sigma2[2] * exp(- pow(l[2] * d[i,j],2) - 2 * pow(l[3] * sin(pi*d[i,j] / per), 2)) k3[i,j] <- sigma2[3] * pow(1+pow(l[4] * d[i,j],2)/alpha, -alpha) k4[i,j] <- sigma2[4] * exp(- pow(l[5] * d[i,j],2)) Sigma[i,j] <- k1[i,j] + k2[i,j] + k3[i,j] + k4[i,j] Sigma[j,i] <- Sigma[i,j] } } for (i in 1:length(y)) { Sigma[i,i] <- sigma2[1] + sigma2[2] + sigma2[3] + sigma2[4] + sigma2[5] } for(i in 1:5){ sigma2[i] ~ dt(0, 2.5, 1) T(0,) l[i] ~ dt(0, 2.5, 1) T(0,) } alpha ~ dt(0, 2.5, 1) T(0,) }” 28

Diagnostics

alpha l[1] l[2] l[3] l[4] l[5] sigma2[1] sigma2[2] sigma2[3] sigma2[4] sigma2[5] 0 250500750 1000 0 250500750 1000 0 250500750 1000 0 250500750 1000 0 250500750 1000 0.02 0.03 0.04 0.0 0.3 0.6 0.9 0.0 0.2 0.4 0.6 0.8 2 4 6 0.0 0.5 1.0 1.5 2.0 0.6 0.8 1.0 1.2 10 20 30 40 0.005 0.010 0.015 0.020 2000 4000 6000 0.02 0.04 0.06 2 4 6 8

.iteration estimate 29

Fit Components

Sigma_3 Sigma_4 Sigma_1 Sigma_2 1960 1970 1980 1990 1960 1970 1980 1990 −4 4 8 −1.5 −1.0 −0.5 0.0 0.5 1.0 −20 −10 10 20 30 −2 −1 1 2

x post_mean cov

Sigma_1 Sigma_2 Sigma_3 Sigma_4

30

Forecasting

330 360 390 1960 1980 2000

x post_mean 31

Forecasting (zoom)

360 370 380 390 400 410 2000 2005 2010 2015

x post_mean 32

Forecasting ARIMA (auto)

330 360 390 1960 1980 2000

Time co2 level

80 95

Forecasts from ARIMA(1,1,1)(1,1,2)[12]

33

Forecasting RMSE

dates RMSE (arima) RMSE (gp) Jan 1998 - Jan 2003 1.103 1.911 Jan 1998 - Jan 2008 2.506 4.575 Jan 1998 - Jan 2013 3.824 7.706 Jan 1998 - Mar 2017 5.461 11.395

34

Forecasting Components

Sigma_3 Sigma_4 Sigma_1 Sigma_2 2000 2005 2010 2015 2000 2005 2010 2015 −5 5 −1.0 −0.5 0.0 0.5 1.0 30 40 50 60 −1 1

x post_mean cov

Sigma_1 Sigma_2 Sigma_3 Sigma_4

35