###### tags: `Computational Statistics` `Spatial Statistics` `Variogram` `Kriging`

打這篇筆記主要是上了台大土木系電輔組汪立本老師的課，很喜歡老師的教學方式，學到超多空間統計相關的知識跟實際應用，也因為該門課我的研究很順利的產出（以下很多圖片是來自老師的講義、內容就是我的碩士論文有用到的概念，很認真的打了這個筆記，希望台大學弟妹都能幸運選到該門課。）

# What is Spatial Statistics?

在現今社會上，很多data都是來自於偵測地球表面上的狀態，像是雨量、溫度或是礦石等等。在這些資料分析上，會發現到data具備空間上的特徵，並非是random，我們可以透過地理或區域上探討資料之間的空間幾何資訊，因此運用空間統計分析，我們可以進一步瞭解到資料的特徵。

:::info
舉例來說，當有個地方採集到礦石，有很高的機率其周圍地區也會擁有礦石的痕跡(Clustering)，因此透過空間統計，可以量化出資料之間的空間特徵。
:::

進一步探討空間統計(Spatial Statistics)，要先認識變異圖**Variogram**。

---
## 變異圖 (Variogram)
* 能夠透過Variogram抓到空間中資料之間的特徵，例如彼此之間的關聯性。

### Stationary and Intrinsic Hypotheses
***各變數意義***
>$X$ - 二維空間上的某一位置，包含參數$x,y$。
$Z(X)$ - 對應在空間上的value
$h$ - 兩個data在空間上的距離
$C(h)$ - Covariance with h，只跟距離$h$有關，方向無關
$m$ - 期望值
$\gamma(h)$ - Variogram

:::success
- **Stationary Random Fields (較狹隘的)**:

1. $E(Z(X)) = m(x) = m$ 表示每個變數之期望值都相同。
2. $E(Z(X)Z(X+h))-m^2 = C(h)$

- **Intrinsic Random Fields (較廣泛的)**:
在這個假設下，我們無法得知數值之期望值$m$，但能夠知道在某一距離下$h$，兩個變數之期望值會相等。

1. $E(Z(X)-Z(X+h)) = 0$
2. $Var(Z(X+h)-Z(X)) = 2\gamma(h)$

<font color="orange">**以上兩個定義會作為後續延伸應用到Kriging的推導公式**</font>

:::

### Definition of Variogram
根據第一節的Intrinsic random fields，我們可以得到下列公式。

$$
\begin{align}
\gamma(h) = 0.5Var(Z(X+h)-Z(X))\\
= 0.5E[(Z(X+h)-Z(X))^2]
\end{align}
$$

==從下圖中，可以看到Variogram plot裡面有幾個專有名詞，他們所隱含的物理意義是，當有個事物距離你所在地在range以內，你們彼此之間是有相關性的(correlated)，若在距離為0的情況，在真實世界上會有個Nugget effect，就是本身資料含有一些量測誤差而導致的，而超過一個距離之後(sill)，彼此之間不相關的(uncorrelated)，在圖上形成一條近似水平線的樣態。==

