利用 gplearn 進行特徵工程

07-28

利用 gplearn 進行特徵工程

來自專欄猴子聊數據分析4 人贊了文章

featuretools有專門的 transfrom 類的特徵處理方式，但是我在實際使用中沒有使用，因為featuretool的思路是凡是能夠進行特徵變換的特徵都要應用一遍，所以應用的模式基本上是生成大量的性能不強的特徵，下一步必須進行嚴格的。

所以今天嘗試用gplearn進行一下特徵工程，本人能力有限，如有錯誤還請指正，也歡迎交流。

gplearn是Python內最成熟的符號回歸演算法實現，作為一種一種監督學習方法，符號回歸（symbolic regression）試圖發現某種隱藏的數學公式，以此利用特徵變數預測目標變數。

符號回歸的具體實現方式是遺傳演算法（genetic algorithm）。一開始，一群天真而未經歷選擇的公式會被隨機生成。此後的每一代中，最「合適」的公式們將被選中。隨著迭代次數的增長，它們不斷繁殖、變異、進化，從而不斷逼近數據分布的真相。

漢語介紹
官方文檔

當然我這裡不是想嘗試符號回歸的機器學習方式，而是符號回歸的機器學習方式提供了另外一種生成特徵的思路，在庫中通過Symbolic Transformer 類進行有關特徵工程相關操作。

本文還是以 Home Credit Default Risk這個項目進行說明。

在手工特徵工程中，我們會對兩個或者多個特徵進行一些加減乘除的操作，來生成一些特徵，希望能夠生成一些

例如在start-here-a-gentle-introduction的這篇kernel中，

根據領域的先驗知識，對金額特徵, 日期特徵進行比值操作生成一些特徵，這些特徵經常能夠提升驗證集和測試集的分數，在模型中也有很高的重要程度。

app_test_domain[DAYS_EMPLOYED_PERCENT] = app_test_domain[DAYS_EMPLOYED] / app_test_domain[DAYS_BIRTH]

但畢竟領域的先驗知識是有限的，我們有時候也沒有足夠的先驗知識，同時自己進行特徵組合也經常費時費力，自己生成的大量特徵不是都有足夠的重要性，還需要進一步篩選。

gplearn這個庫提供了解決的思路：

- 隨機化生成大量特徵組合方式，解決了沒有先驗知識，手工生成特徵費時費力的問題

- 通過遺傳演算法進行特徵組合的迭代，而且這種迭代是有監督的迭代，存留的特徵和label相關性是也來越高的，大量低相關特徵組合會在迭代中被淘汰掉，用決策樹演算法做個類比的話，我們自己組合特徵然後篩選，好比是後剪枝過程，遺傳演算法進行的則是預剪枝的方式。

變異的方式有興趣的可以看看文檔變異，具體遺傳演算法的方式我只是看看思路，並沒有仔細推導。

官方給出的例子是用 sklearn中自帶的波士頓房價數據來進行特徵組合，但是只是給出了一個新生成的特徵缺失提升了測試集的效果的結論。

所以，嘗試還是用官方的數據，探索一下新生成的特徵到底表現如何。

從一下幾個方面考察：

特徵和label的相關性
特徵之間的共線性
特徵在模型中的重要程度
特徵對模型性能的提升

官方例子中，使用的是Ridge回歸，使用 r2_score 評價回歸效果。在生成10個特徵後，0.7591提升至0.8417。

我在測試中使用mse損失作為模型評價標準。

1 控制實驗

分析一下原有的特徵情況

def data_prepare(): boston = load_boston() boston_feature = pd.DataFrame(boston.data, columns=boston.feature_names) boston_label = pd.Series(boston.target).to_frame("TARGET") boston = pd.concat([boston_label, boston_feature], axis=1) return bostondata = data_prepare()print(data.shape)(506, 14)# 查看現有特徵分布plot_dist(data, data.columns)

# 使用log1p變換將特徵基本拉到一個尺度進行建模for col in data.columns.drop(TARGET): data[col] = np.log1p(data[col])# 觀察新的特徵的分布plot_dist(data, data.columns)

corr = data.corr(spearman)corrsns.heatmap(corr)

threshold = 0.85correlated_pairs= {}for col in corr: # Find correlations above the threshold above_threshold_vars = [x for x in list(corr.index[abs(corr[col]) > threshold]) if x != col] correlated_pairs[col] = above_threshold_varscorrelated_pairsplt.scatter(data[DIS], data[NOX])

