第 8 章：整合開源

此章節主要介紹 SAS 平台如何整合開源程式碼 R 和 Python。

基本介紹

首先 SAS Cloud Analytic Service (CAS) 雲端分析服務為在 SAS Viya 分析平台的高效能處理資料管理和分析的軟體和硬體，簡稱為 CAS 伺服器，並且我們能夠透過 SAS Viya 分析平台進行資料管理，CAS 伺服器能夠多種資料來源載入資料至記憶體中，而當刪除記憶體中的資料表時，其不會影響資料來源中的資料檔案，所以 CAS 伺服器僅會將需要使用的資料自動載入記憶體中，以利後續進行高效能處理，此外 CAS 伺服器具有以下特點，分別為：

支援開放源始碼 Python 和 R 授權存取資料。
針對資料統一 UTF-8 編碼格式。
針對資料身份驗證和存取控管。
針對資料進行安全傳輸和加密儲存。
在多台伺服器記憶體中進行高效能處理資料。

接著 CAS 伺服器主要支援 Python 支援 3.4.0 以上版本和 R 支援 3.1.0 以上版本，並且我們主要皆能夠使用 SAS Scripting Wrapper for Analytics Transfer (SWAT) 介面連接至 CAS 伺服器，以及對於熟悉開放源始碼 Python 和 R 的資料科學家將能夠透過 Python 和 R 所對應的 SWAT 套件函式庫授權存取 SAS Viya 分析平台中已被管理的資料。此外 CAS 伺服器僅支援 UTF-8 統一編碼格式，所以我們必須將讀取至 CAS 伺服器中資料轉碼為 UTF-8 編碼格式，預設 CAS Session 設定為 UTF-8 能夠避免資料轉碼問題。至於所有在 CAS 伺服器中的資料皆能夠透過 CAS 資料館進行所有操作，同時 CAS 資料館針對多種資料來源，像是資料庫、檔案系統目錄和記憶體中檔案的存取，以及針對 CAS 資料館進行關聯以管理資料的身份驗證、存取控管、安全傳輸和加密儲存，當然 CAS 伺服器中的所有資料皆會被載入至記憶體中進行高效能的處理。

再來在 SAS Viya 分析平台中的資料，主要是儲存至 CAS 伺服器記憶體中的資料表，若載入至 Python 或 R 程式語言中使用則為 CASTable 物件，我們除了能夠透過 SWAT 套件函式庫授權針對 CASTable 物件進行資料操作和資料探索，更能夠將 CASTable 物件轉換為 Python 和 R 的 DataFrame 物件，以利 Python 和 R 的資料科學家能夠將 DataFrame 物件繼續用於訓練 Python 和 R 所提供最新的機器學習模型，至於要如何讓 Python 和 R 的資料科學家能夠透過 SWAT 套件函式庫授權存取 SAS Viya 分析平台中已被管理的資料，請參考下指令對照表。

最後我們更能夠透過 SAS Viya 分析平台為企業組織落實資料治理，而 SAS 資料治理框架主要提供了全面的功能，其將能夠符合各種規模企業組織結構，並且建立和維持有效的資料治理計劃，以利能夠為企業組織中的所有使用者提供可被信任，即時性和高品質的資料，帶來更高的商業價值，至於更多有關 SAS 資料治理框架的詳細資訊，請參考官方白皮書。

總結我們將能夠以 SAS Viya 分析平台將資料載入至 CAS 伺服器中進行資料管理，並且透過 Python 和 R 授權存取 CAS 伺服器中的資料，以利用於訓練 Python 和 R 所提供最新的機器學習模型，以及透過 SAS Viya 分析平台協助企業組織達到資料治理的願景。

整合應用

首先透過 SAS Viya 平台中的應用程式介面 (API) 將能夠幫助企業客戶在核心業務應用程式中使用人工智慧的強大分析功能，而為什麼 API 很重要呢？主要有三個重點，分別為：

結果很重要：核心應用程式使用人工智慧模型需要有效地解決問題。
生產很重要：核心應用程式自動化取得人工智慧模型分析結果。
開放很重要：核心應用程式從開放的分析平台取得分析結果。

接著我們主要有三種方式使用 SAS Viya 平台的 API，分別為：

SAS Job Execution：針對應用程式快速建立基於 SAS 程式的網路服務。
Full stack REST API apps：與應用程式無逢整合 SAS Viya 分析平台的功能。
Flask App with SWAT package：應用程式整合基於網站伺服器介面與 CAS 伺服器互動。

