陰謀還是悲劇？- 基於機器學習假設檢驗視角，看泰坦尼克號事件
- 2019 年 10 月 3 日
- 筆記
1. 引言
0x1：故事背景
泰坦尼克號（RMS Titanic），又譯作鐵達尼號，是英國白星航運公司下轄的一艘奧林匹克級郵輪，排水量46000噸，於1909年3月31日在北愛爾蘭貝爾法斯特港的哈蘭德與沃爾夫造船廠動工建造，1911年5月31日下水，1912年4月2日完工試航。
泰坦尼克號是當時世界上體積最龐大、內部設施最豪華的客運輪船，有“永不沉沒”的美譽 。然而不幸的是，在它的處女航中，泰坦尼克號便遭厄運——它從英國南安普敦出發，途經法國瑟堡-奧克特維爾以及愛爾蘭科夫(Cobh)，駛向美國紐約。1912年4月14日23時40分左右，泰坦尼克號與一座冰山相撞，造成右舷船艏至船中部破裂，五間水密艙進水。次日凌晨2時20分左右，泰坦尼克船體斷裂成兩截後沉入大西洋底3700米處。2224名船員及乘客中，逾1500人喪生，其中僅333具罹難者遺體被尋回。泰坦尼克號沉沒事故為和平時期死傷人數最為慘重的一次海難，其殘骸直至1985年才被再度發現，目前受到聯合國教育、科學及文化組織的保護。
這場聳人聽聞的悲劇震驚了國際社會，並為船舶制定了更好的安全規定。造成海難失事的原因之一是乘客和機組人員沒有足夠的救生艇。儘管倖存下沉有一些運氣因素，但有些人比其他人更容易生存，例如婦女，兒童和上流社會。
0x2：陰謀還是悲劇
在事件發生後，當時有很多陰謀論者懷疑泰坦尼克號事件是一個世紀大陰謀，是航運公司的一次騙保行為，一度爭論了很久。
當然，筆者並不是歷史學家，所以我們這次不從證據學的角度來分析這次事件。相反，我們嘗試運用統計學的思維，從數據挖掘與分析層面來看一下，這次事故是真實的可能性有多少。
基本假設是這樣的，我們首先建立一個二元對立假設，這是本文的最大的基礎：
- 如果是一次騙保行為，則航運公司的目標是將 船給摧毀，除此之外不會有任何其他的行為，那麼船長也不應該表現出任何的紳士和道德楷模風度，或者說，即使是作秀，表現出的道德程度也不應該很強。
- 如果是一個真實事故，則所有人（包括船長）表現出的行為都是按照人的第一反應作出的，則按照事後倖存者的採訪記錄，副船長愛德華·約翰·史密斯表現出了驚人的紳士風度和道德楷模風範，他說出了那句著名的話，『 lady and kid first！』，而自己則堅守到最後，非常令人感動。
接下來，我們將通過對數據進行建模，從而看某些特徵指標（婦女、小孩）和生還之間是否存在強相關關係，以此推斷出上述二元對立的結果。
Relevant Link:
https://baike.baidu.com/item/%E6%B3%B0%E5%9D%A6%E5%B0%BC%E5%85%8B%E5%8F%B7/5677
2. 數據建模
0x1：數據集
原始數據集可以從kaggle這裡下載得到，kaggle的原題是做predict的，所以將數據集分為了train/test set，但是筆者這裡並不做預測任務，因此將其合併起來作為一個統一的有標數據集，百度網盤鏈接在這裡。
鏈接: https://pan.baidu.com/s/18Ushn_Oms8Q0RLw6TGFMnQ 提取碼: u5gt
0x2：數據集欄位探查
各個特徵欄位意義如下：
0x3：數據缺失值補充
在很多統計項目中，由於數據採集手段上的限制，缺失值是非常常見的。因此，這裡我們先進行缺失值處理。
# -*- coding: utf-8 -*- import pandas as pd import numpy as np if __name__ == '__main__': dataset_file = "full_dataset.csv" dataset = pd.read_csv(dataset_file) print(dataset.info())
可以看到，【Sex，Age，Ticket，Fare，Cabin，Embarked】這幾個欄位都存在缺失值。
缺失值處理的思路有以下兩種：
- 1. 默認填充值的範圍[(mean – std) ，(mean + std)]。
- 2. 將缺失的Age當做label，將其他列的屬性當做特徵，通過已有的Age的記錄訓練模型，來預測缺失的Age值。
這裡採取方法1，以Age為例：
for dataset in full_data: age_avg = dataset['Age'].mean() age_std = dataset['Age'].std() age_null_count = dataset['Age'].isnull().sum() age_default_list = np.random.randint(low=age_avg-age_std,high=age_avg+age_std,size=age_null_count) dataset['Age'][np.isnan(dataset['Age'])] = age_default_list dataset['Age'] = dataset['Age'].astype(int)
其他欄位以此類推，這裡不再贅述，文末會給出完整程式碼。
0x4：離散特徵數值化處理
對於有限離散的字元串特徵，還需要進行映射處理或者one-hot處理。
def Passenger_Embarked(x): Embarked = {'S':0, 'C':1, 'Q':2} return Embarked[x]
0x5：特徵二次處理 – 專家經驗特徵
除了原始數據集的基礎特徵之外，還可以基於專家經驗進行組合特徵，當然這可以通過AutoEncoder進行自動完成，之所以不這麼做的原因是因為AE的內部機制不容易對外暴露，可解釋性工作會變得很複雜，因為我們這裡使用人工特徵組合，並使用XGB進行模型擬合。
1. 家族名特徵
在所有特徵中注意到有一個Name欄位，這是一個字元串特徵，字元串本身是無法輸入模型的。顯然，我們不能直接將Name進行ascii向量化，因為名字本身是沒有任何意義的，我們需要將其groupby處理後得到某種統計意義上的特徵。
我們發現Name的title name是存在類別的關係的。於是可以對Name進行提取出稱呼這一類別title name。
import re def get_title_name(name): title_s = re.search(' ([A-Za-z]+).', name) if title_s: return title_s.group(1) return "" for dataset in full_data: dataset['TitleName'] = dataset['Name'].apply(get_title_name) print(pd.crosstab(train['TitleName'],train['Sex'])) ###### out Sex female male TitleName Capt 0 1 Col 0 2 Countess 1 0 Don 0 1 Dr 1 6 Jonkheer 0 1 Lady 1 0 Major 0 2 Master 0 40 Miss 182 0 Mlle 2 0 Mme 1 0 Mr 0 517 Mrs 125 0 Ms 1 0 Rev 0 6 Sir 0 1 ####可以看出不同的titlename中男女還是有區別的。進一步探索titlename對Survived的影響。 ####看出上面的離散取值範圍還是比較多，所以可以將較少的幾類歸為一個類別。 train['TitleName'] = train['TitleName'].replace('Mme', 'Mrs') train['TitleName'] = train['TitleName'].replace('Mlle', 'Miss') train['TitleName'] = train['TitleName'].replace('Ms', 'Miss') train['TitleName'] = train['TitleName'].replace(['Lady', 'Countess','Capt', 'Col', 'Don', 'Dr', 'Major', 'Rev', 'Sir', 'Jonkheer', 'Dona'], 'Other') print (train[['TitleName', 'Survived']].groupby('TitleName', as_index=False).mean()) ###### out1 TitleName Survived 0 Master 0.575000 1 Miss 0.702703 2 Mr 0.156673 3 Mrs 0.793651 4 Other 0.347826
可以看出TitleName對存活率還是有影響差異的，這就說明了TitleName中蘊含了資訊，有資訊就可以為模型所用。
從資訊增益的角度出發，最後將TitleName總共為了5個類別：Mrs,Miss,Master,Mr,Other。
2. 家庭大小
SibSp和Parch分別為同船的兄弟姐妹和父母子女數，離散數據，沒有缺失值。於是可以根據該人的家庭情況組合出不同的特徵。
並且，根據FamilySize大小增加新特徵IsAlone。
0x6：模型訓練
# -*- coding: utf-8 -*- import pandas as pd import numpy as np import re from sklearn.preprocessing import LabelBinarizer from sklearn.ensemble import RandomForestRegressor import xgboost as xgb from xgboost import plot_importance import matplotlib.pyplot as plt def Passenger_sex(x): sex = {'female': 0, 'male': 1} return sex[x] def Passenger_Embarked(x): Embarked = {'S': 0, 'C': 1, 'Q': 2} return Embarked[x] def get_title_name(name): title_s = re.search(' ([A-Za-z]+).', name) if title_s: return title_s.group(1) return "" def Passenger_TitleName(x): TitleName = {'Mr':0,'Miss':1, 'Mrs':2, 'Master':3, 'Other':4} return TitleName[x] def data_feature_engineering(full_data, age_default_avg=True, one_hot=True): """ :param full_data:全部數據集，包括train,test :param age_default_avg: age默認填充方式，是否使用平均值進行填充 :param one_hot: Embarked字元處理是否是one_hot編碼還是映射處理 :return:處理好的數據集 """ for dataset in full_data: # Pclass,Parch,SibSp不需要處理 # Sex用 0，1做映射，{'female':0, 'male':1} dataset['Sex'] = dataset['Sex'].map(Passenger_sex).astype(int) # SibSp和Parch均沒有空值，組合特徵FamilySize dataset['FamilySize'] = dataset['SibSp'] + dataset['Parch'] + 1 # 根據FamilySize大小增加新特徵IsAlone dataset['IsAlone'] = 0 isAlone_mask = dataset['FamilySize'] == 1 dataset.loc[isAlone_mask, 'IsAlone'] = 1 # Fare數據存在空值，先填充 fare_median = dataset['Fare'].median() dataset['Fare'] = dataset['Fare'].fillna(fare_median) # Embarked存在空值，眾數填充，映射處理或者one-hot 編碼 dataset['Embarked'] = dataset['Embarked'].fillna('S') dataset['Embarked'] = dataset['Embarked'].astype(str) if one_hot: # 因為OneHotEncoder只能編碼數值型，所以此處使用LabelBinarizer進行獨熱編碼 Embarked_arr = LabelBinarizer().fit_transform(dataset['Embarked']) dataset['Embarked_0'] = Embarked_arr[:, 0] dataset['Embarked_1'] = Embarked_arr[:, 1] dataset['Embarked_2'] = Embarked_arr[:, 2] dataset.drop('Embarked', axis=1, inplace=True) else: # 字元映射處理 dataset['Embarked'] = dataset['Embarked'].map(Passenger_Embarked).astype(int) # Name選取稱呼Title_name，{'Mr':0,'Miss':1, 'Mrs':2, 'Master':3, 'Other':4} dataset['TitleName'] = dataset['Name'].apply(get_title_name) dataset['TitleName'] = dataset['TitleName'].replace('Mme','Mrs') dataset['TitleName'] = dataset['TitleName'].replace('Mlle', 'Miss') dataset['TitleName'] = dataset['TitleName'].replace('Ms', 'Miss') dataset['TitleName'] = dataset['TitleName'].replace(['Lady', 'Countess', 'Capt', 'Col', 'Don', 'Dr', 'Major', 'Rev', 'Sir', 'Jonkheer', 'Dona'], 'Other') dataset['TitleName'] = dataset['TitleName'].map(Passenger_TitleName).astype(int) # age 缺失值，可以使用隨機值填充，也可以用RF預測 if age_default_avg: # 缺失值使用avg處理 age_avg = dataset['Age'].mean() age_std = dataset['Age'].std() age_null_count = dataset['Age'].isnull().sum() age_default_list = np.random.randint(low=age_avg-age_std, high=age_avg+age_std, size=age_null_count) dataset.loc[np.isnan(dataset['Age']),'Age'] = age_default_list dataset['Age'] = dataset['Age'].astype(int) else: # 將age作為label，使用RF的age # 特徵為TitleName,Sex,pclass,SibSP,Parch,IsAlone,CategoricalFare,FamileSize,Embarked feature_list = ['TitleName', 'Sex', 'Pclass', 'SibSp', 'Parch', 'IsAlone', 'CategoricalFare', 'FamilySize', 'Embarked', 'Age'] if one_hot: feature_list.append('Embarked_0') feature_list.append('Embarked_1') feature_list.append('Embarked_2') feature_list.remove('Embarked') Age_data = dataset.loc[:, feature_list] un_Age_mask = np.isnan(Age_data['Age']) Age_train = Age_data[~un_Age_mask] #要訓練的Age feature_list.remove('Age') rf0 = RandomForestRegressor(n_estimators=60, oob_score=True, min_samples_split=10, min_samples_leaf=2, max_depth=7, random_state=10) rf0.fit(Age_train[feature_list], Age_train['Age']) def set_default_age(age): if np.isnan(age['Age']): data_x = np.array(age.loc[feature_list]).reshape(1,-1) age_v = np.round(rf0.predict(data_x))[0] print(age_v) return age_v return age['Age'] dataset['Age'] = dataset.apply(set_default_age, axis=1) return full_data # 特徵選擇 def data_feature_select(full_data): """ :param full_data:全部數據集 :return: """ drop_list = ['PassengerId', 'Name', 'Ticket', 'Cabin'] for data_set in full_data: # remove meaningless feature data_set.drop(drop_list, axis=1, inplace=True) train_y = np.array(full_data[0]['Survived']) train = full_data[0].drop('Survived', axis=1, inplace=False) return train, train_y def modelfit(alg,dtrain_x,dtrain_y,useTrainCV=True,cv_flods=5,early_stopping_rounds=50): """ :param alg: 初始模型 :param dtrain_x:訓練數據X :param dtrain_y:訓練數據y（label） :param useTrainCV: 是否使用cv函數來確定最佳n_estimators :param cv_flods:交叉驗證的cv數 :param early_stopping_rounds:在該數迭代次數之前，eval_metric都沒有提升的話則停止 """ if useTrainCV: xgb_param = alg.get_xgb_params() xgtrain = xgb.DMatrix(dtrain_x,dtrain_y) print(alg.get_params()['n_estimators']) cv_result = xgb.cv(xgb_param,xgtrain,num_boost_round = alg.get_params()['n_estimators'], nfold = cv_flods, metrics = 'auc', early_stopping_rounds=early_stopping_rounds) print(cv_result) print(cv_result.shape[0]) alg.set_params(n_estimators=cv_result.shape[0]) # train data alg.fit(train_X,train_y,eval_metric='auc') #predict train data train_y_pre = alg.predict(train_X) print ("nModel Report") print ("Accuracy : %.4g" % metrics.accuracy_score(train_y, train_y_pre)) feat_imp = pd.Series(alg.get_booster().get_fscore()).sort_values(ascending=False) feat_imp.plot(kind = 'bar',title='Feature Importance') plt.ylabel('Feature Importance Score') plt.show() print(alg) if __name__ == '__main__': dataset_file = "./full_dataset.csv" full_data = pd.read_csv(dataset_file) full_data = [full_data] # feature engineering full_data = data_feature_engineering(full_data, age_default_avg=True, one_hot=False) # select feature train_X, train_y = data_feature_select(full_data) print train_X print train_y # train the model # 1. set booster param param = {'max_depth': 7, 'eta': 1, 'silent': 1, 'objective': 'reg:linear'} param['nthread'] = 4 param['seed'] = 100 param['eval_metric'] = 'auc' num_round = 10 dtrain = xgb.DMatrix(data=train_X, label=train_y, missing=-999.0) bst_xgboost_model = xgb.train(param, dtrain, num_round) # print the feature importance table ### plot feature importance fig, ax = plt.subplots(figsize=(15, 15)) plot_importance(bst_xgboost_model, height = 0.5, ax=ax, max_num_features=64) plt.show()
3. 結果分析
Xgboost訓練後得到的特徵重要性評估表如下：
從結果中可以得到幾點資訊：
- Age年齡確實是和是否生還強相關的，佔據了主導地位
- 緊隨年齡之後，排名第二的生還因素是Fare（船票價格），這從一定程度上說明，不管是出於富人逃生意識更強，富人有更多機會上到頂層夾板玩了透氣，還是說逃生的時候有意識地偏向了富人，財富在其中的作用一定程度上是存在的
- 排名第三的是Sex（性別），這也副船長的lady fisrt的說法也是一致的
- 往後的幾個特徵，從各個維度都說明了生還幾率和財富有關
從數據和統計層面可以解讀出很多資訊，更深層次的猜測，讀者朋友可以參閱下面一些鏈接：
https://baijiahao.baidu.com/s?id=1623250627137421869&wfr=spider&for=pc https://baijiahao.baidu.com/s?id=1621654492865556041&wfr=spider&for=pc https://baijiahao.baidu.com/s?id=1597619726579770526&wfr=spider&for=pc http://www.ufo-1.cn/article/201604/896.html
同樣的分析方法，也可以用於二戰珍珠港事件以及911事件的分析，讀者朋友可以Google之。
筆者這裡想說的是，無論什麼時代，什麼時候，總會有噴子無腦地拋出陰謀論的說法，好像所有事情背後總有一個看不見的手在操縱之類的，但是，歷史總是在上演各種各樣的大事件，所有事件的發生與否從單個事件來看也許是偶然離奇事件，但放到具體歷史背景下，合理性是不容置疑的。
相反，當某個巨大災難發生時，當民族遇到危機，當國家遭到危難時，站出來的總是看起來沉默寡言的英雄，正是這些英雄的存在，慰藉了我們的心靈，給了我們溫暖與前進的力量，向英雄致敬。
