# 統計與大氣科學 HW1 相關補充資料

## 季節性週期 (Seasonal cycle)

### 季節性週期的意義

- 季節性週期是指將長期氣候資料以年為間隔切分，並把每一年中相同時間點的資料做平均。因此，若原始資料為月資料，則計算出來的季節性週期會是一個長度 12 的數列（每年的一月平均、二月平均......一直至十二月平均）。若原始資料為日資料，則季節性週期的長度為 365。
- 除了季節性週期 (Seasonal cycle) 外，其他常用的還包括日夜週期 (diurnal cycle)。這些時間長度都是某個外在的地球系統驅動力的週期，例如季節性週期就是太陽軌跡的周年運動，日夜週期就是太陽軌跡的周日運動，這些運動會影響太陽給予地球系統的太陽輻射，影響大氣所接收到的能量。因此，使用這些時間來計算週期循環則可以得知長期下來的氣候系統對於這些外在驅動力的反應。

### 去除季節性週期的計算方式

- 標準做法為先計算出季節性週期，再將原始資料扣除該季節性週期。
- 在時間序列分析中，有其他扣除季節性的方式，例如使用虛擬變數 (dummy variable) 或季節性差分 (seasonal differencing)。但這些方法是為了直接解決時間序列資料非定態 (non-stationary) 的問題，加快自相關函數的截尾，以利使用各種時間序列模型（如 ARMA、ARCH 等）。在氣候資料當中，計算季節性週期是比較主流的方式。因為可以先得到一個氣候長期平均的季節性週期，以原始資料扣除則代表每一筆資料相對於季節性週期的距平值（Anomaly，亦即與平均值的差距）。因此，處理後所得到的資料具有氣候上的意義。本次作業有些同學使用季節性差分來去除季節性，雖然處理後的結果與季節性週期的方法相似，且去除季節性後的資料也服從常態分布，但其數值就沒有距平值的含意。因此還是建議各位同學使用計算季節性週期的方式來去除季節性。
- 有些同學特別使用三個月作為切分進行資料平均。雖然這樣也可以達到相同的效果，可以得到具有季節特徵的震盪，但是季節性週期的長度會只剩 4 個季節，能保留下來的資訊就會變少。
- 有些同學把原始資料直接用第二張圖中的四個資料集中的區域來分組，然後將四組分別扣掉自己的平均。這樣雖然可以把季節的效應去除，但這樣的做法比較沒有說服力。因為用人工的方式去判斷四個區域，標準不夠明確。因此在實務上，扣除季節性週期的做法會比較標準，建議同學採用這樣的程序。

### 去除季節性週期的語法

在 Python 中，可以使用 reshape 函式來重新排列資料，計算季節性週期的語法可以更簡潔。

import numpy as np