corr[TARGET].sort_values().plot.barh()

特徵LSTAT 和 target 相關性較高， DIS 和 NOX 兩個特徵相關性較高,可視化展示一下

原始特徵的建模

features = data.columns.drop(TARGET)x_train, x_test, y_train, y_test = train_test_split(data[features], data[TARGET].to_frame(), test_size=0.3, shuffle=True)alphas = [0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1]for alpha in alphas: metrics, fi = Preds(x_train, y_train, x_test, y_test, alpha) print( {}.format(alpha)) print(metrics)# 驗證集上 0.01最優metrics, fi = Preds(x_train, y_train, x_test, y_test, 0.01)0.01 fold train valid test0 0 18.865634 15.696148 20.5579631 1 15.828226 24.284637 19.9558492 2 17.190880 20.304150 18.6935563 3 17.687572 18.251201 18.2418854 overall 17.393078 19.636047 19.362313fifi.plot.barh()

原特徵的模型貢獻

2 實驗一

構建新特徵

function_set = [add, sub, mul, div, log, sqrt, abs, neg, max, min]gp1 = SymbolicTransformer(generations=10, population_size=1000, hall_of_fame=100, n_components=10, function_set=function_set, parsimony_coefficient=0.0005, max_samples=0.9, verbose=1, random_state=0, n_jobs=3)print(gp1)[div(min(X5, X5), add(X10, X12)), sub(-0.988, div(min(X5, X5), log(div(X5, add(X10, X12))))), div(X5, add(X10, add(X10, X12))), div(div(X5, add(X10, X12)), add(X10, X12)), div(log(log(div(X5, add(X10, X12)))), add(X10, X12)), div(X5, log(sqrt(add(X10, X12)))), div(sqrt(div(X5, add(X10, X12))), add(X10, X12)), div(X5, add(add(X10, X12), X12)), log(div(X5, add(X10, X12))), div(div(X5, add(X10, X12)), log(div(sub(div(X5, add(X10, X12)), min(X5, X5)), add(X10, X12))))]from IPython.display import Imageimport pydotplusgraph = gp1._best_programs[0].export_graphviz()graph = pydotplus.graphviz.graph_from_dot_data(graph)Image(graph.create_png())

可視化特徵組合的方式，可以看到隨機性是很強的，這個圖示左子樹其實沒什麼用。。

gp_train_feature = gp1.transform(x_train)gp_test_feature = gp1.transform(x_test)new_feature_name = [str(i)+V for i in range(1, 11)]train_new_feature = pd.DataFrame(gp_train_feature, columns=new_feature_name, index=train_idx)test_new_feature = pd.DataFrame(gp_test_feature, columns=new_feature_name, index=test_idx)x_train_0 = pd.concat([x_train, train_new_feature], axis=1)x_test_0 = pd.concat([x_test, test_new_feature], axis=1)new_x_data = pd.concat([x_train_0, x_test_0], axis=0)new_data = pd.concat([data[TARGET], new_x_data], axis=1)new_data.columnsIndex([TARGET, CRIM, ZN, INDUS, CHAS, NOX, RM, AGE, DIS, RAD, TAX, PTRATIO, B, LSTAT, 1V, 2V, 3V, 4V, 5V, 6V, 7V, 8V, 9V, 10V], dtype=object)

新特徵的分布

特徵之間的相關性

新特徵與label的相關性都很高

特徵與標籤之間的相關性

新特徵之間的相關性很高, 都在0.95以上

alphas = [1e-4, 5e-4, 0.001, 0.005, 0.01, 0.05, 0.1, 0.3, 0.5, 0.7, 1]for alpha in alphas: metrics, fi = Preds(x_train_0, y_train, x_test_0, y_test, alpha) print( {}.format(alpha)) print(metrics)# 驗證集上 0.001最優metrics, fi = Preds(x_train_0, y_train, x_test_0, y_test, 0.001)0.001 fold train valid test0 0 13.716057 10.485317 19.6846151 1 10.505336 20.059070 19.6797002 2 11.941672 15.593777 18.4747923 3 13.531627 10.720574 18.7300454 overall 12.423673 14.220659 19.142288

特徵貢獻

對比元特徵訓練的 mse分數，val 19.636047 test 19.362313 性能提升有限

結論