其中 SAS Job Execution 主要是將 SAS 程式封裝在網路服務中，並且能夠在 CAS 或SPRE 中執行。Full stack REST API apps 主要是採用標準 REST API 透過 javascript 或 python 進行整合應用。Flask App with SWAT package 主要是網站服務器和 API 腳本皆使用 Python 通用語言，一旦開發工作完成，就能夠在不同的路徑下分配所需的腳本和使用 SWAT 套件與 CAS 伺服器進行互動。

再來企業應用 API 的案例主要有三種類型，分別為：

產生基於規則或評分的商業決策。
建立推薦或警示通知系統。
整合文字分析實務應用。

其中基於規則或評分的商業決策所面臨的挑戰主要有直接來自客戶的輸入、需要快速決定以及可能涉及多層商業邏輯，此時我們就能夠透過 API 整合 SAS Viya 平台中 Micro Analytics Services 的功能解決所面臨的問題。建立推薦或警示通知系統所面臨的挑戰主要有直接來自交易作為輸入、許多不同的模式以及快速持續變化，此時我們就能夠透過 API 整合 SAS Viya 平台中 ASTORE 模型的功能解決所面臨的問題。文字分析整合應用所面臨的挑戰主要有提取和分類文字資料和許多不同文字的複雜解析，此時我們就能夠透過 API 整合 SAS Viya 平台中 Natural Language Processing模型的功能解決所面臨的問題。

最後當我們選擇以 API 的方式整合 SAS Viya 平台之後，更需針對授權、合作和管理進行控管，其中授權是以 OAuth 2.0 的方式為主，合作則建議參考 SAS Software 的 Github 開放源始碼，以及管理除了透過 SAS Viya 平台管理資料和內容之外，更進一步建議透過企業專門的 API 管理工具針對 SAS Viya 平台的 API 進行控管。

開始使用

整合 SAS Visual Data Mining and Machine Learning 中的管線

開源程式碼 R

登入至 SAS Viya 分析平台。

在畫面左上方點選【顯示應用程式功能表】鈕。

選取【建置模型】。

點選【新增專案】鈕。

輸入【名稱】為【HMEQ_R】。

選取【類型】為【資料採礦和機器學習】。

按下【瀏覽】鈕，選擇【資料】。

點選【HMEQ】資料，按下【確定】鈕。

按下【進階】鈕。

選取【分割資料】。

將訓練、驗證和測試，分別修改為【70】、【20】和【10】，按下【儲存】。

按下【儲存】。

自動識別【BAD】變數名稱為【目標】角色，其它變數為【輸入】角色。

點選【管線】頁籤。

點選【節點】。

展開【資料採礦前置處理】節點。

在【補值】節點上按住左鍵，拖曳至【資料】節點上，請注意要呈現【+ 補值】。

在【補值】節點上按下右鍵，點選【增加子系節點】，點選【其他】，點選【開源程式碼】。

將【開源程式碼】節點中的【語言】改為【R】。

按下【開啟程式碼編輯器】鈕。

貼上以下【R 程式碼】。

library(randomForest)

# RandomForest
dm_model <- randomForest(dm_model_formula, ntree=100, mtry=5, data=dm_traindf, importance=TRUE)

# Score
pred <- predict(dm_model, dm_inputdf, type="prob")
dm_scoreddf <- data.frame(pred)
colnames(dm_scoreddf) <- c("P_BAD0", "P_BAD1")

# Print/plot model output
png("rpt_forestMsePlot.png")
plot(dm_model, main='randomForest MSE Plot')
dev.off()

# Print Variable Importance
write.csv(importance(dm_model), file="rpt_forestIMP.csv", row.names=TRUE)

R 程式碼，請參考官方網站。

按下【儲存】鈕，按下【關閉】鈕。

按下【執行管線】鈕。

已成功完成管線執行。

在【開源程式碼】節點上按右鍵點選【結果】。

查看【開源程式碼】節點執行的結果【摘要】。

在【R 程式碼】上按下【展開】鈕，查看自動產生的【R 程式碼】，按下【關閉】鈕。

#-------------------------------------------------------------------------------
# 語言: R
#-------------------------------------------------------------------------------
# 將 R 工作設為節點工作目錄
dm_nodedir <- '/opt/sas/viya/config/var/tmp/compsrv/default/f51e1c4c-5ce7-44f8-88b1-cac06c0a6014/SAS_workB9A100001364_sasserver.demo.sas.com/81aed726-83c9-4bf0-95fc-417b6121a8c7'
setwd(dm_nodedir)

# 變數宣告
dm_dec_target <- 'BAD'
dm_partitionvar <- '_PartInd_'
dm_partition_train_val <- 1

