我想在我的即将推出的模型中使用python中的xgboost。但是,由于我们的生产系统是在SAS中,我试图从xgboost中提取决策规则,然后编写SAS评分代码以在SAS环境中实现此模型。
我已经通过了多个链接。以下是其中一些:
How to extract decision rules (features splits) from xgboost model in python3?
以上两个链接特别有助于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)我发现了与伽玛(正则化术语)相关的东西。它是否会影响xgboost使用叶值计算分数的方式。
3)转储文件给出的叶值是否会有任何四舍五入,从而产生此问题
此外,除了解析转储文件之外,我还要感谢执行此任务的任何其他方法。
P.S。:我只有SAS EG,无法访问SAS EM或SAS IML。
我在获得匹配分数方面有类似的经验。
我的理解是,除非你修复ntree_limit
选项以匹配你在模型拟合期间使用的n_estimators
,否则得分可能会提前停止。
df['score']= xgclfpkl.predict(df[xg_features], ntree_limit=500)
在我开始使用ntree_limit
之后,我开始获得匹配的分数。