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

我想在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] + ' thenn    ' +
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)我发现了与伽马(正则化项)有关的东西。它是否会影响 xgboost 使用叶值计算分数的方式。

3)转储文件给出的叶值是否会有任何舍入，从而产生此问题

另外，除了解析转储文件之外，我将不胜感激任何其他方法来执行此任务。

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

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

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

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

我有类似的经历，需要将 xgboost 评分代码从 R 提取到 SAS。

最初，我遇到了和你在这里遇到的同样的问题，即在较小的树中，R 和 SAS 中的分数之间没有太大差异，一旦树的数量上升到 100 或更多，我就开始观察差异。

我做了 3 件事来缩小差异：

确保失踪的小组朝着正确的方向前进，您需要明确。否则，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;

我使用了一个名为float的 R 包来使分数具有更多的小数位。as.numeric(float::fl(Quality))
确保 SAS 数据的形状与在 Python 中训练的数据相同。

希望以上内容对您有所帮助。

我稍微考虑过将其合并到我自己的代码中。

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

它似乎在你有这样的逻辑的地方工作正常

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 (.)会得到什么状态？好吧，在这里它仍然进入状态 5(而不是预期的状态 6)，因为在 SAS 中，缺失被归类为小于任何数字。

要解决此问题，您可以将所有缺失值分配给999999999999999(选择一个较大的数字，因为 XGBoost 格式始终使用小于 (<))，然后替换

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

跟

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

在你的string_parser.

几点-

首先，匹配叶返回值的正则表达式不会捕获转储中的"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。然后软最大化它以匹配来自 xgb 的预测。

下面是代码片段，它打印了从 xgboost 模型中的助推器树中提取的所有规则。下面的代码假定你已经将模型打包到一个 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

相关内容

最新更新

热门标签：