【統計】尤度と最尤推定量
尤度推定のお勉強
【目次】
尤度と最尤推定量
最尤推定法は、標本から母集団のパラメータを点推定する方法（母集団の分布は既知）。
得られた標本を
（独立同一分布に従う）、推定したいパラメータを 
とすると、尤度関数は次と表される。
また、自然対数をとったものを対数尤度という。確率の足し算となり、値が小さくならないため計算機で扱いやすい。
この尤度関数を最大化するパラメータ を最尤推定量と呼ぶ。
なお、最小化問題などでは、対数尤度にマイナスを乗じて 0 を最小値とした負の対数尤度を用いることがある。
得られた標本を
また、自然対数をとったものを対数尤度という。確率の足し算となり、値が小さくならないため計算機で扱いやすい。
この尤度関数を最大化するパラメータ
なお、最小化問題などでは、対数尤度にマイナスを乗じて 0 を最小値とした負の対数尤度を用いることがある。
尤度と確率
尤度と確率では分布パラメータθと得られる値 x のどちらが変数か定数かが違う。
見た目が同じなのでややこしい。
■ 尤度
θで積分しても１にならない（１になるとは限らない）。比較に意味がある。
■ 確率
θで積分すると１になる。確率の合計は１。
見た目が同じなのでややこしい。
■ 尤度
θで積分しても１にならない（１になるとは限らない）。比較に意味がある。
■ 確率
θで積分すると１になる。確率の合計は１。
簡単な図解
尤度と最尤推定量のざっくりしたイメージ。
例としてパラメータθの候補に a, b, c の３つの値（変数）を用意し、各パラメータの下での尤度を計算。
この中で尤度が最大なのは b なので、最尤推定量は b になる。
参考
統計学（出版：東京図書）, 日本統計学会編
【数分解説】尤度(尤度関数): あるデータが与えられる時そのデータが出やすいパラメータを求める評価値が欲しい【Likelihood Function】 - YouTube
【統計学】尤度って何？をグラフィカルに説明してみる。 - Qiita
【統計】一般化線形混合モデル（GLMM）
一般化線形混合モデル (Generalized Linear Mixed Model, GLMM) の備忘録。
一般化線形モデル (Generalized Linear Model, GLM) との違いだったり、SAS や R での使い方だったり。
GLM と GLMM
GLMでは固定効果しか扱えなかったが、GLMMではランダム効果を扱えるようになった。
ランダム効果は次のようなもの。ランダム効果を投入することで、測定データ間の関係を細かく検討できるようになった。
- 被験者によるズレ（同じ痛みでも人により反応が違う）
- 評価者によるズレ（甘い採点／辛い採点）
- 施設・地域によるズレ（組織・地域的な性格・性質）
GLM
■ モデル、期待値・分散
GLMM
■ モデル、期待値・分散
式をみると、GLMの分散にランダム効果に関わる項（Z' G Z）が加わった形となっている。
SAS (proc mixed)
SASの mixed プロシジャでの指定の仕方。固定効果、ランダム効果、反復効果をそれぞれ Model、Random、Repeatedステートメントで指定する。
モデルにランダム効果や反復効果が要らなければ、そのステートメント書かなくて良い。
ランダム効果と反復効果の分散共分散構造を指定したい場合、各ステートメントのオプションで指定できる。詳細はプロシジャガイドや参考資料を参照。
オプションに g , r を指定すると、それぞれ係数γの分散共分散行列Gと誤差εの分散共分散行列Rがアウトプットされる。
ランダム効果として、ランダム切片を入れたければ intercept と書く。ランダム傾きを入れたければ任意の変数を指定する。
推定方法のデフォルトは制限付き最尤法 (REstricted Maximum Likelihood, REML) で、最尤法 (Maximum Likelihood, ML) への切り替えもできる。
R (lmer, nlme)
r の lmer パッケージでランダム効果を加えたモデルを作成できる。
library(tidyverse) library(lme4) library(lmerTest) #-- data ADS <- sleepstudy #-- model fit (mod1 <- lmer(Reaction ~ Days + (1 | Subject), ADS)) #-- ランダム切片 (mod2 <- lmer(Reaction ~ Days + (1 + Days | Subject), ADS)) #-- ランダム切片＋ランダム傾き (mod2 <- lmer(Reaction ~ Days + (Days | Subject), ADS)) #-- 別記法 #-- random effect lme4::ranef(mod1) #-- fixed effect lme4::fixef(mod1) summary(mod1)
nlme パッケージでも対応できるみたい。ここは未勉強なので使い方の詳細はよく分からない。
参考
SAS Help Center
http://www.sigmath.es.osaka-u.ac.jp/~kano/research/application/gasshuku02/sas1/MIXED.pdf
https://content.sph.harvard.edu/fitzmaur/ala/lectures.pdf
https://www.mwsug.org/proceedings/2007/stats/MWSUG-2007-S02.pdf
Introduction to SAS® Proc Mixed - ppt download
統計手法アラカルト Mixed Model ～混合モデル～ - ppt download
第4章 MIXED Model 4.1 MIXED Model とは 4.2 反復測定データの分析１ 分割法タイプのデータ - ppt download
【SAS】SASで実装できる機械学習
まとめ用。随時更新。
教師あり学習
最小二乗回帰（OLS：Ordinary Least Squares regression、線形回帰）
【プロシジャ】HPREG, GLMSELECT, HPGENSELECT(dist=gaussian)
【SAS】glmselectプロシジャ_Lasso回帰、ElasticNet回帰 - こちにぃるの日記
(例：10-fold ElasticNet回帰)
ods graphics on;
proc glmselect data=訓練データ plots=all seed=19;
model ラベル = 特徴量1, ..., 特徴量p
/ details=all stats=all
selection=elasticnet(choose=cv l2search=grid showl2search stop=none) cvmethod=random(10) stb;
score out=score_elasticnet predicted residual;
code file="ファイル名.sas";
run;
ods graphics off;
*-- 予測;
data _Pred;
set 予測データ;
%inc "ファイル名.sas"; *-- 新しいデータで予測;
run;
ロジスティック回帰（Logistic regression）
【プロシジャ】HPLOGISTIC, HPGENSELECT(dist=binary link=logit)
【SAS】hpgenselect プロシジャ_罰則付きロジスティック回帰（LASSO） - こちにぃるの日記
(例：Logistic regression with Lasso)
proc hpgenselect data=訓練データ lassorho=0.8 lassosteps=20;
partition roleVar=selected(train='0' validate='1');
model ラベル(event="1")=特徴量1, ..., 特徴量p / distribution=binary link=logit;
selection method=lasso(choose=validate stop=none) details=all;
performance details;
code file="ファイル名.sas" residual; *-- 予測用のPGM生成;
ods output ParameterEstimates=_Coef_lasso; *-- 回帰係数;
output data=_Pred_lasso role=ROLE pred=PRED resid=RESID;
run;
*-- 予測;
data _Pred;
set 予測データ;
%inc "ファイル名.sas"; *-- 新しいデータで予測;
run;
四分位回帰（Quantile regression）
【プロシジャ】HPQUANTSELECT
サポートベクターマシーン（SVM: Support vector machine）
【プロシジャ】HPSVM
【SAS】hpsvmプロシジャ_サポートベクターマシーン - こちにぃるの日記
proc hpsvm data = 訓練データ;
target ラベル / order = 順序 level = 変数のタイプ;
input 名義尺度の特徴量 / level=nominal;
input 順序尺度の特徴量 / level=ordinal order=ascending;
input 間隔尺度の特徴量 / level=interval;
code file = "ファイル名.sas" ; /* 将来データ用に予測アルゴリズム出力 */
run;
*-- 予測 ;
data _Pred;
set 予測データ;
%inc "ファイル名.sas"; *-- 新しいデータで予測;
run;
決定木（Decision tree）
【プロシジャ】HPSPLIT
proc hpsplit data=train;
class ラベル カテゴリの特徴量 ;
model ラベル(event="1") = 特徴量1, ..., 特徴量p ;
partition roleVar=selected(valid='0' train='1');
prune costcomplexity;
code file= "ファイル名.sas" ; /* 将来データ用に予測アルゴリズム出力 */
output out=_SCORE;
run;
*-- 予測 ;
data _Pred;
set 予測データ;
%inc "ファイル名.sas"; *-- 新しいデータで予測;
run;
ランダムフォレスト（Random Forest）
【プロシジャ】HPFOREST
【SAS】hpforestプロシジャ_ランダムフォレスト - こちにぃるの日記
proc hpforest data = 訓練データ seed = 19 scoreprole = valid
maxtrees = 100 /* pythonではn_estimators */
maxdepth = 6 /* pythonではmax_depth */
;
partition roleVar=selected(train='0' valid='1');
performance details;
target ラベル / order = 順序 level = 変数のタイプ;
input 名義尺度の特徴量 / level=nominal;
input 順序尺度の特徴量 / level=ordinal order=ascending;
input 間隔尺度の特徴量 / level=interval;
ods output FitStatistics = _Res_FitStatistics; /* 予測値の出力 */
ods output VariableImportance = _Res_VariableImportance; /* 変数重要度（特徴量重要度）の出力 */
save file = "hpforest_fit.bin"; /* 将来データ用に予測アルゴリズム出力 */
run;
*-- 予測 ;
proc hp4score data = 予測データ;
score file="ファイル名.bin" out=_Pred_RandomForest;
id 表示したい特徴量;
run;
勾配ブースティング木（Gradient Boosting Tree）
【プロシジャ】TREEBOOST
【SAS】treeboostプロシジャ_Gradient Boosting Tree（勾配ブースティング木） - こちにぃるの日記
proc treeboost data=訓練データ(where=(selected=0))
iterations = 1000 /* pythonではn_estimators */
maxdepth = 6 /* pythonではmax_depth */
seed = 19
;
target ラベル / order = 順序 level = 変数のタイプ;
input 名義尺度の特徴量 / level=nominal;
input 順序尺度の特徴量 / level=ordinal order=ascending;
input 間隔尺度の特徴量 / level=interval;
assess validata = valid_dataset;
subseries best;
code file="ファイル名.sas";
save fit=_gb_fit importance=_gb_importance; /* 予測値と特徴量重要度 */
quit;
run;
*-- 予測;
data _Pred;
set 予測データ;
%inc "ファイル名.sas"; *-- 新しいデータで予測;
run;
k近傍法（k-NN: k-nearest neighbor algorithm）
【プロシジャ】
教師なし学習
主成分分析（PCA: principal component analysis）
【プロシジャ】
k-means法
【プロシジャ】
【HTML】別HTMLを埋め込む（iframe）
作成したHTMLテーブルを別のHTMLファイルに埋め込んで表示させたい！と思ったので、その方法をメモ。
HTMLの埋め込みは、iframe タグで実装できる。
<iframe src='パス/ファイル名.html'></iframe>
試しに前回記事で作ったHTMLテーブルを埋め込んでみる。
【VBA & HTML】ExcelのテーブルをHTML形式に変換する - こちにぃるの日記
HTMLは以下。
<!DOCTYPE html> <html lang='ja'> <head> <meta charset='utf8'> <title>Top</title> <link rel='stylesheet' href='https://stackpath.bootstrapcdn.com/bootstrap/4.3.1/css/bootstrap.min.css' integrity='sha384-ggOyR0iXCbMQv3Xipma34MD+dH/1fQ784/j6cY/iJTQUOhcWr7x9JvoRxT2MZw1T' crossorigin='anonymous'> <script type='text/javascript' src='https://ajax.googleapis.com/ajax/libs/jquery/3.6.0/jquery.min.js'></script> <style type='text/css'> #header{ background-color: #CBFFD3; position: fixed; top: 0px; left: 0px; width: 100%; height: 75px; border: 5px solid #fff; } #main_body{ padding:75px 0 75px; top: 75px; } .side { background-color:#FFFFEE; width:10%; overflow-y: scroll; border: 5px solid #fff; } .fixed{ position:fixed; top:75px; width: 300px; } .content { width:90%; margin-left:20px; } article { display:flex; } </style> <script type='text/javascript'> $(function($) { var tab = $('.sidebar'), offset = tab.offset(); $(window).scroll(function () { if($(window).scrollTop() > offset.top) { tab.addClass('fixed'); } else { tab.removeClass('fixed'); } }); }); </script> </head> <body> <div id='header' style='font-size:30px;'>タイトルXXX</div> <article id='main_body'> <div class='side'> <div class='sidebar'> <p>目次</p> <ol> <li><a href='#Table1'>テーブル１</a></li> <li><a href='#Table2'>テーブル２</a></li> <li><a href='#Table3'>テーブル３</a></li> </ol> </div> </div> <div class='content'> <hr/> <div id='Table1'>■テーブル１</div> <iframe width='90%' height='500px' src='./SampleHTMLtable.html' frameborder='0'></iframe> <br/> <hr/> <div id='Table2'>■テーブル２</div> <iframe width='40%' height='200px' src='./SampleHTMLtable.html' frameborder='0'></iframe> <br/> <hr/> <div id='Table3'>■テーブル３</div> <iframe width='70%' height='500px' src='./SampleHTMLtable.html' frameborder='0'></iframe> <br/> </div> </article> </body> </html>
うまくHTMLテーブルを表示できた。
サイドメニューの設定は以下サイト様を参考にさせていただきました。
直帰率改善に！スクロールしてもサイドバーを固定して記事一覧を表示する方法 | イズクル
【VBA & HTML】ExcelのテーブルをHTML形式に変換する
Excelを開いて表を見るのは面倒だったりする。
（Excelの起動を待ったり、別作業でExcelを使っていたりするとき）
そんなときブラウザで表形式データが見れたら便利だ、というモチベーションでExcelで作成したテーブルをHTMLに変換するプログラムを作ってみた。
Excel の準備
こんな感じでExcelを準備。セル色：グレーに設定情報、白に実データが入ります。
* 3行目：HTMLタイトル入力。
* 5行目：Table対象にする列に「1」を設定し、該当列のみHTMLに出力。
* 6行目：aタグ用にURLを入力した列の番号を指定。例だと7列目にPDFのURLを指定。
* 7行目：imgタグ用にURLを入力した列の番号を指定。例だと8列目にjpegファイルのURLを指定。
* 8行目：当該列の列幅を指定。
* 9行目：ヘッダー名を入力。
* 10行目以降：実データを入力。
VBAコード
このエクセルにHTML生成用のVBAコードを仕込む。
Option Explicit Sub CreateHTMLTable() '--------------------------------------- ' 設定 '--------------------------------------- On Error GoTo ERROR1 Application.DisplayAlerts = False Application.ScreenUpdating = False Dim myWB As Workbook, myWS As Worksheet Dim strow, enrow, stcol, encol, t_row, a_row, i_row, w_row, h_row, r, c As Long Dim outfile As String Dim adoSt As Object Set adoSt = CreateObject("ADODB.Stream") Const adTypeBinary = 1 Const adTypeText = 2 Const adCR = 13 Const adLF = 10 Const adCRLF = -1 Const adWriteChar = 0 Const adWriteLine = 1 Const adSaveCreateNotExist = 1 Const adSaveCreateOverWrite = 2 Set myWB = ThisWorkbook Set myWS = myWB.Sheets(ActiveSheet.Name) t_row = 5 a_row = 6 i_row = 7 w_row = 8 h_row = 9 strow = 10 enrow = myWS.UsedRange.Find("*", , xlFormulas, , xlByRows, xlPrevious).Row stcol = 2 encol = myWS.UsedRange.Find("*", , xlFormulas, , xlByColumns, xlPrevious).Column '--------------------------------------- ' HTMLテーブル生成 '--------------------------------------- outfile = Replace(myWS.Cells(3, 2), "\", "_") outfile = Replace(outfile, "/", "_") outfile = Replace(outfile, ":", "_") outfile = Replace(outfile, "*", "_") outfile = Replace(outfile, "?", "_") outfile = Replace(outfile, """", "_") outfile = Replace(outfile, "<", "_") outfile = Replace(outfile, ">", "_") outfile = Replace(outfile, "|", "_") outfile = Replace(outfile, " ", "_") outfile = Replace(outfile, " ", "_") With adoSt .Charset = "UTF-8" .Type = adTypeText .LineSeparator = adLF .Open .WriteText "<!DOCTYPE html>", adWriteLine .WriteText "<html lang='ja'>", adWriteLine .WriteText "<head>", adWriteLine .WriteText " <meta charset='utf8'>", adWriteLine .WriteText " <title>SampleHTMLtable</title>", adWriteLine .WriteText " <!-- Bootstrap -->", adWriteLine .WriteText " <link rel='stylesheet' href='https://stackpath.bootstrapcdn.com/bootstrap/4.3.1/css/bootstrap.min.css' integrity='sha384-ggOyR0iXCbMQv3Xipma34MD+dH/1fQ784/j6cY/iJTQUOhcWr7x9JvoRxT2MZw1T' crossorigin='anonymous'>", adWriteLine .WriteText " ", adWriteLine .WriteText " <!-- jQuery -->", adWriteLine .WriteText " <script type='text/javascript' src='https://ajax.googleapis.com/ajax/libs/jquery/3.6.0/jquery.min.js'></script>", adWriteLine .WriteText " ", adWriteLine .WriteText " <!-- Datatables -->", adWriteLine .WriteText " <link rel='stylesheet' href='https://cdn.datatables.net/t/bs-3.3.6/jqc-1.12.0,dt-1.10.11/datatables.min.css'/> ", adWriteLine .WriteText " <script src='https://cdn.datatables.net/t/bs-3.3.6/jqc-1.12.0,dt-1.10.11/datatables.min.js'></script>", adWriteLine .WriteText " <script type='text/javascript'>", adWriteLine .WriteText " jQuery(function($){", adWriteLine .WriteText " $.extend( $.fn.dataTable.defaults, {", adWriteLine .WriteText " language: {", adWriteLine .WriteText " url: 'http://cdn.datatables.net/plug-ins/9dcbecd42ad/i18n/Japanese.json'", adWriteLine .WriteText " } ", adWriteLine .WriteText " }); ", adWriteLine .WriteText " $('#TableLayout1').DataTable({", adWriteLine .WriteText " scrollX: true", adWriteLine .WriteText " ,scrollY: '75vh'", adWriteLine .WriteText " ,displayLength: 10", adWriteLine .WriteText " });", adWriteLine .WriteText " });", adWriteLine .WriteText " </script>", adWriteLine .WriteText "</head>", adWriteLine .WriteText "<body>", adWriteLine .WriteText " <table id='TableLayout1' border='1'>", adWriteLine .WriteText " <thead style='background-color: #64F9C1;'>", adWriteLine .WriteText " <tr>", adWriteLine For c = stcol To encol If myWS.Cells(t_row, c) = 1 Then .WriteText " <th width='" & myWS.Cells(w_row, c) & "px'>" & EscapeTxt(myWS.Cells(h_row, c)) & "</th>", adWriteLine End If Next c .WriteText " </tr>", adWriteLine .WriteText " </thead>", adWriteLine .WriteText " <tbody>", adWriteLine For r = strow To enrow .WriteText " <tr>", adWriteLine For c = stcol To encol If myWS.Cells(t_row, c) = 1 Then If myWS.Cells(a_row, c) <> "" Then If myWS.Cells(r, myWS.Cells(a_row, c)) <> "" Then .WriteText " <td><a target='_blank' href='" & myWS.Cells(r, myWS.Cells(a_row, c)) & "'>" & EscapeTxt(myWS.Cells(r, c)) & "</a></td>", adWriteLine Else .WriteText " <td>" & EscapeTxt(myWS.Cells(r, c)) & "</td>", adWriteLine End If ElseIf myWS.Cells(i_row, c) <> "" Then If myWS.Cells(r, myWS.Cells(i_row, c)) <> "" Then .WriteText " <td><img src='" & myWS.Cells(r, myWS.Cells(i_row, c)) & "' title='" & EscapeTxt(myWS.Cells(r, c)) & "' width='" & myWS.Cells(w_row, c) & "px'></td>", adWriteLine Else .WriteText " <td>" & EscapeTxt(myWS.Cells(r, c)) & "</td>", adWriteLine End If Else .WriteText " <td>" & EscapeTxt(myWS.Cells(r, c)) & "</td>", adWriteLine End If End If Next c .WriteText " </tr>", adWriteLine Next r .WriteText " </tbody>", adWriteLine .WriteText " </table>", adWriteLine .WriteText "</body>", adWriteLine .WriteText "</html>", adWriteLine .SaveToFile ThisWorkbook.Path & "\" & outfile & ".html", adSaveCreateOverWrite .Close End With myWS.Cells(1, 1).Select MsgBox "出力しました" GoTo END1 '--------------------------------------- ' 後始末 '--------------------------------------- ERROR1: MsgBox Err.Number & ":" & Err.Description END1: Application.DisplayAlerts = False Application.ScreenUpdating = False Set adoSt = Nothing Set myWB = Nothing Set myWS = Nothing End Sub Function EscapeTxt(InTxt As String) As String EscapeTxt = Replace(InTxt, "&", "&") EscapeTxt = Replace(EscapeTxt, "<", "<") EscapeTxt = Replace(EscapeTxt, ">", ">") EscapeTxt = Replace(EscapeTxt, "'", "'") EscapeTxt = Replace(EscapeTxt, """", """) EscapeTxt = Replace(EscapeTxt, " ", " ") EscapeTxt = Replace(EscapeTxt, vbNewLine, "<br/>") EscapeTxt = Replace(EscapeTxt, vbCrLf, "<br/>") EscapeTxt = Replace(EscapeTxt, vbCr, "<br/>") EscapeTxt = Replace(EscapeTxt, vbLf, "<br/>") End Function
HTMLは、Stream オブジェクト（ ADODB.Stream ）を使って UTF 形式で出力している。
Stream オブジェクト (ADO) - ActiveX Data Objects (ADO) | Microsoft Docs
Stream オブジェクトのプロパティ、メソッド、およびイベント - ActiveX Data Objects (ADO) | Microsoft Docs
また、HTMLテーブルにソート・フィルター機能等を付けたかったので DataTables を使用している。
DataTables を初めて使ったが、実装は簡単で見栄えも良い。表の幅とか高さとかは使う環境で微調整が必要。
DataTables | Table plug-in for jQuery
Options
HTMLファイルを生成
表に値を入れて、HTMLファイルを生成してみるとこんな感じ。
PDFのリンクも開け、画像も表示され、最低限の機能は揃えられたかな？
一元管理ができそうで、何より何より。
今回作ったHTMLテーブルを別のHTMLファイルに埋め込む方法は以下の記事で。
【HTML】別HTMLを埋め込む（iframe） - こちにぃるの日記
【VBA】パラメタを変えてスクリプトを大量生成
久々にVBA。
R、Python、SAS ... なんでも良いが、
共有ロジックを外部ファンクションや外部マクロにしておいて、
パラメタを変えてぐるぐる回すことがある。
例えば、疾患名を変えて患者数を推移図にするなど。
そんなとき、スクリプトをいちいち作っていると面倒だし、ミスにもつながる。
だからVBAで一括生成してしまおうというお話。
こんな感じでエクセルに情報を埋めて、
VBAコードを仕込んで実行すると、
Sub OutputScript() '--------------------------------------- ' 設定 '--------------------------------------- On Error GoTo ERROR1 Application.DisplayAlerts = False Application.ScreenUpdating = False Dim myWB As Workbook, myWS As Worksheet Dim strow, enrow, encol, fnumber, i As Long Dim outpath As String Set myWB = ThisWorkbook Set myWS = myWB.Sheets(ActiveSheet.Name) strow = 7 enrow = myWS.UsedRange.Find("*", , xlFormulas, , xlByRows, xlPrevious).Row outpath = myWS.Cells(3, 3) If outpath = "" Then outpath = myWB.Path '--------------------------------------- ' テキスト生成 '--------------------------------------- For i = strow To enrow fnumber = FreeFile Open outpath & "\" & myWS.Cells(i, 2) For Output As #fnumber Print #fnumber, "#===============================================================================" Print #fnumber, "# ファイル名 : " & myWS.Cells(i, 2) Print #fnumber, "# 作 成 日 : " & Format(Now(), "yyyy/mm/dd") Print #fnumber, "# 作 成 者 : " & myWS.Cells(4, 3) Print #fnumber, "#===============================================================================" Print #fnumber, "#-------------------------------------------------------------------------------" Print #fnumber, "# 設定" Print #fnumber, "#-------------------------------------------------------------------------------" Print #fnumber, "library(tidyverse)" Print #fnumber, "setwd(パスを入力)" Print #fnumber, "" Print #fnumber, "#-------------------" Print #fnumber, "# ロード" Print #fnumber, "#-------------------" Print #fnumber, "Source('./myFunctionX.r')" Print #fnumber, "" Print #fnumber, "#-------------------" Print #fnumber, "# 実行" Print #fnumber, "#-------------------" Print #fnumber, "myFunctionX(PARAM1='" & myWS.Cells(i, 3) & "', PARAM2=" & myWS.Cells(i, 4) & ")" Print #fnumber, "" Print #fnumber, "#-------------------------------------------------------------------------------" Print #fnumber, "# End of File" Print #fnumber, "#-------------------------------------------------------------------------------" Close #fnumber Next i myWS.Cells(1, 1).Select MsgBox "出力しました" GoTo END1 '--------------------------------------- ' 後始末 '--------------------------------------- ERROR1: MsgBox Err.Number & ":" & Err.Description END1: Application.DisplayAlerts = False Application.ScreenUpdating = False Set myWB = Nothing Set myWS = Nothing End Sub
一気にスクリプトが生成される（今回はRスクリプト）。
スクリプトの中身はこんな感じ。外部ファンクションへの指定パラメタをスクリプト毎に変えて生成している。
これで少し楽ができる。めでたし、めでたし。
【Python】Kerasでロジスティック回帰
Kerasの勉強用にロジスティック回帰と同じ構造を作って実行してみる。
【目次】
基本の整理
ニューラルネットワークは入力層ー中間層ー出力層からなり、図示するとこんな感じ。
このうち、ロジスティック回帰は中間層がない単純な形。
活性化関数にはシグモイド関数が使われる。
また、損失関数には交差エントロピー誤差が使われる。
使用データ
import numpy as np import pandas as pd from scipy import stats from sklearn import model_selection #-- import train = pd.read_csv("../input/titanic/train.csv") test = pd.read_csv("../input/titanic/test.csv") ###--- 欠損処理 train["Age"] = train["Age"].fillna(train["Age"].median()) #--中央値で補完 train["Embarked"] = train["Embarked"].fillna(train["Embarked"].mode().iloc(0)) #-- 最頻値で補完 test["Age"] = test["Age"].fillna(test["Age"].median()) #--中央値で補完 test["Fare"] = test["Fare"].fillna(test["Fare"].median()) #--中央値で補完 ###--- カテゴリの数値化 train["Sex"] = train["Sex"].map({'male':0, 'female':1}) train = pd.get_dummies(train, columns=['Embarked']) test["Sex"] = test["Sex"].map({'male':0, 'female':1}) test = pd.get_dummies(test, columns=['Embarked']) ###--- データ分割 Xval = ["Pclass", "Sex", "Age", "Fare", "Embarked_C", "Embarked_Q", "Embarked_S"] df_y = train["Survived"].values df_X = train[Xval].values train_X, valid_X, train_y, valid_y = model_selection.train_test_split(df_X, df_y, random_state=19) test_X = test[Xval].values
Keras での試し書き（モデル作成）
先ず下準備として乱数シードを固定する。
また、データを標準化した方が精度が良かったので標準化。
import pandas as pd import numpy as np import tensorflow as tf from keras.layers import Activation, Dense, Dropout, BatchNormalization from keras.models import Sequential from keras import optimizers, callbacks from keras.utils.np_utils import to_categorical from sklearn.preprocessing import StandardScaler #-- 乱数シード np.random.seed(seed=19) tf.random.set_seed(19) #-- 標準化 sc = StandardScaler() train_X_std = sc.fit_transform(train_X) valid_X_std = sc.transform(valid_X)
モデルを作成する。
Sequential() でインスタンス化し、add メゾッドでモデルを定義する。
units がユニット数。ロジスティック回帰だと層が一つなのでここが出力層の定義になる。
input_shape が入力層で特徴量の数。
activation が活性化関数。ロジスティック回帰だとシグモイド関数。
#-- モデル model = Sequential() model.add(Dense(units=1 ,input_shape=(len(Xval),) ,activation="sigmoid"))
作成モデルをコンパイルする。
#-- コンパイル model.compile(optimizer=optimizers.SGD(lr=0.1) ,loss='binary_crossentropy' ,metrics=["accuracy"])
コンパイルしたモデルを最適化する。
callbacks.EarlyStopping で Early Stopping を実装できる。
callbacks.ModelCheckpoint('best_model.h5', save_best_only=True) でベストモデルを保存できる（後にモデルの読み込み方法を記載）。
#-- Fitting callbacks = [ callbacks.ModelCheckpoint('best_model.h5', save_best_only=True) ,callbacks.EarlyStopping(patience=10) ] history = model.fit(train_X_std, train_y ,batch_size=64 ,epochs=1000 ,verbose=1 ,callbacks=callbacks ,validation_data=(valid_X_std, valid_y))
Keras での試し書き（バリデーションデータで評価）
作成モデルで予測する。
#-- 予測 pred_y = model.predict(valid_X_std) pred_y = np.where(pred_y < 0.5, 0, 1)
重み係数を見たい場合、次で取得できる。
model.layers[0].get_weights()
今回の作成モデルを評価してみる。
ついでにsklearnで作ったロジスティック回帰を並べてみる。
from sklearn.linear_model import LogisticRegression model2 = LogisticRegression(solver='liblinear', random_state=19) model2.fit(train_X_std, train_y) pred_y2 = model2.predict(valid_X_std) pred_y2 = np.where(pred_y2 < 0.5, 0, 1)
from sklearn.metrics import classification_report, confusion_matrix print("-- keras --") print(confusion_matrix(valid_y, pred_y)) print(classification_report(valid_y, pred_y)) print() print("-- sklearn.linear_model.LogisticRegression --") print(confusion_matrix(valid_y, pred_y2)) print(classification_report(valid_y, pred_y2))
アルゴリズムが違うだろうが、だいたい同じ？
Keras での試し書き（テストデータで予測）
別セッションで作成モデルを使用したい場合、最初からPGMを再実行するのは手間。
モデルを保存してあれば load_model(パス/ファイル名) で読み込んで再利用できる。
from keras.models import load_model model3 = load_model('./best_model.h5') pred_y3 = model3.predict(test_X) pred_y3 = np.where(pred_y3 < 0.5, 0, 1)
参考
Home - Keras Documentation
Kerasによる2クラスロジスティック回帰 - 人工知能に関する断創録
本ブログで参考になりそうなもの
【統計】ロジスティック回帰分析 - こちにぃるの日記