从Python的XGBoost中提取决策规则

6

我想在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)dump文件给出的叶值是否会进行四舍五入,从而创建此问题
此外,我会感激除解析转储文件之外的任何其他方法来完成此任务。
附注:我只有SAS EG,没有访问SAS EM或SAS IML的权限。
5个回答

1
我在得到匹配分数方面有过类似的经历。
我的理解是,如果你没有将ntree_limit选项设置为与模型拟合期间使用的n_estimators相匹配,则评分可能会提前停止。
df['score']= xgclfpkl.predict(df[xg_features], ntree_limit=500)

在我开始使用ntree_limit之后,我开始得到匹配分数。


嗨KKane,非常感谢您的评论,因为我现在陷入了困境。然而,我没有理解您所说的话。您是指XGboost会自动停止树的数量吗?您能否帮助我理解这一点。另外,看起来您已经解决了这个问题。您能否通过在此处发布代码来帮助我。非常感谢您的帮助。 - Ved
嗨KKane,你能否回复我上面的问题? - Ved

0
以下是代码片段,用于打印从xgboost模型的booster树中提取的所有规则。下面的代码假定您已经将一个模型打包成了一个pickle文件。
import pandas as pd
import numpy as np
import pickle
import networkx as nx

_model = pickle.load(open(MODEL_FILE, "rb"))

df = _model._Booster.trees_to_dataframe()
df['_missing'] = df.apply(
    lambda x: 'Yes' if pd.notnull(x['Missing']) and pd.notnull(x['Yes']) and pd.notnull(x['No']) and x['Missing'] == x[
        'Yes'] else 'No', axis = 1)

G = nx.DiGraph()
G.add_nodes_from(df.ID.tolist())

yes_edges = df[['ID', 'Yes', 'Feature', 'Split', '_missing']].dropna()
yes_edges['label'] = yes_edges.apply(
    lambda x: "({feature} < {value:.4f} or {feature} is null)".format(feature = x['Feature'], value = x['Split']) if x['_missing'] == 'Yes'
    else "({feature} < {value:.4f})".format(feature = x['Feature'], value = x['Split']),
    axis = 1
)

no_edges = df[['ID', 'No', 'Feature', 'Split', '_missing']].dropna()
no_edges['label'] = no_edges.apply(
    lambda x: "({feature} >= {value:.4f} or {feature} is null)".format(feature = x['Feature'], value = x['Split']) if x['_missing'] == 'No'
    else "({feature} >= {value:.4f})".format(feature = x['Feature'], value = x['Split']),
    axis = 1
)

for v in yes_edges.values:
    G.add_edge(v[0], v[1], feature = v[2], expr = v[5])

for v in no_edges.values:
    G.add_edge(v[0], v[1], feature = v[2], expr = v[5])

leaf_node_score_values = {i[0]: i[1] for i in df[df.Feature == 'Leaf'][['ID', 'Gain']].values}
nodeID_to_tree_map = {i[1]: i[0] for i in df[['Tree', 'ID']].values}

roots = []
leaves = []
for node in G.nodes:
    if G.in_degree(node) == 0:  # it's a root
        roots.append(node)
    elif G.out_degree(node) == 0:  # it's a leaf
        leaves.append(node)

paths = []
for root in roots:
    for leaf in leaves:
        for path in nx.all_simple_paths(G, root, leaf):
            paths.append(path)

rules = []
temp = []
for path in paths:
    parts = []
    for i in range(len(path) - 1):
        parts.append(G[path[i]][path[i + 1]]['expr'])
    rules.append(" and ".join(parts))
    temp.append((
        path[0],
        nodeID_to_tree_map.get(path[0]),
        " and ".join(parts),
        leaf_node_score_values.get(path[-1])
    ))

rules_df = pd.DataFrame.from_records(temp, columns = ['node', 'tree', 'rule', 'score'])
rules_df['prob'] = rules_df.apply(lambda x: 1 / (1 + np.exp(-1 * x['score'])), axis = 1)
rules_df['rule_idx'] = rules_df.index
rules_df = rules_df.drop(['node'], axis = 1)

