線性模型練習題
1. 學習理解如何用最小二乘法的矩陣公式來獲得線性迴歸的解，並使用numpy庫來實現該算法。
在研究一個問題時，從某種理論或假定出發，獲得一個模型。根據這個模型，咱們感興趣的某個量有其理論值，同時能夠對這個量進行實際觀測，而得出其觀測值。因爲種種緣由，如模型不徹底正確以及觀測有偏差，理論值與觀測值會有差距，這差距的平方和python
\[J(\theta)=\sum (理論值-觀測值)^{2}\]git
能夠做爲理論與實測符合程度的度量。求和是針對若干次不一樣的觀測，一般理論中含有未知參數或參數向量 \(\theta\)。最小二乘法（Least Squares Method，如下簡記爲 LSE）要求選擇使 \(J(\theta)\) 達到最小的 \(\theta\) 值 \(\widehat{\theta}\)，所以，LSE 的直接意義，是做爲一種估計未知參數的方法。算法
尋找 \(\theta\) 有兩種方法，一種是梯度降低法，是搜索算法，先給 \(\theta\) 賦個初值，而後再根據使 \(J(\theta)\) 更小的原則對 \(\theta\) 進行修改，直到\(\theta\)收斂，\(J(\theta)\) 達到最小，也就是不斷嘗試；另一種是矩陣求導法，要使 \(J(\theta)\) 最小，就對 \(J(\theta)\) 求導，使導數等於 0，求得 \(\theta\)。api
上一週介紹過梯度降低法，這裏詳細推導下最小二乘法的矩陣公式。app
假定有 m 個例子的訓練集 {\((x^{(i)},y^{(i)});i=1,...,m\)}，x 是 n 維向量，是自變量，y 是實際觀測值，是因變量，假定 y 是 x 的線性函數，h(x) 是對 y 的預測值。dom
\[h(x)=\sum_{i=0}^{n}\theta_{i}x_{i}=\theta^{T}x\]機器學習
x 和 \(\theta\) 都是長度爲 n+1 的向量，其中 \(x_{0}=1\)，\(\theta_{0}\) 是截距。ide
獲得了參數 \(\theta\)，就能根據 x 計算出 y 值。那麼，怎麼根據當前的訓練數據集獲得 \(\theta\) 呢？ \(\theta\) 值確定不是隨便給，那要知足什麼要求呢？函數
毫無疑問，要想準確預測新數據，得出的\(\theta\)最起碼應該保證在訓練集上預測準確些，怎麼纔算準確些呢？計算出的價格與真實的價格應儘可能接近。計算價格和實際價格的差距以損失函數來表示。學習
\[J(\theta)=\frac{1}{2}\sum_{i=1}^{m}(h_{\theta}(x^{(i)}-y^{(i)}))^{2}\]
咱們最終選擇的\(\theta\)要使\(J(\theta)\)最小，對\(\theta\)就這個要求 。
好，如今咱們的目標就是要找到使\(J(\theta)\)最小的\(\theta\)。梯度降低法不說了，下面咱們使用矩陣方法。
首先定義：
\[X= \left[ \begin{matrix} —(x^{(1)})^{T}— \\ —(x^{(2)})^{T}— \\ . \\ . \\ . \\ —(x^{(m)})^{T}— \end{matrix} \right] \]
\[\vec{y}= \left[ \begin{matrix} (y^{(1)}) \\ (y^{(2)}) \\ . \\ . \\ . \\ (y^{(m)}) \end{matrix} \right] \]
由 \(h_{\theta}(x^{(i)})=(x^{(i)})^{T}\theta\)，計算預測值和實際值的差向量：
\[X \theta -\vec{y}=\left[ \begin{matrix} (x^{(1)})^{T}\theta \\ . \\ . \\ . \\ (x^{(m)})^{T}\theta \end{matrix} \right]-\left[ \begin{matrix} (y^{(1)}) \\ . \\ . \\ . \\ (y^{(m)}) \end{matrix} \right]=\left[ \begin{matrix} h_{\theta}(x^{(1)})-(y^{(1)}) \\ . \\ . \\ . \\ h_{\theta}(x^{(m)})-(y^{(m)}) \end{matrix} \right]\]
由 \(z^{T}z=\sum_{i} z^{2}_{i}\)，得：
\[\frac{1}{2}(X\theta -\vec{y})^{T}(X\theta-\vec{y})=\frac{1}{2}\sum_{i=1}^{m}(h_{\theta}(x^{(i)}-y^{(i)}))^{2}=J(\theta)\]
因此
\[J(\theta)=\frac{1}{2}(X\theta -\vec{y})^{T}(X\theta-\vec{y})\]
而後就是求\(J(\theta)\)關於\(\theta\)的導數，並使其等於 0，最後就能夠求出\(\theta\)。\(J(\theta)\) 求導過程以下。
上面的推導過程看着很複雜，但稍微耐心看下每一步，就會以爲瓜熟蒂落了。具體每一步的推導依據可參考機器學習課件第 11 頁。
爲最小化 \(J(\theta)\)，咱們置公式(6)爲 0，得：
\[X^{T}X\theta=X^{T}\vec{y}\]
由此得
\[\theta =(X^{T}X)^{-1}X^{T}\vec{y}\]
至此，求解完畢。
代碼以下。
import statsmodels.api as sm import statsmodels.formula.api as smf import statsmodels.graphics.api as smg import patsy %matplotlib inline import matplotlib.pyplot as plt import numpy as np import pandas as pd from scipy import stats import seaborn as sns
先生成訓練數據集。
np.random.seed(123456789)
N = 100
x1 = np.random.randn(N) # 自變量
x2 = np.random.randn(N)
def y_true(x1, x2):
return 1 + 2 * x1 + 3 * x2
ytrue = y_true(x1, x2)
y = ytrue + e
X = np.vstack([np.ones(N), x1, x2]).T
先用 numpy 自帶的最小二乘法包 np.linalg.lstsq 求 \(\theta\)值。
beta, res, rank, sval = np.linalg.lstsq(X, y) print beta
[ 0.86896993 1.98195932 2.96540985]
再使用咱們上面的矩陣方法 \(\theta =(X^{T}X)^{-1}X^{T}\vec{y}\) 求解。
def matrixLSE(X, y): # 矩陣法求最小二乘
xtx = np.dot(X.T, X)
inverse = np.linalg.inv(xtx)
xxx = np.dot(inverse, X.T)
thetalse = np.dot(xxx, y)
return thetalse
thetalse = matrixLSE(X, y) print thetalse
[ 0.86896993 1.98195932 2.96540985]
求出的 \(\theta\) 值跟 numpy.linalg.lstsq 求的解同樣，說明咱們的代碼是正確的~
2. 使用pandas庫中的函數下載上證綜指和任一成份股票數據，計算日收益率，對這兩組數據創建迴歸模型，將上證綜指的收益率做爲解釋變量，說明這個模型的用處。
from datetime import datetime from pandas_datareader import data, wb # 須要先安裝 pandas_datareader
獲取上證綜指數據。
start = datetime(2015, 6, 3)
end = datetime(2016, 6, 3)
Shanghai = data.DataReader("000001.SS", 'yahoo', start, end) # 000001.SS 表示上證綜指，返回 DataFrame
Shanghai.head() # 縱軸是日期，橫軸是開盤價、最高價、最低價、收盤價、成交量、復權收盤。因上證綜指並不是具體某支股票，因此交易量爲 0。
Shanghai = Shanghai.drop('Open',axis=1)
Shanghai = Shanghai.drop('High',axis=1)
Shanghai = Shanghai.drop('Low',axis=1)
Shanghai = Shanghai.drop('Volume',axis=1)
Shanghai = Shanghai.drop('Adj Close',axis=1)
Shanghai["retRate"] = (Shanghai.Close-Shanghai.Close.shift(1))/Shanghai.Close.shift(1) # 日收益率，固然在作實證時都用對數收益率
Shanghai = Shanghai.ix[1:]
Shanghai.head()
| Close | retRate | |
|---|---|---|
| Date | ||
| 2015-06-04 | 4947.10 | 0.007560 |
| 2015-06-05 | 5023.10 | 0.015363 |
| 2015-06-08 | 5131.88 | 0.021656 |
| 2015-06-09 | 5113.53 | -0.003576 |
| 2015-06-10 | 5106.04 | -0.001465 |
獲取宇通客車的數據。
start = datetime(2015, 6, 3)
end = datetime(2016, 6, 3)
Yutong = data.DataReader("600066.SS", 'yahoo', start, end) # 宇通客車，返回 DataFrame
Yutong.head() # 縱軸是日期，橫軸是開盤價、最高價、最低價、收盤價、成交量、復權收盤。因上證綜指並不是具體某支股票，因此交易量爲 0。
Yutong = Yutong.drop('Open',axis=1)
Yutong = Yutong.drop('High',axis=1)
Yutong = Yutong.drop('Low',axis=1)
Yutong = Yutong.drop('Volume',axis=1)
Yutong = Yutong.drop('Adj Close',axis=1)
Yutong["retRate"] = (Yutong.Close-Yutong.Close.shift(1))/Yutong.Close.shift(1) # 日收益率，固然在作實證時都用對數收益率
Yutong = Yutong.ix[1:]
Yutong.head()
| Close | retRate | |
|---|---|---|
| Date | ||
| 2015-06-04 | 23.17 | -0.022775 |
| 2015-06-05 | 23.41 | 0.010358 |
| 2015-06-08 | 23.11 | -0.012815 |
| 2015-06-09 | 23.07 | -0.001731 |
| 2015-06-10 | 22.83 | -0.010403 |
使用上面介紹的 smf.ols 作線性迴歸。
x = Shanghai["retRate"].values
y = Yutong.ix[Shanghai.index].retRate.values
data = pd.DataFrame({"x": x, "y": y})
model = smf.ols("y ~ x", data) # 這裏沒有寫截距，截距是默認存在的
result = model.fit()
print(result.summary())
OLS Regression Results
==============================================================================
Dep. Variable: y R-squared: 0.542
Model: OLS Adj. R-squared: 0.540
Method: Least Squares F-statistic: 270.3
Date: Sat, 18 Jun 2016 Prob (F-statistic): 1.38e-40
Time: 14:54:28 Log-Likelihood: 578.89
No. Observations: 230 AIC: -1154.
Df Residuals: 228 BIC: -1147.
Df Model: 1
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [95.0% Conf. Int.]
------------------------------------------------------------------------------
Intercept 0.0011 0.001 0.845 0.399 -0.001 0.004
x 0.8362 0.051 16.442 0.000 0.736 0.936
==============================================================================
Omnibus: 24.051 Durbin-Watson: 1.863
Prob(Omnibus): 0.000 Jarque-Bera (JB): 42.013
Skew: 0.578 Prob(JB): 7.53e-10
Kurtosis: 4.746 Cond. No. 39.3
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
看自變量的 P 值，Intercept 能夠去除，只由 x 預測。
R-squared 只有 0.542，數據的擬合不夠好，說明用上證綜指預測宇通客車的日收益率不靠譜。
線性迴歸，是利用數理統計中迴歸分析，來肯定兩種或兩種以上變量間相互依賴的定量關係的一種統計分析方法，運用十分普遍。線性迴歸假設特徵和結果知足線性關係。線性關係的表達能力很是強大，每一個特徵對結果的影響強弱能夠由前面的參數體現，並且每一個特徵變量能夠首先映射到一個函數，而後再參與線性計算。這樣就能夠表達特徵與結果之間的非線性關係。
3. 在kaggle網站上找到titanic數據，並使用 logistic 迴歸來建模，研究分析每一個因素對應生存的重要性。
數據包含的字段以下：
- PassengerID
- Survived(存活與否)
- Pclass（客艙等級）
- Name（姓名）
- Sex（性別）
- Age（年齡）
- SibSp（親戚和配偶在船數量）
- Parch（父母孩子的在船數量）
- Ticket（票編號）
- Fare（價格）
- Cabin（客艙位置）
- Embarked（上船的港口編號）
數據讀取
input_df = pd.read_csv('Data/train.csv', header=0)
input_df.head()
| PassengerId | Survived | Pclass | Name | Sex | Age | SibSp | Parch | Ticket | Fare | Cabin | Embarked | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 0 | 3 | Braund, Mr. Owen Harris | male | 22.0 | 1 | 0 | A/5 21171 | 7.2500 | NaN | S |
| 1 | 2 | 1 | 1 | Cumings, Mrs. John Bradley (Florence Briggs Th... | female | 38.0 | 1 | 0 | PC 17599 | 71.2833 | C85 | C |
| 2 | 3 | 1 | 3 | Heikkinen, Miss. Laina | female | 26.0 | 0 | 0 | STON/O2. 3101282 | 7.9250 | NaN | S |
| 3 | 4 | 1 | 1 | Futrelle, Mrs. Jacques Heath (Lily May Peel) | female | 35.0 | 1 | 0 | 113803 | 53.1000 | C123 | S |
| 4 | 5 | 0 | 3 | Allen, Mr. William Henry | male | 35.0 | 0 | 0 | 373450 | 8.0500 | NaN | S |
數據預處理
input_df.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 891 entries, 0 to 890 Data columns (total 12 columns): PassengerId 891 non-null int64 Survived 891 non-null int64 Pclass 891 non-null int64 Name 891 non-null object Sex 891 non-null object Age 714 non-null float64 SibSp 891 non-null int64 Parch 891 non-null int64 Ticket 891 non-null object Fare 891 non-null float64 Cabin 204 non-null object Embarked 889 non-null object dtypes: float64(2), int64(5), object(5) memory usage: 83.6+ KB
能夠看到 Age、Cabin、Embarked 幾個字段有空值，填充下。其中客艙位置 Cabin 空值最多，沒什麼用，乾脆去掉。
input_df = input_df.drop('Cabin', axis=1) # 去除沒用的字段客艙位置，空值還多
input_df.ix[input_df['Age'].isnull(),'Age'] = input_df['Age'].median() # 年齡取箇中間值
input_df.ix[input_df['Embarked'].isnull(),'Embarked'] = input_df['Embarked'].mode().values # 上船的港口編號
Name 也是一個沒用的字段，不信什麼名字能影響生死。固然有人會說，有些姓名是貴族纔有的，那還有價格 Fare 和客艙等級 Pclass 體現呢，從名字來識別身份判斷存亡最不靠譜，去掉。
input_df = input_df.drop('Name', axis=1)
票編號 Ticket 也去掉，跟名字 Name 同樣沒用。
input_df = input_df.drop('Ticket', axis=1)
input_df.Embarked.unique() # 上船的港口編號，沒用的字段，果斷去掉
array(['S', 'C', 'Q'], dtype=object)
input_df = input_df.drop('Embarked', axis=1)
input_df = input_df.drop('PassengerId', axis=1)
input_df.Sex = input_df.Sex.map({"male": 1, "female": 0}) # 作個映射
input_df.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 891 entries, 0 to 890 Data columns (total 7 columns): Survived 891 non-null int64 Pclass 891 non-null int64 Sex 891 non-null int64 Age 891 non-null float64 SibSp 891 non-null int64 Parch 891 non-null int64 Fare 891 non-null float64 dtypes: float64(2), int64(5) memory usage: 48.8 KB
input_df.head()
| Survived | Pclass | Sex | Age | SibSp | Parch | Fare | |
|---|---|---|---|---|---|---|---|
| 0 | 0 | 3 | 1 | 22.0 | 1 | 0 | 7.2500 |
| 1 | 1 | 1 | 0 | 38.0 | 1 | 0 | 71.2833 |
| 2 | 1 | 3 | 0 | 26.0 | 0 | 0 | 7.9250 |
| 3 | 1 | 1 | 0 | 35.0 | 1 | 0 | 53.1000 |
| 4 | 0 | 3 | 1 | 35.0 | 0 | 0 | 8.0500 |
由於邏輯迴歸建模時，須要輸入的特徵都是數值型特徵，一般會先對類目型的特徵因子化。
dummies_Sex = pd.get_dummies(input_df['Sex'], prefix= 'Sex') dummies_Pclass = pd.get_dummies(input_df['Pclass'], prefix= 'Pclass') input_df = pd.concat([input_df, dummies_Sex, dummies_Pclass], axis=1) input_df.drop(['Pclass','Sex'], axis=1, inplace=True) input_df.head()
| Survived | Age | SibSp | Parch | Fare | Sex_0 | Sex_1 | Pclass_1 | Pclass_2 | Pclass_3 | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 22.0 | 1 | 0 | 7.2500 | 0.0 | 1.0 | 0.0 | 0.0 | 1.0 |
| 1 | 1 | 38.0 | 1 | 0 | 71.2833 | 1.0 | 0.0 | 1.0 | 0.0 | 0.0 |
| 2 | 1 | 26.0 | 0 | 0 | 7.9250 | 1.0 | 0.0 | 0.0 | 0.0 | 1.0 |
| 3 | 1 | 35.0 | 1 | 0 | 53.1000 | 1.0 | 0.0 | 1.0 | 0.0 | 0.0 |
| 4 | 0 | 35.0 | 0 | 0 | 8.0500 | 0.0 | 1.0 | 0.0 | 0.0 | 1.0 |
Age 和Fare 兩個屬性幅度變化太大了吧，對於邏輯迴歸與梯度降低來講，各屬性值之間scale差距太大，將對收斂速度形成嚴重損害，甚至不收斂。這裏用scikit-learn 裏面的 preprocessing 模塊對這兩個屬性作一個縮放 scaling ，所謂 scaling，其實就是將一些變化幅度較大的特徵化到[-1,1]以內。
import warnings
warnings.filterwarnings("ignore") # 屏蔽掉討厭的警告
import sklearn.preprocessing as preprocessing
scaler = preprocessing.StandardScaler()
age_scale_param = scaler.fit(input_df['Age'])
input_df['Age_scaled'] = scaler.fit_transform(input_df['Age'], age_scale_param)
fare_scale_param = scaler.fit(input_df['Fare'])
input_df['Fare_scaled'] = scaler.fit_transform(input_df['Fare'], fare_scale_param)
input_df.head()
| Survived | Age | SibSp | Parch | Fare | Sex_0 | Sex_1 | Pclass_1 | Pclass_2 | Pclass_3 | Age_scaled | Fare_scaled | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 22.0 | 1 | 0 | 7.2500 | 0.0 | 1.0 | 0.0 | 0.0 | 1.0 | -0.565736 | -0.502445 |
| 1 | 1 | 38.0 | 1 | 0 | 71.2833 | 1.0 | 0.0 | 1.0 | 0.0 | 0.0 | 0.663861 | 0.786845 |
| 2 | 1 | 26.0 | 0 | 0 | 7.9250 | 1.0 | 0.0 | 0.0 | 0.0 | 1.0 | -0.258337 | -0.488854 |
| 3 | 1 | 35.0 | 1 | 0 | 53.1000 | 1.0 | 0.0 | 1.0 | 0.0 | 0.0 | 0.433312 | 0.420730 |
| 4 | 0 | 35.0 | 0 | 0 | 8.0500 | 0.0 | 1.0 | 0.0 | 0.0 | 1.0 | 0.433312 | -0.486337 |
邏輯迴歸
使用 statsmodels.formula.api 的 logit 接口進行邏輯斯特迴歸。
model = smf.logit("Survived ~ Pclass_1 + Pclass_2 + Pclass_1 + Sex_0 + Sex_1 + Age_scaled + SibSp + Parch + Fare_scaled", data=input_df)
result = model.fit() # 內部使用極大似然估計，因此會有結果返回，表示成功收斂
print(result.summary())
Optimization terminated successfully.
Current function value: 0.442813
Iterations 8
Logit Regression Results
==============================================================================
Dep. Variable: Survived No. Observations: 891
Model: Logit Df Residuals: 883
Method: MLE Df Model: 7
Date: Mon, 20 Jun 2016 Pseudo R-squ.: 0.3350
Time: 02:03:31 Log-Likelihood: -394.55
converged: True LL-Null: -593.33
LLR p-value: 7.950e-82
===============================================================================
coef std err z P>|z| [95.0% Conf. Int.]
-------------------------------------------------------------------------------
Intercept -0.5152 7.34e+06 -7.02e-08 1.000 -1.44e+07 1.44e+07
Pclass_1 2.1519 0.290 7.420 0.000 1.584 2.720
Pclass_2 1.1400 0.228 5.006 0.000 0.694 1.586
Sex_0 1.1216 7.34e+06 1.53e-07 1.000 -1.44e+07 1.44e+07
Sex_1 -1.6368 7.34e+06 -2.23e-07 1.000 -1.44e+07 1.44e+07
Age_scaled -0.5096 0.102 -4.997 0.000 -0.709 -0.310
SibSp -0.3481 0.109 -3.192 0.001 -0.562 -0.134
Parch -0.1084 0.117 -0.923 0.356 -0.339 0.122
Fare_scaled 0.1499 0.122 1.234 0.217 -0.088 0.388
===============================================================================
僞 R 方 Pseudo R-squ. 只有 0.3350，效果不行。
看各自變量的 P 值，客艙等級 Pclass、年齡 Age 和 親戚和配偶在船數量 SibSp 要影響要顯著些。性別 Sex 居然沒影響，不是女士優先嗎？價格 Fare 和 父母孩子的在船數量 Parch 的影響能夠忽略。
用這個 model 作下預測。
讀取 test.csv，和 train.csv 作同樣的預處理。
testAll_df = pd.read_csv('Data/test.csv', header=0)
test_df = testAll_df.drop('Cabin', axis=1) # 去除沒用的字段客艙位置，空值還多
test_df.ix[test_df['Age'].isnull(),'Age'] = test_df['Age'].median() # 年齡取箇中間值
test_df.ix[test_df['Fare'].isnull(),'Fare'] = test_df['Fare'].median() # 票價取箇中間值
test_df.ix[test_df['Embarked'].isnull(),'Embarked'] = test_df['Embarked'].mode().values # 上船的港口編號
test_df = test_df.drop('Name', axis=1)
test_df = test_df.drop('Ticket', axis=1)
test_df = test_df.drop('Embarked', axis=1)
test_df.Sex = test_df.Sex.map({"male": 1, "female": 0}) # 作個映射
dummies_Sex = pd.get_dummies(test_df['Sex'], prefix= 'Sex')
dummies_Pclass = pd.get_dummies(test_df['Pclass'], prefix= 'Pclass')
test_df = pd.concat([test_df, dummies_Sex, dummies_Pclass], axis=1)
test_df.drop(['Pclass','Sex'], axis=1, inplace=True)
import warnings
warnings.filterwarnings("ignore") # 屏蔽掉討厭的警告
import sklearn.preprocessing as preprocessing
scaler = preprocessing.StandardScaler()
age_scale_param = scaler.fit(test_df['Age'])
test_df['Age_scaled'] = scaler.fit_transform(test_df['Age'], age_scale_param)
fare_scale_param = scaler.fit(test_df['Fare'])
test_df['Fare_scaled'] = scaler.fit_transform(test_df['Fare'], fare_scale_param)
test_df.head()
| PassengerId | Age | SibSp | Parch | Fare | Sex_0 | Sex_1 | Pclass_1 | Pclass_2 | Pclass_3 | Age_scaled | Fare_scaled | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 892 | 34.5 | 0 | 0 | 7.8292 | 0.0 | 1.0 | 0.0 | 0.0 | 1.0 | 0.386231 | -0.497413 |
| 1 | 893 | 47.0 | 1 | 0 | 7.0000 | 1.0 | 0.0 | 0.0 | 0.0 | 1.0 | 1.371370 | -0.512278 |
| 2 | 894 | 62.0 | 0 | 0 | 9.6875 | 0.0 | 1.0 | 0.0 | 1.0 | 0.0 | 2.553537 | -0.464100 |
| 3 | 895 | 27.0 | 0 | 0 | 8.6625 | 0.0 | 1.0 | 0.0 | 0.0 | 1.0 | -0.204852 | -0.482475 |
| 4 | 896 | 22.0 | 1 | 1 | 12.2875 | 1.0 | 0.0 | 0.0 | 0.0 | 1.0 | -0.598908 | -0.417492 |
predictions = result.predict(test_df)
Survivedresult = pd.DataFrame({'PassengerId':testAll_df['PassengerId'].as_matrix(), 'Survived':predictions.astype(np.int32)})
Survivedresult.to_csv("Data/logistic_regression_predictions.csv", index=False)
最後在 Kaggle 的 Make a submission 頁面，提交上結果。
4. 蒐集某個城市過去1個月的PM2.5小時級數據，根據時間序列預測方法，來預測將來三天的PM2.5狀況。
# 數據讀取並清理
def readpmfile(filepath):
beijing_aq = pd.read_csv(filepath, header=0) # header 指明用做列名的行號
beijing_pm25 = beijing_aq.ix[beijing_aq['type']=='PM2.5',[0,1,3]] # 抽取 PM2.5 數據子集，就選東四的 PM2.5 數據
beijing_pm25.rename(columns={"東四": "Dongsi"}, inplace=True) # 列重命名
beijing_pm25['datetime'] = pd.to_datetime(beijing_pm25.date, format='%Y%m%d')
beijing_pm25.index = range(24) # 修改索引，以便後面把 datetime 加上 hour
length = beijing_pm25.shape[0] # 獲取數據個數
for i in range(0, length): # 由於 datetime.timedelta 不能接受 Series 參數，因此只能一個一個加
beijing_pm25.ix[i,'datetime'] = beijing_pm25.ix[i,'datetime'] + datetime.timedelta(hours = beijing_pm25.ix[i,'hour'])
beijing_pm25 = beijing_pm25.drop(['date','hour'], axis = 1) # 丟掉這兩個已經沒用的字段
beijing_pm25 = beijing_pm25.set_index("datetime") # 用 datetime 來做爲索引
return beijing_pm25
import datetime
beijing_pm = pd.DataFrame() # 每次先清空，省得屢次執行形成重複加載
starttime= datetime.datetime(2016, 5, 1) # 讀取從 5 月 1 日開始的 34 天數據，5 月 31 天，6 月 3 天
for i in range(0, 34): # 天天一個 pm2.5 數據文件，文件須要一個一個讀取
curtime = starttime + datetime.timedelta(days = i)
filepath = 'Data/beijing_20160101-20160611/'
filename = 'beijing_all_' + curtime.strftime('%Y%m%d') + '.csv'
beijing_pm25 = readpmfile(filepath + filename)
beijing_pm = pd.concat([beijing_pm, beijing_pm25])
beijing_pm.info() # 天天 24 條 pm2.5 數據，34 天，應有 24*34=816 條數據，這裏只有 792 條是非空的
<class 'pandas.core.frame.DataFrame'> DatetimeIndex: 816 entries, 2016-05-01 00:00:00 to 2016-06-03 23:00:00 Data columns (total 1 columns): Dongsi 792 non-null float64 dtypes: float64(1) memory usage: 12.8 KB
處理缺失數據，注意這裏不能把數據刪除，要填充，不然後面 predict 會出問題。血的教訓，切記~！
#beijing_pm = beijing_pm.dropna() # 缺失數據直接過濾掉 beijing_pm = beijing_pm.fillna(method = 'pad') # 用前一個數據代替NaN beijing_pm.head()
| Dongsi | |
|---|---|
| datetime | |
| 2016-05-01 00:00:00 | 213.0 |
| 2016-05-01 01:00:00 | 202.0 |
| 2016-05-01 02:00:00 | 172.0 |
| 2016-05-01 03:00:00 | 172.0 |
| 2016-05-01 04:00:00 | 172.0 |
pm_May = beijing_pm[beijing_pm.index.month == 5] # 取 5 月份的 PM2.5 數據作訓練 pm_June = beijing_pm[beijing_pm.index.month == 6] # 預測 6 月份的 PM2.5 含量
pm_May.plot(figsize=(12, 4)); # 看看 5 月份的數據
可見 PM2.5 數據仍是有明顯的規律，有規律才能夠建模型。規律體如今什麼地方，畫個散點圖，看必定時間先後的 PM2.5 有沒有相關性。
plt.scatter(pm_May[1:], pm_May[:-1]); # 滯後一小時，看前一小時的 PM2.5 數據對後一小時有沒有影響
從圖中看，是呈線性相關的，可見一個小時仍是有必定影響的。滯後項的相關性稱爲自相關，本身的歷史數據會影響後來的數據，這種規律稱爲自相關。
plt.scatter(pm_May[2:], pm_May[:-2]); # 間隔兩小時還有沒有影響
影響比一個小時弱點，很容易理解。
再來相鄰兩天的數據有沒有相關性。
plt.scatter(pm_May[24:], pm_May[:-24]); # 滯後 24 小時
影響不太明顯了。
通常來講，自相關就是越近的數據影響越大。能夠畫相關係數圖。建模時要把這種自相關性去掉，使用差分來去掉，第二張圖就使用了差分，差分以後，自相關性有所減弱，第三張圖是差分的差分，第四張圖是三階差分，差分階次越高，自相關性就越小。
fig, axes = plt.subplots(1, 4, figsize=(12, 3)) smg.tsa.plot_acf(pm_May.Dongsi, lags=24, ax=axes[0]) # 時間序列子模塊，自相關，今天的跟歷史 24 個小時的相關性 smg.tsa.plot_acf(pm_May.Dongsi.diff().dropna(), lags=24, ax=axes[1]) # 每一個條狀是個相關係數 smg.tsa.plot_acf(pm_May.Dongsi.diff().diff().dropna(), lags=24, ax=axes[2]) smg.tsa.plot_acf(pm_May.Dongsi.diff().diff().diff().dropna(), lags=24, ax=axes[3]) fig.tight_layout()
使用自迴歸模型 sm.tsa.AR 來作時間序列預測。
model = sm.tsa.AR(pm_May.Dongsi) # AR 模型就是自迴歸模型 result = model.fit(24) #認爲過去 24 小時都會影響如今時間點的 PM2.5 濃度 sm.stats.durbin_watson(result.resid) # 把殘差拿出來檢驗下
1.9984097620329431
殘差是不該該有自迴歸的，把殘差拿出來檢驗下，在 2 附近，意味着沒有自迴歸。
fig, ax = plt.subplots(1, 1, figsize=(8, 3)) smg.tsa.plot_acf(result.resid, lags=24, ax=ax) # 畫殘差自相關性圖 fig.tight_layout()
殘差沒有自相關性。
下面作預測，預測 6 月份頭三天的 PM2.5 數據。
fig, ax = plt.subplots(1, 1, figsize=(12, 4))
ax.plot(pm_May.index.values[-72:], pm_May.Dongsi.values[-72:], label="train data")
ax.plot(pm_June.index.values, pm_June.Dongsi.values, label="actual outcome")
ax.plot(pd.date_range("2016-06-01", "2016-06-04", freq="H").values,
result.predict("2016-06-01", "2016-06-04"), label="predicted outcome")
ax.legend()
fig.tight_layout()
綠色是 6 月份真實值，紅色是預測值。預測效果太差了！有待仔細檢查。
- 1. 線性迴歸——練習題
- 2. 線性模型
- 3. 從線性到非線性模型-對數線性模型
- 4. python題型練習
- 5. 泛型練習題
- 6. 【C++練習題】實現單例模型
- 7. Tensorflow學習（練習）—CPU訓練模型
- 8. Chap03 線性模型
- 9. 02線性模型
- 10. chap3-線性模型
- 更多相關文章...
- • ASP.NET MVC - 模型 - ASP.NET 教程
- • SVG 漸變 - 線性 - SVG 教程
- • Kotlin學習（二）基本類型
- • 委託模式
-
每一个你不满意的现在，都有一个你没有努力的曾经。
- 1. 線性迴歸——練習題
- 2. 線性模型
- 3. 從線性到非線性模型-對數線性模型
- 4. python題型練習
- 5. 泛型練習題
- 6. 【C++練習題】實現單例模型
- 7. Tensorflow學習（練習）—CPU訓練模型
- 8. Chap03 線性模型
- 9. 02線性模型
- 10. chap3-線性模型