Interpolation for huge spatial datasets: theory, implementations and - PowerPoint PPT Presentation

Interpolation for huge spatial datasets: theory, implementations and ideas @ReinhardFurrer, I-Math/ICS, UZH NZZ.ch Uppsala, 2019/01/30

Joint work ◮ Emilio Porcu ◮ Florian Gerber ◮ Francois Bachoc and with contributions of several others 2

Outlook ◮ Motivation ◮ Covariance approximations ◮ Gaussian equivalent measures ◮ Prediction and estimation ◮ Implementation and illustration � Questions and ≪ Fast-Forward ≫ appreciated! Slides at http://www.math.uzh.ch/furrer/slides/ 3

Spatial statistics: prediction Observations: y ( s 1 ) , . . . , y ( s n ) y (s ) y (s ) 2 n y (s ) (s ) y s i 1 2 s n s 1 s i Source: wikipedia.org First law of geography (Waldo Tobler): Everything is related to everything else, but near things are more related than distant things. 4

Spatial statistics: prediction Predict the quantity of interest at an arbitrary location. noise signal s 2 Why? s ◮ Fill-in missing data n s 1 s i ◮ Force data onto a regular grid ◮ Smooth out the measurement error How? ◮ By eye ◮ Linear interpolation ◮ The correct way . . . 4

Visual example Source: Gerber et al (2018), TGRS 7

Visual example Source: Gerber et al (2018), TGRS 8

Spatial statistics: prediction Predict Z ( s 0 ) given y ( s 1 ) , . . . , y ( s n ) s Y ( s ) = x T ( s ) β + µ ( s ) + Z ( s ) + ε ( s ) 2 ? s n s 1 s i Minimize mean squared prediction error (over all linear unbiased predictors) s 0 � Best Linear Unbiased Predictor: � � � � − 1 obs BLUP = Cov Z ( s predict ) , Y ( s obs ) Var Y ( s obs ) Z ( s 0 ) = c T Σ − 1 y � (one spatial process, no trend, known covariance structure; otherwise almost the same) 9

Spatial statistics: prediction Predict Z ( s 0 ) given y ( s 1 ) , . . . , y ( s n ) s Y ( s ) = x T ( s ) β + µ ( s ) + Z ( s ) + ε ( s ) 2 s n s 1 s i Covariance C ( dist ( s 1 , s 2 )) ● C ( dist ( s 1 , s n )) ● 0.0 0.2 0.4 0.6 0.8 Distance, lag h � � Covariance matrix Σ contains elements C dist( s i , s j ) . 9

Issues of basic, classical kriging Cov(pred , obs) · Var(obs) − 1 · obs = c Σ − 1 y ◮ “Simple” spatial interpolation . . . . . . on paper or in class! ◮ BUT: 1. Complex mean structure 2. Unknown covariance function, unknown parameters 3. Large spatial fields 4. Non-stationary covariances 5. Multivariate/space-time/data on the sphere Many R packages . . . . . . many black boxes . . . Heaton et al. JABES 10

Outlook ◮ Motivation ◮ Covariance approximations ◮ Gaussian equivalent measures ◮ Prediction and estimation ◮ Implementation and illustration 12

Spatial modeling Lattice model (GMRF): Geostatistical model (GRF): s 2 s n s 1 s i � E( Z i | z − i ) = β z j Covariance j neighbor of i C ( dist ( s 1 , s 2 )) Var( Z i | z − i ) = τ 2 ● C ( dist ( s 1 , s n )) ● 0.0 0.2 0.4 0.6 0.8 Gaussianity and Distance, lag h regularity conditions: Σ = τ 2 ( I − B ) − 1 Covariance matrix: Σ 13

Spatial modeling Geostatistical model (GRF): Lattice model (GMRF): Σ − 1 Σ Σ app Σ 14

Sparseness Using sparse covariance functions for greater computational efficiency. Sparseness is guaranteed when the covariance function has a compact support ◮ a compact support is (artificially) imposed � tapering ◮ Matern ν = 1.5 Wendland Wendland Wendland Matern * Wendland 0 10 20 30 40 0 10 20 30 40 Distance, lag h Distance, lag h 16