dm_class_input <- c("IMP_DELINQ","IMP_DEROG","IMP_JOB","IMP_NINQ","IMP_REASON")
dm_interval_input <- c("IMP_CLAGE","IMP_CLNO","IMP_DEBTINC","IMP_MORTDUE","IMP_VALUE","IMP_YOJ"    ,"LOAN")
dm_input <- c(dm_class_input, dm_interval_input)
dm_model_formula <- as.formula(paste(dm_dec_target, '~', paste(dm_input, collapse='+')))
dm_predictionvar <- c('P_BAD0', 'P_BAD1')
dm_classtarget_intovar <- 'I_BAD'
dm_classtarget_level <- c('0', '1')

#-------------------------------------------------------------------------------
# 產生資料框架: Y
#-------------------------------------------------------------------------------
dm_inputdf <- read.csv(file='node_data.csv', stringsAsFactors=TRUE, check.names=FALSE)

# Change class variable type to R factor
dm_inputdf$BAD <- factor(dm_inputdf$BAD, ordered=FALSE)
dm_inputdf$IMP_DELINQ <- factor(dm_inputdf$IMP_DELINQ, ordered=FALSE)
dm_inputdf$IMP_DEROG <- factor(dm_inputdf$IMP_DEROG, ordered=FALSE)
dm_inputdf$IMP_NINQ <- factor(dm_inputdf$IMP_NINQ, ordered=FALSE)

dm_traindf = subset(dm_inputdf, get(dm_partitionvar) == dm_partition_train_val)

#-------------------------------------------------------------------------------
# 使用者程式碼
#-------------------------------------------------------------------------------
library(randomForest)

# RandomForest
dm_model <- randomForest(dm_model_formula, ntree=100, mtry=5, data=dm_traindf, importance=TRUE)

# Score
pred <- predict(dm_model, dm_inputdf, type="prob")
dm_scoreddf <- data.frame(pred)
colnames(dm_scoreddf) <- c("P_BAD0", "P_BAD1")

# Print/plot model output
png("rpt_forestMsePlot.png")
plot(dm_model, main='randomForest MSE Plot')
dev.off()

# Print Variable Importance
write.csv(importance(dm_model), file="rpt_forestIMP.csv", row.names=TRUE)

選取【輸出資料】頁籤，按下【輸出資料】頁籤，按下【檢視輸出資料】鈕。

按下【檢視輸出資料】鈕。

點選【變數表格】鈕，查看補值變數【IMP_XXX】、分割指標變數【_PartInd_】和索引值【_dmIndex_】，按下【關閉】鈕。

在【資料】節點上按下右鍵點選【增加子系節點】，點選【監督式學習】，點選【梯度提升】。

在【開源程式碼】節點上按下右鍵點選【移動】，點選【監督式學習】。

按下【執行管線】鈕。

按下【模型比較】節點上按下右鍵點選【結果】。

查看【模型比較】結果，其中【最佳模型】為【開源程式碼】節點的演算法。

點選【評估】頁籤，查看【提升報表】、【ROC 報表】和【配適統計】的評估資訊，按下【關閉】。

點選【管線比較】頁籤。

點選【觀點】頁籤。

開放源始碼 Python

登入至 SAS Viya 分析平台。

在畫面左上方點選【顯示應用程式功能表】鈕。

選取【建置模型】。

點選【新增專案】鈕。

輸入【名稱】為【HMEQ_Python】。

選取【類型】為【資料採礦和機器學習】。

按下【瀏覽】鈕，選擇【資料】。

點選【HMEQ】資料，按下【確定】鈕。

按下【進階】鈕。

選取【分割資料】。

將訓練、驗證和測試，分別修改為【70】、【20】和【10】，按下【儲存】。

按下【儲存】。

自動識別【BAD】變數名稱為【目標】角色，其它變數為【輸入】角色。

點選【管線】頁籤。

點選【節點】。

展開【資料採礦前置處理】節點。

在【補值】節點上按住左鍵，拖曳至【資料】節點上，請注意要呈現【+ 補值】。

在【補值】節點上按下右鍵，點選【增加子系節點】，點選【其他】，點選【開源程式碼】。

將【開源程式碼】節點中的【語言】改為【Python】。

按下【開啟程式碼編輯器】鈕。

貼上以下【Python 程式碼】。

from sklearn import ensemble

# Get full data with inputs + partition indicator
dm_input.insert(0, dm_partitionvar)
fullX = dm_inputdf.loc[:, dm_input]

# Dummy encode class variables
fullX_enc = pd.get_dummies(fullX, columns=dm_class_input, drop_first=True)

# Create X (features/inputs); drop partition indicator
X_enc = fullX_enc[fullX_enc[dm_partitionvar] == dm_partition_train_val]
X_enc = X_enc.drop(dm_partitionvar, 1)

# Create y (labels)
y = dm_traindf[dm_dec_target]