print("n_rules -> {}".format(len(rules_df)))

del G, df, roots, leaves, yes_edges, no_edges, temp, rules

上述代码按以下格式打印每个规则:
if x>y and a>b and c<d then e

0

我有一个类似的经验,需要从R转换xgboost评分代码到SAS。

最初,我遇到了你在这里遇到的相同问题,也就是说,在较小的树中,在R和SAS中的得分之间没有太大的差异,但一旦树的数量增加到100或更多,我开始观察到不一致性。

我做了三件事来缩小这些差异:

  1. 确保缺失组朝正确方向前进,您需要明确说明。否则,SAS会将缺失值视为所有数字中最小的值。规则应该像下面这样在SAS中:

if sepal_width > 2.95000005 or missing(sepal_width) then state = 1;else state = 2;
或者
if sepal_width <= 2.95000005 and ~missing(sepal_width) then state = 1;else state = 2;

  1. 我使用了一个名为float的R包,使得分数有更多的小数位。 as.numeric(float::fl(Quality))

  2. 确保SAS数据与您在Python中训练的数据具有相同的形状。

希望以上内容能够帮到您。


嗨DDZR,谢谢你的回复。我已经处理了1和3点。但是我不理解第2点。我正在从转储文件中读取浮点值,如何在此处增加小数点? - Ved

0
几点需要注意:
首先,用于匹配叶子返回值的正则表达式无法捕获转储中的“e-小数”科学计数法(默认)。以下是一个明确的示例(第二个是正确的修改!)-
s = '3:leaf=9.95066429e-09'
out = re.findall(r"[\d.-]+", s)
out2 = re.findall(r"-?[\d.]+(?:e-?\d+)?", s)
out2,out

(容易修复,但很难发现,因为我的模型只受到了一个叶子的影响!)

其次,这个问题是关于二进制的,但在多类目标中,转储文件中有每个类别的单独树,所以总共有T*C棵树,其中T是增强轮数,C是类别数。对于类别c(在{0,1,...,C-1}中),您需要评估(并求和终端叶子)树i*C +c,其中i = 0,...,T-1。然后进行softmax操作以匹配xgb的预测结果。


非常感谢 P.Windridge 对于从转储文件中读取科学记数法所做的更正。我已经进行了必要的更改。然而,在上述示例中并没有与科学记数法相关的错误,问题仍然存在。请问您能告诉我该如何解决吗?附注:我的数据中没有任何缺失值。 - Ved

0

我稍微研究了一下如何将这个功能加入到我的代码中。

我发现有一个小问题,就是缺少处理。

如果你的逻辑类似于这样,它似乎可以正常工作:

if petal_length < 3 or missing(petal_length) then state = 5;
        else state = 6;

但是说缺失的组应该去第6个状态而不是第5个状态。那么你会得到这样的代码:

if petal_length < 3 then state = 5;
        else state = 6;

在这种情况下,petal_length = missing (.)会处于什么状态? 嗯,在SAS中,missing被归类为小于任何数字,因此它仍然进入状态5(而不是预期的状态6)。
要解决这个问题,您可以将所有缺失值分配给999999999999999(选择一个高数字,因为XGBoost格式始终使用小于符号(<)),然后替换。
missing_value_handling = (" or missing(" + out[1] + ")")

使用

missing_value_handling = (" or " + out[1] + "=999999999999999 ")

在您的 string_parser 中。

非常感谢您的建议,David。我已经进行了这些更改。然而,在我的数据中没有缺失的观测值,问题仍然存在,即得分在第四位小数点后仍然不匹配。代码正确读取转储文件,因此只有在转储文件本身存在小数点错误的情况下才会发生这种情况。然而,我肯定不相信这可能是情况,因此不确定该怎么办。您能否请提供一些建议,可能导致这种情况。非常感谢您的帮助。 - Ved

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接