C1 Introduction to Statistical Learning
1.1Statistical Learning介紹:
1.Statistical learning
a vast set of tools ( supervised or unsupervised ) for estimating Model. supervised 為了預測，unsupervised 為了發現。
2.variables
input variables: predictors, independent variables, features, or sometimes just variables
output variables: response or dependent variable
3.\(Y =f(X)+\epsilon\)
Here f is some fixed but unknown function of \(X_1\), . . . , \(X_p\), and ε is a random error term, which is independent of \(X\) and has mean zero. 由於 \(\epsilon\) 代表了所有的錯誤，如模型形式錯誤，遺漏重要predictors，測量錯誤等，所以假設其實難以保證。
4.權衡 Accuracy 和 Interpretability
5.supervised, unsupervised and semi-supervised
6.Regression Versus Classification
1.1.1 估計 \(f\) 的目的：prediction和/或inference。
1.prediction
(1)\(\hat{Y} = \hat{f}(X)\) where \(\hat{f}\) represents our estimate for \(f\) , and \(\hat{Y}\) represents the resulting prediction for \(Y\). \(\hat{f}\) is often treated as a black box. 只要預測准就行，不關心模型的形式。
(2)reducible error (\(f\)) and the irreducible error (\(ε\))：下面式子在\(\hat{f}\)和\(X\)固定的情況下得出
目標是減少reducible error，而irreducible error為模型的准確度設置了我們無法知道的上限。
2.inference
understand how \(Y\) changes as a function of \(X_1\), . . . , \(X_p\) ，需要知道模型的形式
目的是知道哪些predictors與response顯著相關；response和每一個predictor的關系；模型是否能線性解釋
根據估計 \(f\) 的目的選擇相應的模型。
1.1.2 估計 \(f\) 的方法：parametric 或 non-parametric
1.parametric(model-based approach)
步驟：假設 functional form, or shape, of $ f$ (eg. linear), then use the training data to fit or train the model
特點：將估計整個 \(f\) 變為估計參數，問題對模型估計不夠准確。選擇flexible模型的話能提高准確，但可能overfitting
2.Non-parametric
步驟：不假設 functional form，直接訓練模型
特點：沒有簡化問題，提高對模型估計的准確性，但需要更多訓練數據
1.2 評估模型准確性
1.2.1 回歸的評估
以mean squared error (MSE) 為例，假設真正的 \(f\) 是下左圖的黑線，處於某個適中的flexibility，模型的flexibility從最小值開始，test set遠大於training set，則隨着flexibility的增加，train MSE會不斷減少，test MSE會先減少后增加，最小值最接近\(Var(ε)\)。
上圖橙和綠分別代表underfitting和overfitting。由於training set不夠大，overfitting使得模型發現了並不存在於test set的規律（noise部分）。
1.2.2 Bias-Variance的平衡
test MSE的期望可分解如下：
所以，為了最小化expected test MSE，需要減少variance and bias。
前者表示利用不同training set評估模型時產生的差異。一般情況下，越flexible的模型，variance越大（如1.2.1右圖中，綠色模型的test set和training set的MSE的差距很大，橙色很小）
后者表示模型擬合不足
正如1.2.1右圖中的test MSE，隨着flexibility的增加，variance / bias先減少后增加。
1.2.3 分類的情況
error rate
1.The Bayes Classifier
基於predictor值，將每一個observation推測為最有可能的class。即使得 $\mathrm{Pr}(Y = j\ |\ X = x_0) $ 這個概率最大。例如，在binary classification問題上，如果 \(\mathrm{Pr}(Y = 1\ |\ X = x_0) > 0.5\) 則預測class one，否則預測class two。
Bayes Classifier 可以產生 lowest possible test error rate, Bayes error rate: \(1−E( \mathrm{max}_j\mathrm{Pr}(Y =j|X) )\)
2.K-Nearest Neighbors (注意要scale)
在現實中，由於無法得知在 \(X\) 條件下 \(Y\) 的分布，所以無法計算Bayes Classifier。而KNN classifier 是一個很接近Bayes Classifier的方法。
如下左圖所示，假設K=3，藍色和橙色圈分別代表class 1和class 2，交叉為 \(x_0\) ，這里使得Pr最大的是j = class 1。右圖為K=3時的貝葉斯邊界。
KNN Classifier 和 Bayes Classifier的接近程度部分取決於K的選擇，K越小，flexibility越大，反之越小。和回歸情況一樣，當K=1時會overfitting，training error rate等於0，但test error rate很高。
C2 Linear Regression
本章利用Advertising數據，如下所示，三個媒體的數值為廣告投入費用：
| row_id | TV | Radio | Newspaper | Sales |
|---|---|---|---|---|
| 0 | 230.1 | 37.8 | 69.2 | 22.1 |
要回答的問題：
- Sales和廣告投入是否有關系：多元回歸 + F-statistics
- 關系有多大：RES/mean_response or R^2^
- 哪些媒體跟sales有關系：t-statistics or p-value
- 每個媒體的單位投入對應多少sales的變化，多准確：\(\hat{\beta_j}\) 的置信區間
- 如何准確預測：prediction interval（預測individual response）或confidence intervals（預測average response）
- 關系是否為線性：residual plots
2.1 簡單線性回歸
$Y ≈ β_0 + β_1X $ 表示regressing Y on X (or Y onto X)
\(\hat{y} = \hat{β}_0 + \hat{β}_1 x\) least squares line, where \(\hat{y}\) indicates a prediction of \(Y\) on the basis of \(X = x\). 目的是讓 \(y_i ≈ \hat{y}_i\) ，而從全部數據來看，就是要最小化RSS。
\(\mathrm{RSS} = e^2_1 + e^2_2 + · · · + e^2_n\) residual sum of squares，其中\(e_i =y_i - \hat{y}_i\)
通過計算 \(\hat{β}_0\) 和 \(\hat{β}_1\) 使得RSS取最小值。
#計算，進行中心化，但不標准化；reshape(-1,1)就是把X變為n*1的形狀，由於n不知道，可用-1讓python自己算
X = scale(advertising.TV, with_mean=False, with_std=False).reshape(-1,1)
y = advertising.Sales
regr = skl_lm.LinearRegression()
regr.fit(X,y)
print(regr.intercept_) #7.032593549127693
print(regr.coef_) #[0.04753664]
#畫圖，order表示flexibility，ci表示置信水平[0,100]（數據量大時建議None），s表示點的大小
sns.regplot(advertising.TV, advertising.Sales, order=1, ci=None, scatter_kws={'color':'r', 's':9})
plt.xlim(-10,310)
plt.ylim(ymin=0)
population regression line: $Y = β_0 + β_1X + ε $ 假設Y和X是線性關系，該線是最好的線性回歸。下圖的紅線代表population regression line，各點就是根據它加上一個 ε（服從正態分布）得出的，而藍線是上面提到的least squares line，即通過那些數據點估計出的線。
由於根據特定的數據集來推斷，所以難免會有bias。但在模型設定正確的情況下，通過大量數據，所推斷出的least squares line會十分接近 population regression line，會是unbiased。如上右圖所示，如果平均各線，結果不會與真線相差太大。
2.1.1 評估系數的准確性
估計值與真實值的平均差異（標准差，但無法知道多小才算小）；在某個置信水平下，系數值的范圍是多少（置信區間）；關系是否顯著（null hypothesis）。
1.類比“利用sample mean估計population mean的准確性”：
$\mathrm{Var}(\hat{μ}) = \mathrm{SE}(\hat{μ})^2 = \frac{\sigma^2}{n} $ where σ is the standard deviation of each of the realizations \(y_i\) of Y，要求n個觀測值各不相關。
the standard error \(\mathrm{SE}(\hat{μ})\) 表示估計值與真實值的平均差異，當n越大，差異越小，即估計值越准確。
2.類似的，系數的式子如下：
\(\mathrm{SE}(\hat{β}_0)^2 =σ^2[\frac{1}{n}+\frac{\bar{x}^2}{\sum_{i=1}^{n}(x_i-\bar{x})^2}]\) , \(\mathrm{SE}(\hat{β}_1)^2 =\frac{\sigma^2}{\sum_{i=1}^{n}(x_i-\bar{x})^2}\) 其中\(σ^2 = \mathrm{Var}(ε)\) 。這兩個式子嚴格成立的話，\(\epsilon_i\) 和 \(\sigma^2\) 是uncorrelated的。雖然這規定不符合，但是結果還是可以接受的。從式子可以知道，如果x很分散，那么需要更關注斜率；如果x的均值等於0，那么 \(\hat{β}_0 = \bar{y}\) 。通常來說，\(\sigma\) 是要通過樣本推斷的，式子為 \(\mathrm{RSE} = \sqrt{\mathrm{RSS}/(n − 2)}\) residual standard error
利用標准差可以計算置信區間。例如 \(\hat{β}_1\) 的置信區間為\(\hat{β}_1±2·\mathrm{SE}(\hat{β_1})\) ，意思大概是通過100次樣本估計，有95次這些區間包含了未知的真值。這個式子假設error為高斯分布，而2會隨n的變化而變化，更加准確應該用t分布。
上面回歸的系數，在95%置信水平下，區間為 \(β_0\) [6.130, 7.935] 和 \(β_1\) [0.042, 0.053] 。在沒有廣告的情況下，sales的范圍是[6.130, 7.935]，而每增加一單位的電視廣告投入，sales增加的范圍是[0.042, 0.053]。
利用標准差還可以進行假設檢驗，例如下面是針對 \(\hat{\beta_1}\) 的null hypothesis：\(H_0: \beta_1=0\) , \(t=\frac{\hat{\beta_1}-0}{\mathrm{SE}(\hat{\beta}_1)}\) 。從式子可以看出，如果系數值越大，且它的標准差越小，那么該系數就越不可能為0，即該系數的predictors跟response的關系大且顯著。或者說p值越小（拒絕原假設），關系越顯著。
import statsmodels.formula.api as smf
est = smf.ols('Sales ~ TV', advertising).fit()
est.summary().tables[1]
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 7.0326 | 0.458 | 15.360 | 0.000 | 6.130 | 7.935 |
| TV | 0.0475 | 0.003 | 17.668 | 0.000 | 0.042 | 0.053 |
2.1.2 評估模型的准確性
指標為the residual standard error (RSE) 和 the R2
# RSE的計算
import math
n = advertising.TV.count()
rss = ((advertising.Sales - est.predict())**2).sum()
rse = math.sqrt(rss/(n-2)) # 3.26
#R2和F-test
est.summary().tables[0] # R2 0.612, F-test 312.1
1.RSE測量模型the lack of fit部分
RSE=3.26意味着sales有3.26的波動，大小可通過RSE/(sales mean)判斷，此數據結果為23%。
2.R^2^
式子為 $1− \frac{\mathrm{RSS}}{\mathrm{TSS}} $ ，它是一個測量X和Y線性關系的指標。在復雜的問題里，R^2^小於0.1也是很正常的。當單變量進行線性回歸時，R^2^是等於Cor^2^(X, Y)的，但后者不適用於多變量。
2.2 多元線性回歸
簡單線性回歸當中，如果不同回歸中的predictors是相關的，這會產生誤導。
當把TV，Radio，Newspaper三個predictors合在一起回歸\(\hat{y} = \hat{β}_0 + \hat{β}_1x_1 + \hat{β}_2x_2 + \hat{β}_3x_3\)，結果如下
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 2.9389 | 0.312 | 9.422 | 0.000 | 2.324 | 3.554 |
| TV | 0.0458 | 0.001 | 32.809 | 0.000 | 0.043 | 0.049 |
| Radio | 0.1885 | 0.009 | 21.893 | 0.000 | 0.172 | 0.206 |
| Newspaper | -0.0010 | 0.006 | -0.177 | 0.860 | -0.013 | 0.011 |
此處Newspaper並不明顯，但實際上，在簡單回歸當中，Newspaper還是相當顯著的，但在多元線性回歸中卻不顯著。通過分析下面的相關系數我們可以發現，Newspaper和Radio的相關性是挺大的。這反映了，如果在某個市場通過Radio投放廣告，我們會習慣同時增加一些Newspaper廣告的投放，所以歸根到底其實是Radio的影響。sales，radio和newspaper三者的關系類似於shark_attack，temperature和ice-cream_sales。當引入真正重要的變量（radio或temperature）時，另一個與之相關的變量（newspaper或ice-cream_sales）就會不再顯著了。
| TV | Radio | Newspaper | Sales | |
|---|---|---|---|---|
| TV | 1.000000 | 0.054809 | 0.056648 | 0.782224 |
| Radio | 0.054809 | 1.000000 | 0.354104 | 0.576223 |
| Newspaper | 0.056648 | 0.354104 | 1.000000 | 0.228299 |
| Sales | 0.782224 | 0.576223 | 0.228299 | 1.000000 |
#計算結果
est = smf.ols('Sales ~ TV + Radio + Newspaper', advertising).fit()
est.summary().table[1]
#相關系數
advertising.corr()
2.2.2 一些重要問題
1.response和predictors是否相關？
F-test: $H_0 :β_1 =β_2 =···=β_p =0 $ , \(F=\frac{(\mathrm{TSS-RSS)}/p}{\mathrm{RSS}/(n-p-1)}\) 。當H~0~ 為真時，分子等於 \(\sigma^2\) ，而分母的期望就是等於 \(\sigma^2\) ，所以此時F接近1。否則，分子會大於 \(\sigma^2\) ，F也就大於1。但是否拒絕H~0~取決於n和p，當n足夠大，F-statistic會服從F分布，然后查表的結果，但最方便的方法是通過軟件看p值。
另一種F-test: $H_0 : β_{p−q+1} =β_{p−q+2} =...=β_p =0 $ , \(F=\frac{(\mathrm{RSS_0-RSS)}/q}{\mathrm{RSS}/(n-p-1)}\) 。判斷某predictors子集（q個元素）是否與response有關。RSS~0~是去掉q個predictors后回歸並計算RSS得到的。在極端情況下，即q=1時，F-statistics相當於(t-statistics)^2^ 。但通過t-statistics來判斷整個模型的顯著性是有很大風險的，因為每個t-statistics都是有一定置信水平。例如有5%錯過真值，當t檢驗次數一多，就很難保證全部檢驗都不會錯過真值。
2.變量選擇（C5詳細講）
Forward selection：從簡單回歸開始，選擇一個對減少RSS效果最好的變量，然后在一元回歸基礎上選擇下一個效果最大的變量……該策略問題：可能一開始引入一些后來變得不顯著的變量。
Backward selection：先減去p值最大的變量，直到所有變量的p值都很小。
Mixed selection：按照Forward開始，當加入新變量使得之前的某個變量p值不顯著時，去掉不顯著的……直到所有模型中的變量的p值很低，模型外的變量一加進來該變量的p值就很大。
3.模型擬合
RSE and R^2^，后者在多元回歸中等於 $\mathrm{Cor}(Y, \hat{Y})^2 $ 。就目前而言，看RSE比較好，能防止加入太多無關變量。另外，通過比較不同變量值情況下，模型是傾向於高估還是低估還能發現模型是否有忽略多變量交互作用的情況。
4.預測
目前有三個不確定性：
least squares plane \(\hat{Y}\) 和 population regression plane \(f(X)\) (注意這個不含 \(\epsilon\) )的差別用confidence interval衡量。例如95%的置信區間表示這些區間有95%是包含真值的。
模型bias
irreducible error
prediction intervals包含上面三種不確定性，例如95%的預測區間表示這些區間有95%包含某特定 x 對應的 y。自然，這個區間比置信區間大不少。
2.3 回歸模型的其他思考
2.3.1 Qualitative Predictors
two levels predictors:
下面是 \(x_i\) (indicator or dummy variable)為0/1編碼（也可以1/-1編碼，不同方式只影響解釋，不影響結果）
$β_0 $ 作為male的平均值，\(β_0 + β_1\)為female的均值，\(β_1\)就是均值差異。第二個式子沒有dummy variable成為baseline。
2.3.2 線性模型的拓展
標准線性回歸有兩個重要假設：additive和linear。前者表示不同的predictors對response的影響是獨立的，后者表示predictors一單位的變化對response的影響是constant的。
去除additive假設
additive假設限制了interaction效應，比如本章節數據中，TV和radio的廣告投入可能會相互促進sales。換句話說，向TV和radio分別投入0.5比只向TV或者radio單獨投入1效果要好。下面公式解除了additive假設：
\(Y = β_0 + β_1X_1 + β_2X_2 + β_3X_1X_2 + \epsilon\) ，其中 \(X_1X_2\) 稱為交互項
上式中，其中一個predictor的變化，對response的影響不再是固定的（如果另一個predictor也變化）。再次回歸的結果如下：
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 6.7502 | 0.248 | 27.233 | 0.000 | 6.261 | 7.239 |
| TV | 0.0191 | 0.002 | 12.699 | 0.000 | 0.016 | 0.022 |
| Radio | 0.0289 | 0.009 | 3.241 | 0.001 | 0.011 | 0.046 |
| TV:Radio | 0.0011 | 5.24e-05 | 20.727 | 0.000 | 0.001 | 0.001 |
另外，R^2^由原來的89.7%上升到了96.8%，這表示交互項解釋了之前模型RSS的(96.8 - 89.7)/(100 - 89.7) = 69 % 。根據hierarchical principle，如果引入了交互項，那么組成它的每個變量都要引入模型中，不管這些變量是否顯著。
在qualitative情況下，additive
non-additive
非線性（多項回歸）
\(Y = β_0 + β_1X_1 + β_2X_1^2 + \epsilon\) ，這是一個非線性函數，但仍是線性模型，所以在軟件里的回歸方式一樣。
2.3.3 Potential Problems
1.response-predictor 非線性相關
在多元回歸時，通過residual plot，即 \(e_i\) 和 \(\hat{y}_i\) 的散點圖，可以發現模型是否為non-linearity。正常情況下，散點圖是沒有discernible pattern的，比如U型。如果有，則數據有非線性聯系，可以對predictors進行非線性轉換，如開方，平方等。
2.error項的相關性
線性模型假設他們不相關，但如果相關，那么模型真實的置信水平會更低，結果是95%的confidence and prediction intervals實際上並沒有這么高，要把區間擴大才能達到95%。另外，也會高估parameter的顯著性。
相關的情況通常出現在時間序列數據，可以通過計算 e 和 time_id 的相關系數。而非時間序列中，比如height onto weight，如果一些人是來自相同家庭，或有相同的飲食習慣，或生活在相同的環境，這些因素都會破壞error的獨立性。
3.非固定方差的error項
如果 \(\mathrm{Var}(\epsilon_i) = σ^2\) 不成立，這會影響SE, confidence intervals 和 hypothesis tests。現實中經常不成立，比如在Residual plots中，rss的絕對值會隨response的增大而增大（漏斗狀）。解決這一問題可以通過對response進行轉換，如 \(log Y\) 或 $\sqrt{Y} $ 。另一種情況是，\(i\)th response是 \(n_i\) 個觀察值的均值，而這些觀察值和 \(σ^2\) 不相關，那么\(i\)th response的方差為 \(σ_i^2 = σ^2/n_i\) ，利用weighted least squares修正，此處乘上\(n_i\)
4.Outliers
outlier（下面20號）： \(y_i\) 的值和預測值相差很大很大。即使outlier不對回歸的結果造成大影響，它也會影響增大RSE和降低R^2^。Studentized residuals plots可以identify outlier，當Studentized residuals超過±3時為outliers。這是適用於收集數據時出錯產生的異常，但模型本身的缺陷，如缺失predictor也會出現outlier。
5.High Leverage Points
High Leverage Points（上面41號）比outliers嚴重，會大程度地影響擬合。通過計算leverage statistic來判斷，值介於1/n 和 1，所有observations的平均杠桿為(p + 1)/n，如果遠離該值，就很可能是High Leverage Points。
6.共線性
會降低回歸系數的准確性，導致SE(\(\hat{\beta}_j\))增加，t-statistic降低，從而錯誤地放棄一些重要變量。通過correlation matrix可以發現共線性，但如果是多重共線性，則要計算variance inflation factor (VIF) 。該因子最小值為1，表示完全沒有共線性，一般超過5或10才是有問題。解決方法是刪除一個或者利用共線變量構造一個新變量。
#計算三個變量age, rating, and limit各自的VIF
est_Age = smf.ols('Age ~ Rating + Limit', credit).fit()
est_Rating = smf.ols('Rating ~ Age + Limit', credit).fit()
est_Limit = smf.ols('Limit ~ Age + Rating', credit).fit()
print(1/(1-est_Age.rsquared))
print(1/(1-est_Rating.rsquared))
print(1/(1-est_Limit.rsquared))
3.5 線性回歸和KNN的比較
兩者的選擇一般由數據的線性程度來決定。然而，當predictor較多時，即便是非線性數據，線性回歸依然會比KNN好，這是由於curse of dimensionality。當數據是高維度時，non-parametric需要更多的樣本來擬合。
C3 Classification
3.1 總體介紹
本章節default數據是經過改造的，實際數據不會這么理想。
df['default2'] = df.default.factorize()[0] #[0]為標簽數組[1,0,0,0,1...]，[1]為unique values數組，Index(['No', 'Yes'], dtype='object')
df['student2'] = df.student.factorize()[0]
| default | student | balance | income | default2 | student2 | |
|---|---|---|---|---|---|---|
| 1 | No | No | 729.526495 | 44361.625074 | 0 | 0 |
| 2 | No | Yes | 817.180407 | 12106.134700 | 0 | 1 |
fig = plt.figure(figsize=(12,4)) # 整個結果框的大小
gs = mpl.gridspec.GridSpec(1, 4) # 對上面的框進行划分，這里是1行4列
ax1 = plt.subplot(gs[0,:-2]) # 占第一行直到倒數第二列（不包含）的位置
ax2 = plt.subplot(gs[0,-2])
ax3 = plt.subplot(gs[0,-1])
df_no = df[df.default2 == 0].sample(frac=0.15) # 從沒有違約records中抽樣
df_yes = df[df.default2 == 1]
df_ = df_no.append(df_yes)
ax1.scatter(df_[df_.default == 'Yes'].balance, df_[df_.default == 'Yes'].income, s=40, c='orange', marker='+',
linewidths=1)
ax1.scatter(df_[df_.default == 'No'].balance, df_[df_.default == 'No'].income, s=40, marker='o', linewidths='1',
edgecolors='lightblue', facecolors='white', alpha=.6)
ax1.set_ylim(ymin=0)
ax1.set_ylabel('Income')
ax1.set_xlim(xmin=-100)
ax1.set_xlabel('Balance')
c_palette = {'No':'lightblue', 'Yes':'orange'}
sns.boxplot('default', 'balance', data=df, orient='v', ax=ax2, palette=c_palette)
sns.boxplot('default', 'income', data=df, orient='v', ax=ax3, palette=c_palette)
gs.tight_layout(plt.gcf()) # 調整每張圖的距離
線性回歸最多適用於binary分類問題，但最好把prediction限制成probability，如下面的Logistic Regression。更多的類別在線性回歸中會出現無意義的ordering。
3.2 邏輯回歸
3.2.1模型
把回歸值限制在0～1。下面log叫做log odds
X和p(X)的關系是非直線的，其中一個變量的變化對另一個的影響取決於該變量的當前值。
該模型一般用maximum likelihood來估計，最小二乘法是它的一個特例。這個方法的目的是尋找最優參數組合，使得當default時， \(\hat{p}(x_i)\) 值接近1，反之接近0。公式表達如下，maximize這個likelihood function：
3.2.2模型計算和預測
y = df.default2.ravel() # ravel把series變為ndarray，不調用也可以
x = sm.add_constant(df.balance) # 給df.balance加上一列constant，值為1，返回DF
mult_x = sm.add_constant(df[['balance', 'income', 'student2']])
est = smf.Logit(y, x).fit()
est.summary2().tables[1] # 邏輯回歸要用summary2
est.predict([1,2000]) #0.5858 參數里的順序與上面x一樣，先constant
#sklearn
clf = skl_lm.LogisticRegression(solver='newton-cg')
x = df.balance.values.reshape(-1,1)
clf.fit(x,y) #可接受series
| Coef. | Std.Err. | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| const | -10.651331 | 0.361169 | -29.491287 | 3.723665e-191 | -11.359208 | -9.943453 |
| balance | 0.005499 | 0.000220 | 24.952404 | 2.010855e-137 | 0.005067 | 0.005931 |
上z可以當作t-statistics看待，從p值可以看出，系數是顯著的。
Confounding:
通過簡單回歸df.students2 model_1和多元回歸df[['balance', 'income', 'student2']model_2會發現前者中，學生的系數為正，后者為負。原因是，這個數據里，學生的default率更高（背后原因是因為學生傾向有更高的balance），所以考慮但因素時，model_1中系數為正。然而，當考慮上balance因數后會發現，在相同balance情況下，學生的default率更低（對balance的抵抗力更大），所以學生因素在model_2中成了負因子。現實中，如果不知道一個學生的balance信息，他的違約風險比其他人大。這再次說明單因素回歸的危險性。
邏輯回歸可用於多分類預測，但不常用。而當classes是well-separated時，邏輯回歸的參數估計很不穩定，下面模型能夠克服這個問題。
3.3 Linear Discriminant Analysis
當樣本量小，且predictors在各個classes中的分布接近normal時，此模型比邏輯回歸模型更穩定。
在不同的classes/labels情況下，找出屬於該class的observation的分布，然后根據貝葉斯理論，利用這些分布來推測predictors所屬的classes，哪個概率最大推測哪個。下面貝葉斯公式中，\(\pi_k\) 是class~k~在整個樣本中的比例；\(f_k (x)\) 是x在class~k~中的密度函數/分布，即已知\(Y=k\) ，\(X=x\)的概率是多少；分子不影響結果大小的順序，可不管。
\(\pi_k\) 比較好計算 \(\hat{π}_k = n_k / n\)，問題在於如何得出 \(f_k (x)\)。
3.3.1 LDA計算
LDA假設 \(f_k (x)\) 都是正太分布，且各分布同方差
當p=1，即一個predictor時
密度函數如下
通過樣本推斷出分布的均值和方差（假設同方差，所以用樣本求各自方差后加權平均），然后把x的值代入不同的 \(f_k (x)\) 再乘上相應的 \(\hat{π}_k\) ，最后選取結果最大的class作為預測結果。另一對class的結果相等就能解出兩者間的決策邊界。
當p>1
密度函數如下
這里 \(x\) 是n緯向量，重點是 \(\mu_k\) 和 \(\Sigma\) 。前者是一個向量(\(\mu_1, \mu_2,...\mu_n\))，后者是n x n的協方差矩陣。計算流程和p=1一樣。當n = 2，classes數量為3時，各 \(f_k (x)\) 分布的平面圖如下左圖。
上左圖的一個橢圓的范圍包含了該class 95%的概率。由於兩個predictors有相關性，如上右圖所示，所以分布成了橢圓，即便同方差。虛線划分是貝葉斯決策邊界，實線為LDA計算出的決策邊界。
X = df[['balance', 'income', 'student2']].as_matrix()
y = df.default2
lda = LinearDiscriminantAnalysis(solver='svd')
y_pred = lda.fit(X, y).predict(X)
3.3.2 QDA計算
去除各class同方差的限制，決策邊界不在是直線。
LDA和QDA的選擇在於bias-variance trade-off。前者只有kp個系數要估計，限制更多，variance更低，但bias更高，適合於樣本小和各class分布的方差差別不大的情況；后者有kp(p+1)/2個要估計，其他反之。
如果向LDA添加所有的二次項和交互項，它會和QDA一樣。
3.4 比較分類方法
LR和LDA兩者都是線性模型（式子轉換后很相似），區別在於系數的估計方法，前者用maximum likelihood，后者用來自正態分布的mean和variance。兩者的選擇在於數據是否來自於同方差的正態分布，數據是否相關不影響。
KNN是非參方法，比較適合非線形的決策邊界。
QDA處在上面兩類之間
3.5 Lab
df = pd.read_csv('Data/Smarket.csv', index_col=0)
# 了解下面兩個變量和其他所有變量的相關系數
df.corr()[["Volume","Today"]]
# factoriesDirection，可以用factorize()[0]，但具體編碼取決於數據的順序
df['Direction2'] = df.Direction.map(lambda x: 1 if (x == "Up") else 0)
#LR
# sklearn
feature = ["Lag1", "Lag2", "Lag3", "Lag4", "Lag5", "Volume"]
x = df[feature].values
y = df["Direction2"].values
log_reg = LogisticRegression()
log_reg.fit(x,y)
# statsmodels方便分析
y = df['Direction2']
x = sm.add_constant(df[feature])
est = smf.Logit(y, x).fit()
est.summary2().tables[1]
est.pred_table(threshold=0.5) # col: pre, row: label
#LDA
train_x = df[df.Year < 2005][['Lag1','Lag2']]
train_y = df[df.Year < 2005]['Direction']
test_x = df[df.Year == 2005][['Lag1','Lag2']]
test_y = df[df.Year == 2005]['Direction']
lda = LinearDiscriminantAnalysis()
pred = lda.fit(train_x, train_y).predict(test_x)
# 查看先驗概率
lda.priors_
# 比起confusion_matrix，這個更好用
print(classification_report(y_test, pred, digits=3))
""""
precision recall f1-score support
Down 0.500 0.315 0.387 111
Up 0.582 0.752 0.656 141
avg / total 0.546 0.560 0.538 252
""""
# 調整判斷概率
pred_p = lda.predict_proba(test_x4)
np.unique(pred_p[:,1]>0.7, return_counts=True)
C4 Resampling Methods
model assessment：評估模型表現的過程
model selection ：選擇一個模型合適的flexibility的過程
4.1 Cross-Validation
validation set：training和test set對半分。問題：highly variable，training set不足而高估錯誤率
Leave-one-out cross-validation (LOOCV) 能准確確定針對訓練數據的最優flexibility和錯誤率。問題：如果下面公式不成立（通常是非線性模型），則很耗時，高variance。
上面公式中h~i~為leverage statistic（在之前high leverage points有應用），它所在的分母的值介於1/n和1，反映了某觀測值的leverage對他所fit的模型的影響，即這個分母是用於抵消觀測值的leverage的。
k-Fold Cross-Validation（LOOCV的特例）相對省時，可大致確定最優flexibility，適中的variance。問題：不能較為准確地評估錯誤率。
4.2 The Bootstrap
用於測量某個對象的標准差，比如系數，且使用范圍廣。方法是從原始數據集中，通過重復抽樣產生多組與原始數據集相同大小的數據集，然后估計這些數據集的某個指標，比如均值\(\hat{\mu}_i\)，最后利用\(\hat{\mu}_i\)和他們的均值來求出SE(\(\hat{\mu}\))。
# validation set
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.33, random_state=None)
# LOOCV
loo = LeaveOneOut()
loo.get_n_splits(dataframe)
# k-Fold Cross-Validation
kf_10 = KFold(n_splits=10, shuffle=False, random_state=None)
# 上面兩個可放在下面的cv里
cross_val_score(model, X, y, cv=kf_10).mean()
# Bootstrap
X, y = resample(X, y) #可調節n_samples，默認和X大小一樣。得到的就是resample后的X和對應的y
C5 Linear Model Selection and Regularization
用其他擬合方式替代plain 最小二乘法能夠提高模型的prediction accuracy和model interpretability
Prediction Accuracy：如果response和predictors的關系是線性的，則最小二乘法估計是無偏的。再者，如果n遠大於p，則該估計還是低variance的。但如果n不是遠大於p，則會產生一定的variance。在這個階段，通過限制模型來提高部分bias降低variance。
Model Interpretability ：最小二乘法無法自動排除無關變量，一些輔助方法有：Subset Selection， Shrinkage/regularization，Dimension Reduction
5.1 Subset Selection
選擇方式：best subset 和 stepwise
前者通過對predictors所有的組合進行回歸后選出最優，2^p^種可能。問題：耗時和過度擬合。
后者分（介紹見2.2.2 一些重要問題的變量選擇）：forward, backward, and mix
選擇標准：
RSS和R^2^都傾向與包含所有predictors的模型來得到最低的training error。這並不適合test error的評估，由此本書引出下面兩類方法：
（1）用training errord的修正指標：
C~p~ 用於由最小二乘法得出的模型，該值等於 $ \frac{1}{n}(\mathrm{RSS} + 2d\hat{σ}^2)$ 。此式子中，d表示當前模型predictors的數量， \(\hat{\sigma}\) 由包含所有predictors的模型估計，\(2d\hat{σ}^2\) 用於抵消predictor數量的增加所帶來的error的減少，如果 \(\hat{\sigma}\) 是無偏估計，那么 \(C_p\) 也是test MSE的無偏估計。
AIC用於由極大似然得出的模型，該值等於 $ \frac{1}{n\hat{σ}^2}(\mathrm{RSS} + 2d\hat{σ}^2)$ 。此式子忽略一個定值。如果隨機項是正態分布的，那么極大似然和最小二乘法得出的結果是一樣的。可以看出AIC和前面的C~p~是成比例的。
BIC是從貝葉斯的觀點出發，公式與上面的類似，只是把\(2d\hat{σ}^2\)改為\(log(n)d\hat{σ}^2\)，但對添加predictors的懲罰更大，所以傾向於選擇predictors更少的模型。
Adjusted R^2^ 式子：$1-\frac{\mathrm{RSS}/(n − d − 1)}{\mathrm{TSS}/(n-1)} $ 。理論上而言，修正R^2^最高的模型只包含有用的變量。
（2）用validation：
one-standard-error rule：計算每個模型的test MSE，然后從test error在一個標准差內的模型中選test MSE最低的。
5.2 Shrinkage Methods
除了選predictors的子集來回歸，也可以通過regularizes來篩選變量。使用前要standardizing
Ridge Regression
通過調整ridge regression coefficient estimate $\hat{β}_λ^R $，即RSS里的系數，來最小化 \(\mathrm{RSS} + λ \sum_{j=1}^{p}\beta_j^2\) 。二項稱為shrinkage penalty，或寫作λ||β||~2~ 。式子中 \(\lambda ≥ 0\) ，用於調節懲罰的權重。當\(\lambda = 0\)時，\(||\hat{β}_λ^R||_2 = ||\hat{\beta}||_2\)。隨着\(\lambda\)的增大，各個predictors的系數$\hat{β}_λ^R $會進行調整（總體減少）。
Ridge Regression 優於最小二乘法在於它的bias-variance trade-off，意味着它可以在test set中表現得更好。
Lasso Regression 能夠去除一些不顯著的變量（另它的系數為0）。shrinkage penalty改為 \(λ \sum_{j=1}^{p}|\beta_j|\) 或λ||β||~1~
兩者適合的場景取決於有多少個無用的predictors，由於各個predictors的作用是未知的，所以一般用cross-validation來確定選擇。
補充：
上面兩個回歸等價於: minimize RSS subject to \(\sum_{j=1}^{p}\beta_j^2 ≤ s\) or \(\sum_{j=1}^{p}|\beta_j| ≤ s\) ，其中s對應於\(-\lambda\) ，即lambda越小，s越大，回歸的限制越小。從這個角度來看，兩個回歸在數學上相當於best subset selection的改良版，因為后者受限於\(\sum_{j=1}^{p}I(\beta_j \ne 0) ≤ s\) 。另外，通過這個角度，這兩個回歸就能在圖像中體現：
上圖部分說明了為什么Lasso能夠產生sparse解，因為交點會出現在坐標走上，而Ridge不會。
5.3 Dimension Reduction Methods
記得standardizing（以相同單位測量的除外）
對predictors進行降維（將p個變為m個），轉化后的predictors為 $Z_m = \sum_{j=1}^p \phi_{jm}X_j $ （每一個Z都是所有predictors的線性組合）其中m = 1,...,M，M<p，p為原數據中predictors的數量，\(\phi_{jm}\)為某定值。然后利用降維后的predictors進行回歸：
5.3.1 PCA
以一個包含兩個變量pop和ad的數據為例，轉換后得到的第一個Z~1~稱為第一主成分，其等式如下：
其中loadings φ~11~ = 0.839 和 φ~21~ = 0.544，它們是在 \(\phi_{11}^2 + \phi_{21}^2 = 1\)的條件下，使 $\mathrm{Var}(\phi_{11} × (pop − \bar{pop}) + \phi_{21} × (ad − \bar{ad})) $ 最大化得出的，為了盡量保存數據的variance。loadings越大，說明該predictors在此主成分占有的信息越大。而Z~2~是在與Z~1~正交，即不相關，的情況下得到的，其他Z如此類推。如果pop和ad的相關性很大，那么Z~1~（下圖綠線）能夠很好的保存這兩個變量的信息，而Z~2~（下圖虛線）由於與Z~1~正交，所以不能很好地使得所有的pop和ad都落在附近，所以保留的信息就少。在回歸時只回歸Z~1~就可以了。
其實ridge相當於連續性的PCR。用於回歸時，M的選擇可以根據cross-validation；用於數據分析時，M可以根據以PVE和M為坐標的圖中找拐點選擇，或者開始看第一主成分，如果有感興趣的特征，再加上第二主成分。
使用：PCA通常只在有很多predictors的情況下才表現的更好，機器正常能處理的數據一般不需要用。另外，由於降維所得到的主成分是unsupervised，所以無法保證它們是最適合回歸的。
5.3.2 PLS
PLS在降維得出Z~1~的時候，把更多的權重賦給與response更相關的變量。而Z~2~則從“對Z~1~的每一個變量進行回歸后提取得到的residuals”中采用Z~1~同樣的方法獲得。然而此方法的實際表現並不比PCA好。
5.4 思考高維數據
在數據不足的情況下：最小二乘法沒有最優解，即便有，也會過度擬合，所以MSE和R~2~失效。但其他修正指標在高維數據中估計所需的\(\hat{\sigma}\)也有問題。這些問題同樣會出現在傳統的線形模型中。這時需要的是獨立測試或者cross-validation。而在回歸方面就需要本章提到的一些方法，如ridge、lasso、stepwise selection、PCA等來防治過度擬合。在高維數據中，正確地判斷\(\lambda\)能夠在提高模型的表現，但除非附加特征與response正相關，否則測試誤差始終會隨着問題predictors的增加而增加。
在高維數據中，多重共線性會很嚴重的，這使得我們無法判斷哪個變量才是真正有用的。即便使用stepwise所得到的模型，也只是眾多可行模型中的一個，或許還只在某個特定數據集可行。
5.5 Lab
# 下面三行作閱讀理解，用sklearn的onehot更方便
dummies = pd.get_dummies(df[['League', 'Division', 'NewLeague']])
X_ = df.drop(['Salary', 'League', 'Division', 'NewLeague'], axis=1).astype('float64')
X = pd.concat([X_, dummies[['League_N', 'Division_W', 'NewLeague_N']]], axis=1)
# Ridge（和lasso類似）如果結果想和書本的幾乎一樣，需要安裝一個新模塊python-glmnet
ridge = Ridge(alpha) # alpha是懲罰系數
ridge.fit(scaler.transform(X_train), y_train) # scale相當於簡單化的StandardScaler，沒有fit和transform
pred = ridge.predict(scaler.transform(X_test))
mean_squared_error(y_test, pred)
pd.Series(ridge2.coef_.flatten(), index=X.columns)
# cross-validation訓練ridge
ridgecv = RidgeCV(alphas, scoring='neg_mean_squared_error') #alphas是一個lambda值數組
ridgecv.fit(scale(X_train), y_train)
ridgecv.alpha_ # 得出最優alpha
# PCA
pca = PCA()
X_reduced = pca.fit_transform(scale(X)) # 得到主成分，默認多少predictors就多少個
pca.components_ # 查看loadings
C6 非線性模型
本章介紹的方法能盡量保留interpretability，算是向非線形模型方向前進的一小步。
6.1模型介紹
1～4可以看作是簡單線性的擴展。5則是多元回歸
6.1.1 Basis Functions
多項回歸
一般來說，d很少大於4，因為會過度flexible，產生奇怪的形狀，尤其是在predictors的邊界。左下圖表示多項回歸的結果。虛線代表95%的置信區間。從圖中可以看出，wage可以分以250為邊界分為兩組，由此運用多項邏輯回歸。結果如有下圖，由於wage>250在總體里站很小的比重，所以95%的置信區間會很寬。
Step functions：把numerical變量變為categorical后回歸。
Basis Functions ：就是將上面兩種functions的式子寫到一起。
上面的 $b_K(·) $ 既可以是 $b_j(x_i) = x^j_i $ （多項）也可以是 $b_j(x_i) = I(c_j ≤ x_i < c_{j+1}) $ （step）
6.1.2 Regression Splines
Piecewise Polynomials：將變量進行分區，在不同區域運用相同形式的多項回歸，自由度為所有多項回歸系數的數量。由於過於flexible，各分區之間會出現斷點。
Continuous Piecewise：通過添加限制是函數變成連續。下面以cubic spline with K個斷點結合Basis Functions的模型為例：
有n個分區就加上n-1個系數，其中 \(b(x_i) = h(x_i, ξ)\) 。當 \(x_i > ξ\) 等於 \((x-ξ)^3\) ，否則等於0。這樣回歸的結果就能在二次求導內保持連續的效果。
添加additional boundary constraints 成為natural cubic spline，即在x的兩個邊界部分限制為線性來減少這兩部分的variance。
斷點的選擇：在數據變化劇烈的區域多設斷點，雖然實踐中通常按數據比例進行划分，如25%，50%和75%，或者cross-validation。
與多項回歸的比較：當spline的自由度和多項式的degree水平一樣時，前者通常更好。
6.1.3 Smoothing Splines
最小化目標函數：
第二項是懲罰項，積分的是整個函數在一階導數的總變化。如果g很平滑（直線），那么懲罰項就小，反之。（二階導數就是測量函數變化的劇烈程度的，如果值為0，說明函數在該點平滑/直線）當懲罰系數無窮大時，g就會變成直線，反之g會interpolate。
其實這個函數類似shrunken 版本的piecewise natural cubic polynomial，懲罰系數控制shrinkage 程度，另外所有 unique values of x1,...,xn 都是斷點。
懲罰系數的選擇：用LOOCV
λ increases from 0 to ∞, effective degrees of freedom($df_λ $) decrease from n to 2
6.1.4 Local Regression
計算某個x的target只用該點附近的x來推定，這個附近用比例s確定，可設置一些變量為global，有些為local。然而當變量多於3個時，變現並不好，因為附近沒有那么多observation用於訓練。類似KNN的問題。
6.1.5 Generalized Additive Models
additive model: 一個自然的方法是把所有 \(β_jx_{ij}\) 換成 $f_j(x_{ij}) $ ，得到 $y_i = β_0 +f_1(x_{i1})+f_2(x_{i2})+···+f_p(x_{ip})+ \epsilon_i $ 然后每個 \(f_i\) 可以是上面1～4所介紹的函數。
由於additive，對於每個變量對response的影響還是可以解釋的。這個模型為參數和非參數方法提供了一個折中。
在分類0/1問題中，如果某個 $f_j(x_{ij}) $ 中x是categorical，且里面的某個類別連一個1都沒有，則要去掉。
C7 Tree-Based Methods
盡管decision tree比不上C5和C6的方法，但bagging和random forests在prediction方面可以有很大提升。
7.1 決策樹基礎
回歸樹最小化RSS，目標函數：
其中J代表terminal node的集合，R~j~表示第j個node，\(\hat{y}_{R_j}\) 表示該node的預測值。
通常會設置一個 \(\alpha\) 來作懲罰系數，防止過度擬合。
分類樹最小化分類錯誤率，Gini或entropy：
其中 \(\hat{p}_{mk}\) 表示第m個terminal node中kth class的比例。由於分類錯誤率對樹的成長不敏感，所以通常用后面兩種，但對於樹的判讀正確率，E表現更加理想。
樹的優缺點：易於解釋和運用，但對樣本數據的變化很敏感，且效果不比之前介紹的方法好。
7.2 Bagging, Random Forests, Boosting
樹夠多的話，前兩個不會overfitting
Bagging
bootstrap多組數據集來建樹，對於回歸，取樹的均值；對於分類，最多投票的class。
Out-of-Bag (OOB) Error Estimation
變量重要性測量：某變量的划分對RSS/Gini or entropy減少的總量
Random Forests
划分時考慮隨機選取的部分樣本，為了防止生長成相同的樹。其他和bagging一樣。
Boosting
每棵樹在前一個棵樹的訓練數據的修改版response的residual部分上進行訓練。樹數量，學習速度，樹的深度等需要cross-validation決定
C8 SVM
8.1 Maximal Margin Classifier
Hyperplane：在p維空間里，a flat affine（不過原點） subspace of dimension p − 1。下面式子定義一個超平面：
任何 $X = (X_1,X_2,...,X_p)^T $ 表示超平面上的一個點。如果式子變為大於小於，那么表示該點處於超平面的其中一邊。
support vectors：影響邊界的observation，此方法中是在margin邊界上的observation
最大化邊界M，數學表達如下：
限制是為了每個observation都在正確的一邊，而且距離至少為M。
由於此方法在分界上過度敏感，且有時數據不是完全可分的，所以引出下面方法。
8.2 Support Vector Classifiers(soft margin classifier)
允許observation落在錯誤的一側（超過margin）
上面 \(\epsilon_i\) 是slack variables，等於0是正確划分，大於0小於1是在margin里錯誤的一邊，大於1就是在margin外錯誤的一邊。C代表對這些錯誤的總量。
此方法中的support vectors是邊界上以及錯誤地超越邊界的observation。
然而此方法的分界依然是線性的，為了解決不同class數據非線性分堆的問題，需要下面方法。
8.3 Support Vector Machines
通過enlarging the feature space，即用平方、立方或更高階的多項式實現。利用kernel trick可以簡化這些多項式：
兩觀測值的內積： \(⟨x_i,x_{i^{'}}⟩ = \sum_{j=1}^p x_{ij}x_{i^{'}j}\)
linear support vector classifier 可以簡化為：\(f(x) = \beta_0 + \sum_{i=1}^n a_i⟨x,x_i⟩\) 其中 \(\alpha\) 是每個observation都要估計的系數，但如果該觀測值不是support vector，則該系數為0。
為了估計這個函數，需要計算內積，把上面的內積歸納為： \(K(x_i,x_{i^{'}})\) ，K成為kernel，一個量化兩個observation相似程度的函數。上面的kernel是線性的，也可以改用其他形式的kernel。當support vector classifier 采用其他非線形kernel時，稱為SVM。
radial kernel（rbf）根據附近的observation划分，需要調節gamma，越大越flexible
多個分類以上的SVM：One-Versus-One 和 One-Versus-All
SVM和LR類似，因為8.2中的式子可以改寫如下：
第一部分loss function（hinge loss）和LR的很類似，只是對於非support vector，SVM的loss為0，LR中遠離分界線的接近0。總的來說，SVM在 well-separated classes時表現更好，重疊多時后者更好。
8.4 Lab
svc2 = SVC(C=0.1, kernel='linear')
# 用cross-validation
tuned_parameters = [{'C': [0.001, 0.01, 0.1, 1, 5, 10, 100]}]
clf = GridSearchCV(SVC(kernel='linear'), tuned_parameters, cv=10, scoring='accuracy', return_train_score=True)
clf.fit(X, y)
clf.cv_results_
clf.best_params_
C9 Unsupervised Learning
PCA
9.1 Clustering Methods
K-Means
目標：同類的variation盡量小
C~k~表示第k個cluster，W(.)表示衡量variation的函數，具體實現有很多，例如L2。
K-Means的一個問題是需要預先設置K的個數，而hierarchical clustering可以后來再設K。
hierarchical clustering
相似度根據兩個observations交匯點的高度，越高越不相似。下左圖，2，8，5和7與9的交匯高度是一樣的。
分類就看從哪個高度水平切開，得出多少樹枝就多少類。
這個方法的問題出現在這么一個場景：在一組數據中，如果分成兩組的最優解是男/女，而分成三組的最優解是國籍（三個國家）。這個方法分三組的最優方法會是在男/女的基礎上把男或女的組盡心拆分。在這種情況下，它分類的“准確性”比KMeans差。
實現的過程：一開始n個clusters，然后組合最接近的兩個observation得到n-1個cluster，然后不斷重復。判斷兩個cluster接近程度的方法有：complete，average，single和centroid，前兩者最受歡迎，因結果樹更加平衡。
補充：
Dissimilarity Measure 的選擇：距離或者關系。例如按距離來分，會把購買量類似的顧客分為一類；按關系來分，會把購買種類類似的顧客分為一類。
對於數據的scale取決於應用的場景。
上面兩個方法的問題：會對一些outliers強行分類。
參考：
Introduction to Statistical Learning