# Fit RandomForest model w/ training data
params = {'n_estimators': 100, 'max_depth': 20, 'min_samples_leaf': 5}
dm_model = ensemble.RandomForestClassifier(**params)
dm_model.fit(X_enc, y)
print(dm_model)

# Save VariableImportance to CSV
varimp = pd.DataFrame(list(zip(X_enc, dm_model.feature_importances_)),
columns=['Variable Name', 'Importance'])
varimp.to_csv(dm_nodedir + '/rpt_var_imp.csv', index=False)

# Score full data
fullX_enc = fullX_enc.drop(dm_partitionvar, 1)
dm_scoreddf = pd.DataFrame(dm_model.predict_proba(fullX_enc), columns=['P_BAD0',
'P_BAD1'])

按下【儲存】鈕，按下【關閉】鈕。

按下【執行管線】鈕。

已成功完成管線執行。

在【開源程式碼】節點上按右鍵點選【結果】。

查看【開源程式碼】節點執行的結果【摘要】。

在【Python 程式碼】上按下【展開】鈕，查看自動產生的【Python 程式碼】，按下【關閉】鈕。

#-------------------------------------------------------------------------------
# 語言: PYTHON
#-------------------------------------------------------------------------------
# 變數宣告
dm_nodedir = '/opt/sas/viya/config/var/tmp/compsrv/default/253e3264-f8ca-4208-85ab-d8348e6976ab/SAS_workFCD700005193_sasserver.demo.sas.com/0a03b320-897b-483a-b31b-668e2be5a18e'
dm_dec_target = 'BAD'
dm_partitionvar = '_PartInd_'
dm_partition_train_val = 1

dm_class_input = ["IMP_DELINQ","IMP_DEROG","IMP_JOB","IMP_NINQ","IMP_REASON"]
dm_interval_input = ["IMP_CLAGE","IMP_CLNO","IMP_DEBTINC","IMP_MORTDUE","IMP_VALUE","IMP_YOJ"    ,"LOAN"]
dm_input = dm_class_input + dm_interval_input
dm_predictionvar = ['P_BAD0', 'P_BAD1']
dm_classtarget_intovar = 'I_BAD'
dm_classtarget_level = ['0', '1']

#-------------------------------------------------------------------------------
# 產生資料框架: Y
#-------------------------------------------------------------------------------
import pandas as pd
dm_inputdf = pd.read_csv(dm_nodedir + '/node_data.csv')

dm_traindf = dm_inputdf[dm_inputdf[dm_partitionvar] == dm_partition_train_val]

#-------------------------------------------------------------------------------
# 使用者程式碼
#-------------------------------------------------------------------------------
from sklearn import ensemble 

# Get full data with inputs + partition indicator
dm_input.insert(0, dm_partitionvar)
fullX = dm_inputdf.loc[:, dm_input]

# Dummy encode class variables
fullX_enc = pd.get_dummies(fullX, columns=dm_class_input, drop_first=True)

# Create X (features/inputs); drop partition indicator
X_enc = fullX_enc[fullX_enc[dm_partitionvar] == dm_partition_train_val]
X_enc = X_enc.drop(dm_partitionvar, 1)

# Create y (labels)
y = dm_traindf[dm_dec_target]

# Fit RandomForest model w/ training data
params = {'n_estimators': 100, 'max_depth': 20, 'min_samples_leaf': 5}
dm_model = ensemble.RandomForestClassifier(**params)
dm_model.fit(X_enc, y)
print(dm_model)

# Save VariableImportance to CSV
varimp = pd.DataFrame(list(zip(X_enc, dm_model.feature_importances_)),
columns=['Variable Name', 'Importance'])
varimp.to_csv(dm_nodedir + '/rpt_var_imp.csv', index=False)

# Score full data
fullX_enc = fullX_enc.drop(dm_partitionvar, 1)
dm_scoreddf = pd.DataFrame(dm_model.predict_proba(fullX_enc), columns=['P_BAD0',
'P_BAD1'])

選取【輸出資料】頁籤，按下【輸出資料】頁籤，按下【檢視輸出資料】鈕。

按下【檢視輸出資料】鈕。

點選【變數表格】鈕，查看補值變數【IMP_XXX】、分割指標變數【_PartInd_】和索引值【_dmIndex_】，按下【關閉】鈕。

在【資料】節點上按下右鍵點選【增加子系節點】，點選【監督式學習】，點選【梯度提升】。

在【開源程式碼】節點上按下右鍵點選【移動】，點選【監督式學習】。

按下【執行管線】鈕。

按下【模型比較】節點上按下右鍵點選【結果】。

查看【模型比較】結果，其中【最佳模型】為【梯度提升】節點的演算法。

