我想在python中使用xgboost来构建模型。然而,由于我们的生产系统是基于SAS的,因此我正在尝试从xgboost中提取决策规则,然后编写SAS评分代码以在SAS环境中实现该模型。
我已经查阅了多个链接,以下是其中一些:
如何从Python3中的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的权限。