从python中的xgboost中提取决策规则

Question

我想在 python 中使用 xgboost 来构建我即将推出的模型。然而，由于我们的生产系统是在SAS中，我试图从xgboost中提取决策规则，然后编写SAS评分代码以在SAS环境中实现该模型。

我已经浏览了多个链接。以下是其中一些：

如何从python3中的xgboost模型中提取决策规则（特征分割）？

xgboost部署

上面的两个链接非常有帮助，特别是 Shiutang-Li 给出的用于 xgboost 部署的代码。然而，我的预测分数并不完全匹配。

下面是我迄今为止尝试过的代码：

import numpy as np
import pandas as pd
import xgboost as xgb
from sklearn.grid_search import GridSearchCV
%matplotlib inline
import graphviz
from graphviz import Digraph

#Read the sample iris data:
iris =pd.read_csv("C:\\Users\\XXXX\\Downloads\\Iris.csv")
#Create dependent variable:
iris.loc[iris["class"] != 2,"class"] = 0
iris.loc[iris["class"] == 2,"class"] = 1

#Select independent and dependent variable:
X = iris[["sepal_length","sepal_width","petal_length","petal_width"]]
Y = iris["class"]

xgdmat = xgb.DMatrix(X, Y) # Create our DMatrix to make XGBoost more efficient

#Build the sample xgboost Model:

our_params = {'eta': 0.1, 'seed':0, 'subsample': 0.8, 'colsample_bytree': 0.8, 
             'objective': 'binary:logistic', 'max_depth':3, 'min_child_weight':1} 
Base_Model = xgb.train(our_params, xgdmat, num_boost_round = 10)

#Below code reads the dump file created by xgboost and writes a scoring code in SAS:

import re
def string_parser(s):
    if len(re.findall(r":leaf=", s)) == 0:
        out  = re.findall(r"[\w.-]+", s)
        tabs = re.findall(r"[\t]+", s)
        if (out[4] == out[8]):
            missing_value_handling = (" or missing(" + out[1] + ")")
        else:
            missing_value_handling = ""

        if len(tabs) > 0:
            return (re.findall(r"[\t]+", s)[0].replace('\t', '    ') + 
                    '        if state = ' + out[0] + ' then do;\n' +
                    re.findall(r"[\t]+", s)[0].replace('\t', '    ') +
                    '            if ' + out[1] + ' < ' + out[2] + missing_value_handling +
                    ' then state = ' + out[4] + ';' +  ' else state = ' + out[6] + ';\nend;' ) 
        else:
            return ('        if state = ' + out[0] + ' then do;\n' +
                    '            if ' + out[1] + ' < ' + out[2] + missing_value_handling +
                    ' then state = ' + out[4] + ';' +  ' else state = ' + out[6] + ';\nend;' )
    else:
        out = re.findall(r"[\w.-]+", s)
        return (re.findall(r"[\t]+", s)[0].replace('\t', '    ') + 
                '        if state = ' + out[0] + ' then\n    ' +
                re.findall(r"[\t]+", s)[0].replace('\t', '    ') + 
                '        value = value + (' + out[2] + ') ;\n')

def tree_parser(tree, i):
    return ('state = 0;\n'
             + "".join([string_parser(tree.split('\n')[i]) for i in range(len(tree.split('\n'))-1)]))

def model_to_sas(model, out_file):
    trees = model.get_dump()
    result = ["value = 0;\n"]
    with open(out_file, 'w') as the_file:
        for i in range(len(trees)):
            result.append(tree_parser(trees[i], i))
        the_file.write("".join(result))
        the_file.write("\nY_Pred1 = 1/(1+exp(-value));\n")
        the_file.write("Y_Pred0 = 1 - Y_pred1;") 

调用上述模块创建SAS评分代码：

model_to_sas(Base_Model, 'xgb_scr_code.sas')

不幸的是，我无法提供上述模块生成的完整 SAS 代码。但是，如果我们仅使用一种树代码构建模型，请在下面找到 SAS 代码：

value = 0;
state = 0;
if state = 0 then
    do;
        if sepal_width < 2.95000005 or missing(sepal_width) then state = 1;
        else state = 2;
    end;
if state = 1 then
    do;
        if petal_length < 4.75 or missing(petal_length) then state = 3;
        else state = 4;
    end;

if state = 3 then   value = value + (0.1586207);
if state = 4 then   value = value + (-0.127272725);
if state = 2 then
    do;
        if petal_length < 3 or missing(petal_length) then state = 5;
        else state = 6;
    end;
if state = 5 then   value = value + (-0.180952385);
if state = 6 then
    do;
        if petal_length < 4.75 or missing(petal_length) then state = 7;
        else state = 8;
    end;
if state = 7 then   value = value + (0.142857149);
if state = 8 then   value = value + (-0.161290333);

Y_Pred1 = 1/(1+exp(-value));
Y_Pred0 = 1 - Y_pred1;

以下是第一棵树的转储文件输出：

booster[0]:
    0:[sepal_width<2.95000005] yes=1,no=2,missing=1
        1:[petal_length<4.75] yes=3,no=4,missing=3
            3:leaf=0.1586207
            4:leaf=-0.127272725
        2:[petal_length<3] yes=5,no=6,missing=5
            5:leaf=-0.180952385
            6:[petal_length<4.75] yes=7,no=8,missing=7
                7:leaf=0.142857149
                8:leaf=-0.161290333

所以基本上，我想做的是将节点号保存在变量“state”中，并相应地访问叶节点（这是我从上面链接中提到的 Shiutang-Li 的文章中学到的）。

这是我面临的问题：

对于最多大约 40 棵树，预测分数完全匹配。例如请看下面：

情况1：

使用 python 对 10 棵树进行预测的值：

Y_pred1 = Base_Model.predict(xgdmat)

print("Development- Y_Actual: ",np.mean(Y)," Y predicted: ",np.mean(Y_pred1))

输出：

Average- Y_Actual:  0.3333333333333333  Average Y predicted:  0.4021197

使用 SAS 对 10 棵树的预测值：

Average Y predicted:  0.4021197

案例2：

使用 python 对 100 棵树的预测值：

Y_pred1 = Base_Model.predict(xgdmat)

print("Development- Y_Actual: ",np.mean(Y)," Y predicted: ",np.mean(Y_pred1))

输出：

Average- Y_Actual:  0.3333333333333333  Average Y predicted:  0.33232176

使用 SAS 对 100 棵树的预测值：

Average Y predicted:  0.3323159

正如您所看到的，100 棵树的分数并不完全匹配（最多匹配小数点后 4 位）。另外，我在大文件上尝试过此方法，其中分数差异相当高，即分数偏差超过 10%。

谁能让我指出代码中的任何错误，以便分数可以完全匹配。以下是我的一些疑问：

1）我的分数计算正确吗？

2）我发现了一些与gamma（正则化项）相关的东西。它是否影响 xgboost 使用叶值计算分数的方式。

3）转储文件给出的叶子值是否会进行四舍五入，从而产生此问题

另外，除了解析转储文件之外，我还希望有任何其他方法来完成此任务。

PS：我只有 SAS EG，无权访问 SAS EM 或 SAS IML。

Answer 1

我在获得匹配分数方面也有类似的经历。我的理解是，除非您修复选项以匹配您在模型拟合期间使用的选项，
否则评分可能会提前停止。ntree_limitn_estimators

df['score']= xgclfpkl.predict(df[xg_features], ntree_limit=500)

在我开始使用后ntree_limit，我开始获得匹配的分数。