點選【評估】頁籤，查看【提升報表】、【ROC 報表】和【配適統計】的評估資訊，按下【關閉】。

點選【管線比較】頁籤。

點選【觀點】頁籤。

整合 Jupyter 筆記本中的原生開源程式碼

前往 SAS Viya Programming Trial 官方網站，按照指示進行註冊試用。

登入【Jupyter】，輸入【SAS Profile 電子郵件】和【SAS Profile 密碼】，按下【Sign In】鈕。

開始使用 Jupyter。

開源程式碼 R

這是一個簡單的端到端範例，說明如何使用SAS Viya 分析平台進行預測分析，主要有以下步驟：

載入所需的 R 套件包。
在已執行的 CAS 伺服器上啟動 CAS 工作階段。
載入所需的 CAS 操作集。
將 HMEQ 資料檔從本地文件系統上傳至 CAS 伺服器中。
探索資料。
設算缺漏值。
將資料分為訓練資料和驗證資料。
建立一個決策樹模型。
建立一個神經網路模型
建立一個森林模型。
建立一個漸變提升模型。
評估模型。
建立 ROC 圖表。

# Load necessary packages
library('swat')
library('ggplot2')
library('reshape2')
options(cas.print.messages = FALSE)

conn <- CAS('localhost', port=5570, 'YOUR SAS PROFILE EMAIL', 'YOUR SAS PROFILE PASSWORD', caslib = 'casuser')

actionsets <- c('sampling', 'fedsql', 'decisionTree', 'neuralNet', 'percentile')
for(i in actionsets){
    loadActionSet(conn, i)
}

data_dir <- '../data'
castbl <- cas.read.csv(conn, paste(data_dir, 'hmeq.csv', sep = '/'))

head(castbl)
summary(castbl)

# Bring data locally
df <- to.casDataFrame(castbl, obs = nrow(castbl))
# Use reshape2's melt to help with data formatting
d <- melt(df[sapply(df, is.numeric)], id.vars=NULL)
ggplot(d, aes(x = value)) +
    facet_wrap(~variable,scales = 'free_x') +
    geom_histogram(fill = 'blue', bins = 25)

tbl <- cas.simple.distinct(castbl)$Distinct[,c('Column', 'NMiss')]
tbl

cas.nmiss(castbl)

tbl$PctMiss <- tbl$NMiss/nrow(castbl)
ggplot(tbl, aes(Column, PctMiss)) +
    geom_col(fill = 'blue') +
    ggtitle('Pct Missing Values') +
    theme(plot.title = element_text(hjust = 0.5))

cas.dataPreprocess.impute(castbl,
    methodContinuous = 'MEDIAN',
    methodNominal = 'MODE',
    inputs = colnames(castbl)[-1],
    copyAllVars = TRUE,
    casOut = list(name = 'hmeq', 
                replace = TRUE)
)

cas.sampling.srs(conn,
    table = 'hmeq',
    samppct = 30,
    partind = TRUE,
    output = list(casOut = list(name = 'hmeq', replace = T), copyVars = 'ALL')
)

indata <- 'hmeq'
colinfo <- head(cas.table.columnInfo(conn, table = indata)$ColumnInfo, -1)
target <- colinfo$Column[1]
inputs <- colinfo$Column[-1]
nominals <- c(target, subset(colinfo, Type == 'varchar')$Column)
imp.inputs <- grep('IMP_', inputs, value = T)
imp.nominals <- c(target, grep('IMP_', nominals, value = T))

cas.decisionTree.dtreeTrain(conn,
    table = list(name = indata, where = '_PartInd_ = 0'),
    target = target,
    inputs = inputs,
    nominals = nominals,
    varImp = TRUE,
    casOut = list(name = 'dt_model', replace = TRUE)
)

cas.decisionTree.forestTrain(conn,
    table = list(name = indata, where = '_PartInd_ = 0'),
    target = target,
    inputs = inputs,
    nominals = nominals,
    casOut = list(name = 'rf_model', replace = TRUE)
)

cas.decisionTree.gbtreeTrain(conn,
    table = list(name = indata, where = '_PartInd_ = 0'),
    target = target,
    inputs = inputs,
    nominals = nominals,
    casOut = list(name = 'gbt_model', replace = TRUE)
)

cas.neuralNet.annTrain(conn,
    table = list(name = indata, where = '_PartInd_ = 0'),
    target = target,
    inputs = imp.inputs,
    hidden = 7,
    nominals = imp.nominals,
    casOut = list(name = 'nn_model', replace = TRUE)
)

models <- c('dt','rf','gbt','nn')
scores <- c(cas.decisionTree.dtreeScore, cas.decisionTree.forestScore, 
            cas.decisionTree.gbtreeScore, cas.neuralNet.annScore)