在這個例子中，實際上模型效果沒有提升，主要原因感覺是生成的特徵屬性較為單一，遺傳演算法最終的目標都是擬合label，導致特徵之間相關性太高。

感覺可以適當降低遺傳演算法的輪數，保存特徵的多樣性，降低特徵之間的相關性。

3 實驗二

將遺傳代數減低至一代進行模型構建

function_set = [add, sub, mul, div, log, sqrt, abs, neg, max, min] #gp2 = SymbolicTransformer(generations=1, population_size=1000, hall_of_fame=100, n_components=10, function_set=function_set, parsimony_coefficient=0.0005, max_samples=0.9, verbose=1, random_state=0, n_jobs=3)

特徵之間的相關性

空值遺傳1代(那不就是沒有遺傳) 這樣的相關性才變低一點

特徵與標籤之間的的相關性

相比實驗一，與標籤之間的相關性降低了。

alphas = [0.001, 0.005, 0.01, 0.05, 0.1, 0.3, 0.5, 0.7, 1]for alpha in alphas: metrics, fi = Preds(train_new_feature_1, y_train, test_new_feature_1, y_test, alpha) print( {}.format(alpha)) print(metrics)# 驗證集上 0.01最優metrics, fi = Preds(train_new_feature_1, y_train, test_new_feature_1, y_test, 0.01)0.01 fold train valid test0 0 13.407300 12.753712 17.1503281 1 10.368547 25.864163 17.4952242 2 12.626428 14.408992 16.0773903 3 12.959124 13.426672 16.6949904 overall 12.340350 16.628614 16.854483

特徵貢獻

結論

實驗一和實驗二分別只是調整了遺傳演算法的訓練代數，實驗一訓練了三代，生成的特徵與目標相關性高，但是特徵之間的相關性也很高，實驗二隻訓練了一代其實就是沒有進行遺傳迭代，雖然生成的特徵與目標相關性不強，但是生成的特徵具有更高的複雜性，反而在更好的提升了模型的性能，也有幾個特徵在模型中的貢獻很大。

所以可以兩個策略，一種是訓練多代，輸出少量特徵；另一種是不進行遺傳迭代，隨機產生大量特徵，綜合兩種特徵進行建模。

所以實驗三中，綜合前兩個實驗gplearn模型生成特徵的方式，看是否有模型提升效果。

4 實驗三

gp3 = SymbolicTransformer(generations=1, population_size=1000, hall_of_fame=100, n_components=10, function_set=function_set, parsimony_coefficient=0.0005, max_samples=0.9, verbose=1, random_state=0, n_jobs=3)gp4 = SymbolicTransformer(generations=3, population_size=1000, hall_of_fame=100, n_components=1, # 生成特徵生成減少到2 function_set=function_set, parsimony_coefficient=0.0005, max_samples=0.9, verbose=1, random_state=0, n_jobs=3)

特徵之間的相關性

特徵與標籤之間的相關性

alphas = [0.001, 0.005, 0.01, 0.05, 0.1, 0.3, 0.5, 0.7, 1]for alpha in alphas: metrics, fi = Preds(train_new_feature_34, y_train, test_new_feature_34, y_test, alpha) print( {}.format(alpha)) print(metrics)# 驗證集上 0.01最優metrics, fi = Preds(train_new_feature_34, y_train, test_new_feature_34, y_test, 0.01)0.01 fold train valid test0 0 13.401289 12.729824 17.1477191 1 10.356919 25.810338 17.4734782 2 12.604374 14.398953 16.0520303 3 12.951852 13.393118 16.6772304 overall 12.328608 16.598239 16.837614

特徵貢獻

雖然M類型的特徵與標籤的相關性沒有N類特徵高，但是在模型中有可能貢獻較強。

結論

性能稍微提升了一點，當然還可以繼續增加生成的特徵數量看看能不能繼續提升，這裡就不再繼續了。

當然後續還有空值gplearn模型複雜的的操作方式，這裡就沒有繼續研究了，設想一下，如果控制了生成公式的複雜程度，是可以適當提升一下遺傳迭代的層數來生成一些保存了多樣性的特徵的。

最後再和featuetools對比一下：

ft 的強項還是在於鏈接多個表的時候生成大量特徵組合，但是在一個表的層級中，不能表現出transform 和特徵組合的多樣性。

gplearn 的強項在於在單個表的層級之中形成大量隨機的特徵組合，並且能夠有一定的篩選機制。

具體代碼可以到這裡看看 notebook