[圖片來源](https://vsp.pnnl.gov/help/vsample/Kriging_Variogram.htm)
- 更多參考資料
[Semi-Variogram: Nugget, Range and Sill](https://gisgeography.com/semi-variogram-nugget-range-sill/#semivariogram-nugget-range-and-sill)

---

### Spatial Variogram $\gamma(h)$ & Covariance Functions $C(h)$

$$
\begin{split}
2\gamma(h) & = E[(Z(X+h)-Z(X))^2] \\
& = E[(Z(X+h)-m)^2+(Z(X)-m)^2-2(Z(X+h)-m)(Z(X)-m)] \\
& = 2C(0)-2C(h)
\end{split}
$$

$2C(0)=Variance$

:::info
$C(h)$代表兩者之間的相似性，當距離愈遠，相似性隨之降低；$\gamma(h)$則是看兩者之間的差異性多大。&rarr;**能知道$C(h)$和$\gamma(h)$在物理意義上為相反的關係**

[圖片來源](https://www.researchgate.net/profile/Li-Pen-Wang)
:::

---
### Common Variogram models
- **Nugget Effect model**
該模型顯示資料間沒有任何空間相關性的存在，只有本身資料的不確定性(Nugget effect)。

$$\begin{split}
\gamma(h)=0、h=0、|h|>0
\end{split}
$$


- **Spherical model**
該模型表示隨著距離從0開始，先呈現線性，當距離愈遠會到達一個臨界值(sill)，往後能夠看到線呈現水平樣，表示在到達sill以後，資料之間就無空間上的相關性。公式中，range = effective range $(\alpha = \gamma)$

$$\begin{split}
\gamma(\mathrm{h})=\left\{\begin{array}{c}C\left(\frac{3}{2} \frac{|h|}{a}-\frac{1}{2}\left(\frac{|h|^3}{a^3}\right)\right) \quad|h|<a \\ C \quad|h| \geq a\end{array}\right.
\end{split}
$$


- **Exponential model**
一開始跟spherical model很像，呈現線性關係，但是能夠看到沒有一個很明顯的sill能夠讓線呈現水平。此時該模型就會納入一個terminology叫*effective range*。

$$\begin{split}
\gamma(h)=C(1-e^{-\frac{h}{a}})=0.95C
\\
(1-e^{-\frac{h}{a}})=0.95
\end{split}
$$

$$\begin{split}
\alpha = \frac{r}{3}、\alpha=range、r=effective range
\end{split}
$$


- **Gaussian model**
該模型與exponential相比，更加平滑，同時也有effective range的概念。
$$\begin{split}
\gamma(h)=C(1-e^{-\frac{h^2}{a^2}})=0.95C
\\
(1-e^{-\frac{h^2}{a^2}})=0.95
\end{split}
$$

$$\begin{split}
\alpha = \frac{r}{\sqrt 3}、\alpha=range、r=effectiverange
\end{split}
$$


- **Stable model**
該模型內含變數$S$，當$S=1$則變成exponential、$S=2$則為Gaussian。
$$\begin{split}
\gamma(h)=C(1-e^{-\frac{h^S}{a^S}})=0.95C
\\
(1-e^{-\frac{h^S}{a^S}})=0.95
\end{split}
$$

$$\begin{split}
\alpha = \frac{r}{\sqrt[S]{3}}、\alpha=range、r=effective range
\end{split}
$$

#### Comparison models effect
從這邊可以看到Gaussian得出的模型相對於其他兩個模型平滑，可以從上述的Variogram model驗證該結果，因為gaussian的曲線非呈線性關係。

---
### How to calculate the Variograms?

根據擁有的data數量計算每一個點至其他點之距離，並將二者之間所獲得的值相減，得出殘差之後，計算平方和，最終得出variograms，之所以要除2是因為兩個點之間的距離會被重複計算兩次。

$$\begin{split}
\gamma^*(h)=\frac{1}{2N(h)}\sum_{i=1}^{N(h)}[Z(x_i+h)-Z(x_i)]^2
\end{split}$$


---
### Fitting a Variogram model

>若資料本身可以用first or second moment，以gaussian來表示，則使用kriging會較佳，否則資料的好壞大部分會仰賴於站台數量(原始數量的多寡)

根據前面章節所述，我們必須要了解 **==nugget effect, range, sill以及slope of the model of the origin==**，才能找出最適合的model。

另外在前一章節計算variogram中可以知道，當距離一地區愈遠，其variogram愈大，代表standard deviation(uncertainty)愈高，因此在fitting的過程，需要引入**Weighted Least Square Fitting**。


## 克利金法 (Kriging)

>使用kriging的意義為何？在真實世界中，若要預測一地區的雨量，我們主要是依靠當地的站台，直接量測到的雨量（known value)為主，但在每個地區不可能都有新建相關的站台量測數值，
#### Different types of kriging
- Simple Kriging
If you know the mean of data, you should use simple Kriging, which make use of the covariance from data.

- Ordinary Kriging (OK)
Supposed the data has intrinsic random field, we don't know the mean of data, but we know they are the same(stationary). We use variogram to perform kriging, and in calculation of kriging we know the matrix is created by the variogram in data.

- Universal Kriging (UK)
In the data field, there has underlying trend, such as time series, the mean of data is not stable, if you don't consider the trend, you would get a bad result for your prediction.

- Kriging External Drift (KED)
For example, we all know the temperature and the height of the terrain has a high correlated-relationship, so we can say height is a kind of (additional information) which seen as external drift. In the above condition use the KED for kriging will get a better estimation.

:::danger
在這邊要注意，若data本身契合gaussian model，但是在實際空間上兩筆距離相近的資料差異很大的話，會導致權重有可能會出現極大或極小值，甚至可能出現負值的情況。致使Kriging出來的成果會是很奇怪的情況。
:::