names(scores) <- models

# Function to help automate prediction process on new data
score.params <- function(model){return(list(
    object       = defCasTable(conn, indata),
    modelTable   = list(name = paste0(model, '_model')),
    copyVars     = list(target, '_PartInd_'),
    assessonerow = TRUE,
    casOut       = list(name = paste0(model, '_scored'), replace = T)
))}
lapply(models, function(x) {do.call(scores[[x]], score.params(x))})

loadActionSet(conn, 'percentile')
assess.model <- function(model){
    cas.percentile.assess(conn,
        table    = list(name = paste0(model,'_scored'), 
                        where = '_PartInd_ = 1'),
        inputs   = paste0('_', model, '_P_           1'),
        response = target,
        event    = '1')
}
model.names <- c('Decision Tree', 'Random Forest', 
                 'Gradient Boosting', 'Neural Network')
roc.df <- data.frame()
for (i in 1:length(models)){
    tmp <- (assess.model(models[i]))$ROCInfo
    tmp$Model <- model.names[i] 
    roc.df <- rbind(roc.df, tmp)
}
compare <- subset(roc.df, round(roc.df$CutOff, 2) == 0.5)
rownames(compare) <- NULL
compare[,c('Model','TP','FP','FN','TN')]

compare$Misclassification <- 1 - compare$ACC
miss <- compare[order(compare$Misclassification), c('Model','Misclassification')]
rownames(miss) <- NULL
miss

roc.df$Models <- paste(roc.df$Model, round(roc.df$C, 3), sep = ' - ')
ggplot(data = roc.df[c('FPR', 'Sensitivity', 'Models')],
    aes(x = as.numeric(FPR), y = as.numeric(Sensitivity), colour = Models)) +
    geom_line() +
    labs(x = 'False Positive Rate', y = 'True Positive Rate')

cas.session.endSession(conn)

開源程式碼 Python

這是一個簡單的端到端範例，說明如何使用SAS Viya 分析平台進行預測分析，主要有以下步驟：

載入所需的 Python 套件包。
在已執行的 CAS 伺服器上啟動 CAS 工作階段。
載入所需的 CAS 操作集。
將 HMEQ 資料檔從本地文件系統上傳至 CAS 伺服器中。
探索資料。
設算缺漏值。
將資料分為訓練資料和驗證資料。
建立一個決策樹模型。
建立一個森林模型。
建立一個神經網路模型
建立一個漸變提升模型。
評估模型。
建立 ROC 圖表。

from swat import *
from swat.render import render_html
from pprint import pprint
from matplotlib import pyplot as plt
import pandas as pd
import sys
%matplotlib inline

target          = "bad"
class_inputs    = ["reason", "job"]
class_vars      = [target] + class_inputs
interval_inputs = ["im_clage", "clno", "im_debtinc", "loan", "mortdue", "value", "im_yoj", "im_ninq", "derog", "im_delinq"]
all_inputs      = interval_inputs + class_inputs

indata_dir = '../data'
indata = 'hmeq'

cashost='localhost'
casport=5570
useremail='leoyeh.me@gmail.com'
userpassword='xxx'
casauth='~/.authinfo'
sess = CAS(cashost, casport, useremail, userpassword, caslib="casuser")
sess.loadactionset(actionset="dataStep")
sess.loadactionset(actionset="dataPreprocess")
sess.loadactionset(actionset="cardinality")
sess.loadactionset(actionset="sampling")
sess.loadactionset(actionset="regression")
sess.loadactionset(actionset="decisionTree")
sess.loadactionset(actionset="neuralNet")
sess.loadactionset(actionset="svm")
sess.loadactionset(actionset="astore")
sess.loadactionset(actionset="percentile")

if not sess.table.tableExists(table=indata).exists:
    indata = sess.upload_file(indata_dir+"/"+indata+".csv", casout={"name":indata})

sess.summary(indata)

sess.cardinality.summarize(
  table=indata, 
  cardinality={"name":"data_card", "replace":True}
)

tbl_data_card=sess.CASTable('data_card')
tbl_data_card.where='_NMISS_>0'
print("Data Summary".center(80, '-')) # print title
tbl_data_card.fetch() # print obs

tbl_data_card.vars=['_VARNAME_', '_NMISS_', '_NOBS_']
allRows=20000  # Assuming max rows in data_card table is <= 20,000
df_data_card=tbl_data_card.fetch(to=allRows)['Fetch']
df_data_card['PERCENT_MISSING']=(df_data_card['_NMISS_']/df_data_card['_NOBS_'])*100