Tapering: asymptotic optimality For an isotropic and stationary process with covariance C 0 ( h ), find a taper C ρ ( h ), such that kriging with C 0 ( h ) C ρ ( h ) is asymptotically optimal. MSPE( s 0 , C 0 C ρ ) = MSPE taper → 1 MSPE( s 0 , C 0 ) MSPE true MSPE( s 0 , C 0 ) = MSPE naive ̺ ( s 0 , C 0 C ρ ) ̺ ( s 0 , C ) = C (0) − c T Σ − 1 c → 1 MSPE true 16

Tapering: asymptotic optimality For an isotropic and stationary process with covariance C 0 ( h ), C 1 ( h ) C 1 ( h ) find a taper C ρ ( h ), such that kriging with C 0 ( h ) C ρ ( h ) is asymptotically optimal. MSPE( s 0 , C 0 C ρ ) = MSPE taper → 1 MSPE( s 0 , C 0 ) MSPE true MSPE( s 0 , C 0 ) = MSPE naive ̺ ( s 0 , C 0 C ρ ) ̺ ( s 0 , C ) = C (0) − c T Σ − 1 c → 1 MSPE true Proofs based on infill asymptotics and “misspecified” covariances Conditions on the tail behaviour of the spectrum of the (tapered) covariance Furrer, Genton, Nychka (2006) JCGS Stein (2013) JCGS Bevilacqua et al (2019) AoS 16

Covariance functions Mat´ ern: M ν ( r ) = 2 1 − ν Γ( ν ) r ν K ν ( r ) ν > 0 M ν,α,σ 2 ( r ) = σ 2 M ν ( r/α ) σ 2 > 0 α > 0 Generalized Wendland: ∞ � 1 u ( u 2 − r 2 ) κ − 1 (1 − u ) µ W µ,κ ( r ) = + d u B (2 κ, µ + 1) r µ ≥ ( d + 1) / 2 + κ κ > 0 W µ,κ,β,σ 2 ( r ) = σ 2 W µ ( r/β ) σ 2 > 0 β > 0 17

Outlook ◮ Motivation ◮ Covariance approximations ◮ Gaussian equivalent measures ◮ Prediction and estimation ◮ Implementation and illustration 18

Equivalence of Gaussian measures P 0 , P 1 two probability measures on a measureable space { Ω , S} . Definition . P 0 , P 1 are equivalent if P 0 ( A ) = 1, for any A ∈ S implies P 1 ( A ) = 1 and vice versa. ◮ Restrict A to the σ -algebra generated by { Z ( s ) , s ∈ D } Here: ‘ ≡ ’ equivalent on the paths of { Z ( s ), s ∈ D } . ◮ For Gaussian measures: equivalence or orthogonality first two moments discriminate: P ( µ, C ) ◮ Here: µ = 0 and D ⊂ R d , d = 1 , 2 , 3 19

Equivalence of Gaussian measures Covariance functions C 0 , C 1 , associated spectral densities C F 0 , C F 1 . D ⊂ R d any bounded infinite set, s 0 ∈ D . Theorem. (Stein 2004) 0 ( z ) z a ≤ b u < ∞ , 0 < b ℓ ≤ C F If ∃ a > 0, z → ∞ � � 2 ∞ � C F 1 ( z ) − C F 0 ( z ) z d − 1 and if d z < ∞ , C F 0 ( z ) 0 <b then P ( C 0 ) ≡ P ( C 1 ). 20

Fourier transforms 1 M F ν ( z ) = cte · z ≥ 0 . (1 + z 2 ) ν + d/ 2 , � � 2; − z 2 λ ; λ + µ 2 , λ + µ 2 + 1 W F µ,κ ( z ) = cte · 1 F 2 , λ = ( d + 1) / 2 + κ 4 ≍ z − 2 λ z → ∞ , for 21