# Calculate seasonal cycle of TS and tile the array of seasonal cycle
TS_SeasonalCycle = np.tile(np.nanmean(TS.reshape((-1, 12)), axis=0), TS.shape[0]//12)

# Calculate seasonal anomaly
TS_Anomaly = TS - TS_SeasonalCycle

TS_Anomaly 即為去除季節性週期後的 TS 資料。

- 上述範例中，reshape 會先將時間維度拆成 "年分×12" 兩個維度，nanmean 再對年分的維度取平均，即可算出十二個月分別的長期平均。但經過這樣處理後的時間維度只剩 12，因此需要將時間維度重新擴大到原始資料的尺寸。使用 tile 函式可以把資料複製並排列，其中第二個引數的 TS.shape[0]//12 即代表原始資料的年數，tile 函式會把 12 筆資料重新複製 TS.shape[0]//12 次，並沿著資料方向堆疊，使長度與原始資料的長度相同。關於 tile 函式的用法有下列補充：
- 如果需要對多維資料做 tile，則第二個引述需要傳入 tuple 來指定每一個維度分別要堆疊的個數。例如 np.tile(A, (2, 3)) 就表示把 A 陣列沿著第零維度堆疊兩次、第一維度堆疊三次，而其中 A 陣列必須是二維，否則編譯會報錯。
- tile 函式跟 repeat 函式的差異在於堆疊的方式。[1, 2, 3] 陣列如果 tile 兩次會變成 [1, 2, 3, 1, 2, 3]，如果 repeat 兩次則會變成 [1, 1, 2, 2, 3, 3]。在分析資料時，tile 可以用於把短週期的資料擴展至長時間尺度（例如本題中的季節性週期），而 repeat 可以把低頻資料重複變成高頻資料（例如把日資料 repeat 24 次變成小時資料）。
- 其中，reshape 後方所接的 (-1, 12) 代表重新排列後的陣列尺寸。-1 代表第零維度的長度未定，這時候函式會根據第一維度的大小與資料的總長度自動判斷第零維度的長度。在處理大量的大氣資料時，有時候可能還尚未確定資料的總長度，這時候可以善用這種方式來做 reshape。
- 另外，reshape 指定的維度順序也會影響資料排列。在 reshape 的預設下，order 引數的值為 'C'，重排後的資料會優先從靠後的維度開始排列 (the last axis index changing fastest)。因此每個月的資料會先優先排列在靠後的維度，等到十二筆資料排滿後，才會開始排列第二年的資料。此時，第零維度代表年分，第一維度則代表該年的十二個月的資料。例如 TS.reshape((-1, 12))[4, 6] 則可以取出 TS 資料中第五年的第七個月的資料。此時若要計算季節性週期，對第零維度取平均即可。
- 大部分的同學都使用 numpy.average 來做資料平均。average 函式的特性是可以計算權重，但是無法處理無效值。因此建議各位同學在未來分析氣候資料時，可以使用 numpy.nanmean 來計算平均（除了不能計算權重外，其他的功能皆與 average 類似）。如果資料中具有無效值 nan，nanmean 函式會略過所有 nan 並使用有效的值來計算平均，而 average 函式則會直接輸出無效值 nan。

reshape 函式用法參照：https://numpy.org/doc/stable/reference/generated/numpy.reshape.html

另外，使用迴圈來計算逐月的平均並直接扣除也是簡潔的做法：

import numpy as np

TS_Anomaly = TS

for i_Month in range(12):
    
    TS_Anomaly[i_Month::12] = TS[i_Month::12] - np.nanmean(TS[i_Month::12])

- 由於迴圈中每一次扣除都只會動到 i_Month 的那個月分，因此每一次迴圈所操作的資料不會彼此影響。迴圈迭代 12 次後，即可完成每一筆資料的季節性循環扣除。不過，Python 中的 for 迴圈仍較占用記憶體，若迭代次數較多仍可能占用運算資源。
- Numpy array 中的 index 若為 [a : b : c]，表示選取的陣列為從 a 開始、結束為 b、間距為 c，且最後一個元素不包含 b 的資料（記得 Python 第一個元素的 index 為零）。

## Python 中列表 (List) 與陣列 (Array) 的使用

在大量計算的需求下，Python 使用陣列 (Array) 來計算的效率會快得多。陣列會直接在主機上請求一整塊記憶體的空間，因此存取資料時不需要逐個元素去尋找記憶體。由於陣列的尺寸、資料型態在宣告時就已經確定，無法隨意更改，因此存取速度較快。相反的，列表內元素的彈性比較大，可以在同一個列表裡面放入不同資料型態的元素，而且可以增加、減少陣列的長度。但由於列表存取資料時的位置在記憶體當中是打散的，因此需要另外使用 Pointer 來指向每個元素所對應的記憶體位置，處理的速度會慢得多。另外，陣列也能進行多維度的資料處理、重新排列、計算等工作，列表則較難進行。
在 Python 中，Pandas 套件可以視為列表跟陣列的折衷。Pandas 的 DataFrame 可以在不同的欄存取不同資料型態的資料，而且資料處理的函式也很多。在未來分析時間序列資料、資料分組統計等工作時，建議各位同學可以使用 Pandas 套件。

Pandas 套件的說明書參照：https://pandas.pydata.org/docs/

## 把 txt 檔讀成 Numpy array

用 loadtxt 函式：

import numpy as np

TS = np.loadtxt('TS.txt')

由於 TS.txt 裡面的資料本身就是一維，loadtxt 讀進來的資料也會是一維 Numpy array。

## 在 Matplotlib 套件中使用 Latex 語法

分析科學資料的時候，有時候會需要在圖片的文字中輸入特殊符號，這時候可以善用 Latex 語法。在 Matplotlib 中，如果使用 Latex 語法，則可以在字串內使用 $ 符號將 Latex 語法包住，並在字串前方加上 r。

ax.plot(x, y, label=r'$\lambda$')

如上述程式碼，即可設定 Plot 物件的 label 屬性，在顯示文字時即可顯示為 lambda 符號。字串前的 r 是為了避免字串內的任何特殊符號被 Python 轉義成特殊字元。

---
這次作業有一些同學放了一些迷因圖，這樣很棒。以後再辦一個統計學迷因比賽 : )