tbl_forplot=pd.Series(list(df_data_card['PERCENT_MISSING']), index=list(df_data_card['_VARNAME_']))
ax=tbl_forplot.plot(
  kind='bar', 
  title='Percentage of Missing Values'
)
ax.set_ylabel('Percent Missing')
ax.set_xlabel('Variable Names');

r=sess.dataPreprocess.transform(
  table=indata,
  casOut={"name":"hmeq_prepped", "replace":True},
  copyAllVars=True,
  outVarsNameGlobalPrefix="IM",
  requestPackages=[
    {"impute":{"method":"MEAN"}, "inputs":{"clage"}},
    {"impute":{"method":"MEDIAN"}, "inputs":{"delinq"}},
    {"impute":{"method":"VALUE", "valuesNumeric":{2}}, "inputs":{"ninq"}},
    {"impute":{"method":"VALUE", "valuesNumeric":{35.0, 7, 2}}, "inputs":{"debtinc", "yoj"}}
  ]
)
render_html(r)

sess.sampling.stratified(
  table={"name":"hmeq_prepped", "groupBy":target},
  output={"casOut":{"name":"hmeq_part", "replace":True}, "copyVars":"ALL"},
  samppct=70,
  partind=True
)

sess.decisionTree.dtreeTrain(
  table={
    "name":"hmeq_part",
    "where":"strip(put(_partind_, best.))='1'"
  },
  inputs=all_inputs,
  target="bad",
  nominals=class_vars,
  crit="GAIN",
  prune=True,
  varImp=True,
  missing="USEINSEARCH",
  casOut={"name":"tree_model", "replace":True}
)

# Score 
sess.decisionTree.dtreeScore(
  table={"name":"hmeq_part"},
  modelTable={"name":"tree_model"},
  casOut={"name":"_scored_tree", "replace":True},
  copyVars={"bad", "_partind_"}
)

# Create p_bad0 and p_bad1 as _dt_predp_ is the probability of event in _dt_predname_
sess.dataStep.runCode(code = """
data _scored_tree; 
    length p_bad1 p_bad0 8.;
    set _scored_tree; 
    if _dt_predname_=1 then do; 
        p_bad1=_dt_predp_; 
        p_bad0=1-p_bad1; 
    end; 
    if _dt_predname_=0 then do; 
        p_bad0=_dt_predp_; 
        p_bad1=1-p_bad0; 
    end; 
run;
"""
)

sess.decisionTree.forestTrain(
  table={
    "name":"hmeq_part",
    "where":"strip(put(_partind_, best.))='1'"
  },
  inputs=all_inputs,
  nominals=class_vars,
  target="bad",
  nTree=50,
  nBins=20,
  leafSize=5,
  maxLevel=21,
  crit="GAINRATIO",
  varImp=True,
  missing="USEINSEARCH",
  vote="PROB",
  casOut={"name":"forest_model", "replace":True}
)

# Score 
sess.decisionTree.forestScore(
  table={"name":"hmeq_part"},
  modelTable={"name":"forest_model"},
  casOut={"name":"_scored_rf", "replace":True},
  copyVars={"bad", "_partind_"},
  vote="PROB"
)

# Create p_bad0 and p_bad1 as _rf_predp_ is the probability of event in _rf_predname_
sess.dataStep.runCode(
  code="data _scored_rf; set _scored_rf; if _rf_predname_=1 then do; p_bad1=_rf_predp_; p_bad0=1-p_bad1; end; if _rf_predname_=0 then do; p_bad0=_rf_predp_; p_bad1=1-p_bad0; end; run;"
)

sess.decisionTree.gbtreeTrain(
  table={
    "name":"hmeq_part",
    "where":"strip(put(_partind_, best.))='1'"
  },
  inputs=all_inputs,
  nominals=class_vars,
  target=target,
  nTree=10,
  nBins=20,
  maxLevel=6,
  varImp=True,
  missing="USEINSEARCH",
  casOut={"name":"gb_model", "replace":True}
)

# Score 
sess.decisionTree.gbtreeScore(
  table={"name":"hmeq_part"},
  modelTable={"name":"gb_model"},
  casOut={"name":"_scored_gb", "replace":True},
  copyVars={ target, "_partind_"}
)

# Create p_bad0 and p_bad1 as _gbt_predp_ is the probability of event in _gbt_predname_
sess.dataStep.runCode(
  code="data _scored_gb; set _scored_gb; if _gbt_predname_=1 then do; p_bad1=_gbt_predp_; p_bad0=1-p_bad1; end; if _gbt_predname_=0 then do; p_bad0=_gbt_predp_; p_bad1=1-p_bad0; end; run;"
)