Equivalence of Gaussian measures Theorem. (Zhang 2004) 1 ), D ⊂ R d , d = 1 , 2 , 3. Let P ( M ν,α 0 ,σ 2 0 ), P ( M ν,α 1 ,σ 2 � � σ 2 α 2 ν = σ 2 α 2 ν P ( M ν,α 0 ,σ 2 0 ) ≡ P ( M ν,α 1 ,σ 2 ⇐ ⇒ 1 ) 0 0 1 1 ◮ d ≤ 3: not all parameters consistently estimable microergodic parameter σ 2 � α 2 ν ◮ d ≥ 5 all fine (Anderes 2010) ◮ d = 4 open . . . 22

Equivalence of Gaussian measures Theorem. (BFFP 2019) Let P ( W µ,κ,β 0 ,σ 2 0 ), P ( W µ,κ,β 1 ,σ 2 1 ). For a given κ ≥ 0, let µ > d + 1 / 2 + κ : � � β 2 κ +1 β 2 κ +1 σ 2 = σ 2 P ( W µ,κ,β 0 ,σ 2 0 ) ≡ P ( W µ,κ,β 1 ,σ 2 ⇐ ⇒ 1 ) 0 0 1 1 Theorem. (BFFP 2019) 1 ), D ⊂ R d , d = 1 , 2 , 3, κ > 0. Let P ( M ν,α,σ 2 0 ) and P ( W µ,κ,β,σ 2 If ν = κ + 1 / 2, µ > d + 1 / 2 + κ : � � α 2 ν = cte ν,κ,µ · σ 2 σ 2 β 2 κ +1 P ( M ν,α,σ 2 0 ) ≡ P ( W µ,κ,β,σ 2 ⇐ ⇒ 1 ) 0 1 23

Outlook ◮ Motivation ◮ Covariance approximations ◮ Gaussian equivalent measures ◮ Prediction and estimation ◮ Implementation and illustration 24

Prediction: asymptotic optimality Prediction of Z ( s 0 ) using “BLUP” under misspecification � Z n ( µ, κ, β ). Theorem. (BFFP 2019) P ( W µ,κ,β,σ 2 ): if µ > ( d + 1) / 2 + κ and for any fixed β 1 > 0 � � � Var µ,κ,β,σ 2 Z n ( µ, κ, β 1 ) − Z ( s 0 ) � − → 1 � � Var µ,κ,β,σ 2 Z n ( µ, κ, β ) − Z ( s 0 ) P ( M ν,α,σ 2 ): if ν = κ + 1 / 2 and for any fixed β 1 > 0 � � � Var ν,α,σ 2 Z n ( µ, κ, β 1 ) − Z ( s 0 ) � � � − → 1 Z n ( ν, α ) − Z ( s 0 ) Var ν,α,σ 2 25

Prediction: presumed MSPE Prediction of Z ( s 0 ) using “BLUP” under misspecification � Z n ( µ, κ, β ). Theorem. (BFFP 2019) � � β 2 κ +1 β 2 κ +1 if σ 2 = σ 2 P ( W µ,κ,β,σ 2 ): 0 0 1 1 � � � Var µ,κ,β 1 ,σ 2 Z n ( µ, κ, β 1 ) − Z ( s 0 ) 1 � − → 1 � � Var µ,κ,β,σ 2 Z n ( µ, κ, β 1 ) − Z ( s 0 ) if ν = κ + 1 / 2, σ 2 � � α 2 ν = cte · σ 2 β 2 κ +1 P ( M ν,α,σ 2 ): 1 1 � � � Z n ( µ, κ, β 1 ) − Z ( s 0 ) Var µ,κ,β 1 ,σ 2 1 − → 1 � � � Var ν,α,σ 2 Z n ( µ, κ, β 1 ) − Z ( s 0 ) 26

Estimation: tapering Likelihood: − 1 2 log det( Σ ( θ ) ◦ T ) − 1 2 y T ( Σ ( θ ) ◦ T ) − 1 y 2 y T � � − 1 2 log det( Σ ( θ ) ◦ T ) − 1 ( Σ ( θ ) ◦ T ) − 1 ◦ T y Dimension plays crucial role Kaufman et al. (2008) JASA ◮ Two-taper approach cannot be used in practice ◮ 27

Estimation: Mat´ ern Let P ( M ν,α 0 ,σ 2 0 ), with certain constrains on the parameter space. Theorem. (Shaby Kaufmann 2013) σ 2 Let { ˆ n , ˆ α n } the ML estimator. Then as n → ∞ � � a.s. σ 2 α 2 ν → σ 2 α 2 ν 1. ˆ ˆ − n n 0 0 � 0 ) 2 � √ n (ˆ D σ 2 α 2 ν n − σ 2 0 /α 2 ν 0 , 2( σ 2 0 /α 2 ν − → N 2. n / ˆ 0 ) σ 2 Also holds for ˆ n ( α ) with misspecified α ◮ 28