陰謀仍是悲劇？- 基於機器學習假設檢驗視角，看泰坦尼克號事件

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之。

筆者這裏想說的是，不管什麼時代，何時，總會有噴子無腦地拋出陰謀論的說法，好像全部事情背後總有一個看不見的手在操縱之類的，可是，歷史老是在上演各類各樣的大事件，全部事件的發生與否從單個事件來看也許是偶然離奇事件，但放到具體歷史背景下，合理性是無可置疑的。

相反，當某個巨大災難發生時，當民族遇到危機，當國家遭到危難時，站出來的老是看起來沉默寡言的英雄，正是這些英雄的存在，慰藉了咱們的心靈，給了咱們溫暖與前進的力量，向英雄致敬。