sess.neuralNet.annTrain(
  table={
    "name":"hmeq_part",
    "where":"strip(put(_partind_, best.))='1'"
  },
  validTable={
    "name":"hmeq_part",
    "where":"strip(put(_partind_, best.))='0'"
  },
  inputs=all_inputs,
  nominals=class_vars,
  target="bad",
  hiddens={7},
  acts={"TANH"},
  combs={"LINEAR"},
  targetAct="SOFTMAX",
  errorFunc="ENTROPY",
  std="MIDRANGE",
  randDist="UNIFORM",
  scaleInit=1,
  nloOpts={
    "optmlOpt":{"maxIters":250, "fConv":1e-10}, 
    "lbfgsOpt":{"numCorrections":6},
    "printOpt":{"printLevel":"printDetail"},
    "validate":{"frequency":1}
  },
  casOut={"name":"nnet_model", "replace":True}
)

# Score 
sess.neuralNet.annScore(
  table={"name":"hmeq_part"},
  modelTable={"name":"nnet_model"},
  casOut={"name":"_scored_nn", "replace":True},
  copyVars={"bad", "_partind_"}
)

# Create p_bad0 and p_bad1 as _nn_predp_ is the probability of event in _nn_predname_
sess.dataStep.runCode(
  code="data _scored_nn; set _scored_nn; if _nn_predname_=1 then do; p_bad1=_nn_predp_; p_bad0=1-p_bad1; end; if _nn_predname_=0 then do; p_bad0=_nn_predp_; p_bad1=1-p_bad0; end; run;"
)

def assess_model(prefix):
    return sess.percentile.assess(
      table={
        "name":"_scored_" + prefix, 
        "where": "strip(put(_partind_, best.))='0'"
      },
      inputs=[{"name":"p_bad1"}],      
      response="bad",
      event="1",
      pVar={"p_bad0"},
      pEvent={"0"}      
    )

treeAssess=assess_model(prefix="tree")    
tree_fitstat =treeAssess.FitStat
tree_rocinfo =treeAssess.ROCInfo
tree_liftinfo=treeAssess.LIFTInfo

rfAssess=assess_model(prefix="rf")    
rf_fitstat =rfAssess.FitStat
rf_rocinfo =rfAssess.ROCInfo
rf_liftinfo=rfAssess.LIFTInfo

gbAssess=assess_model(prefix="gb")    
gb_fitstat =gbAssess.FitStat
gb_rocinfo =gbAssess.ROCInfo
gb_liftinfo=gbAssess.LIFTInfo

nnAssess=assess_model(prefix="nn")    
nn_fitstat =nnAssess.FitStat
nn_rocinfo =nnAssess.ROCInfo
nn_liftinfo=nnAssess.LIFTInfo

tree_liftinfo["model"]="DecisionTree"
tree_rocinfo["model"]="DecisionTree"
rf_liftinfo["model"]="Forest"
rf_rocinfo["model"]="Forest"
gb_liftinfo["model"]="GradientBoosting"
gb_rocinfo["model"]="GradientBoosting"
nn_liftinfo["model"]="NeuralNetwork"
nn_rocinfo["model"]="NeuralNetwork"

all_liftinfo=rf_liftinfo.append(gb_liftinfo, ignore_index=True).append(nn_liftinfo, ignore_index=True).append(tree_liftinfo, ignore_index=True)
all_rocinfo=rf_rocinfo.append(gb_rocinfo, ignore_index=True).append(nn_rocinfo, ignore_index=True).append(tree_rocinfo, ignore_index=True)
#all_liftinfo=rf_liftinfo.append(tree_liftinfo, ignore_index=True)
#all_rocinfo=rf_rocinfo.append(tree_rocinfo, ignore_index=True)

print("AUC (using validation data)".center(80, '-'))
all_rocinfo[["model", "C"]].drop_duplicates(keep="first").sort_values(by="C", ascending=False)

plt.figure()
for key, grp in all_rocinfo.groupby(["model"]):
    plt.plot(grp["FPR"], grp["Sensitivity"], label=key)
plt.plot([0,1], [0,1], "k--")
plt.xlabel("False Postivie Rate")
plt.ylabel("True Positive Rate")
plt.grid(True)
plt.legend(loc="best")
plt.title("ROC Curve (using validation data)")
plt.show()

plt.figure()
for key, grp in all_liftinfo.groupby(["model"]):
    plt.plot(grp["Depth"], grp["Lift"], label=key)
plt.xlabel("Depth")
plt.ylabel("Lift")
plt.grid(True)
plt.legend(loc="best")
plt.title("Lift Chart (using validation data)")
plt.show();

sess.close()

Previous第 7 章：人工智慧 Next第 9 章：模型管理

Last updated 4 years ago