【統計】尤度と最尤推定量

尤度推定のお勉強

【目次】


尤度と最尤推定


最尤推定法は、標本から母集団のパラメータを点推定する方法(母集団の分布は既知)。

得られた標本を  x_{i} = x_{1}, x_{2},...,x_{n} (独立同一分布に従う)、推定したいパラメータを \theta とすると、尤度関数は次と表される。

 \quad L\left( \theta | x \right) = \prod _{i=1}^{n} f\left( x_{i} | \theta \right)


また、自然対数をとったものを対数尤度という。確率の足し算となり、値が小さくならないため計算機で扱いやすい。

 \quad \ln L\left( \theta | x \right) = \ln \prod _{i=1}^{n} f\left( x_{i} | \theta \right) = \Sigma _{i=1}^{n} \ln f\left( x_{i} | \theta \right)


この尤度関数を最大化するパラメータ \theta最尤推定量と呼ぶ。

 \quad \widehat{\theta }_{ML}=\arg \max f\left( x| \theta \right)


なお、最小化問題などでは、対数尤度にマイナスを乗じて 0 を最小値とした負の対数尤度を用いることがある。

尤度と確率


尤度と確率では分布パラメータθと得られる値 x のどちらが変数か定数かが違う。
見た目が同じなのでややこしい。

■ 尤度

 \quad L\left( \theta | x \right) = \prod _{i=1}^{n} \Pr \left( x | \theta \right)

得られた値 x(定数)を元にして、分布パラメータθ(変数)を求める。
θで積分しても1にならない(1になるとは限らない)。比較に意味がある。

■ 確率

 \quad \prod _{i=1}^{n} \Pr \left( x | \theta \right)

分布パラメータθ(定数)の下で得られる値 x(変数)の確率を求める。
θで積分すると1になる。確率の合計は1。

簡単な図解


尤度と最尤推定量のざっくりしたイメージ。

f:id:cochineal19:20211107235358p:plain

例としてパラメータθの候補に a, b, c の3つの値(変数)を用意し、各パラメータの下での尤度を計算。
この中で尤度が最大なのは b なので、最尤推定量は b になる。

参考


統計学(出版:東京図書), 日本統計学会編
【数分解説】尤度(尤度関数): あるデータが与えられる時そのデータが出やすいパラメータを求める評価値が欲しい【Likelihood Function】 - YouTube
【統計学】尤度って何?をグラフィカルに説明してみる。 - Qiita

【統計】一般化線形混合モデル(GLMM)

一般化線形混合モデル (Generalized Linear Mixed Model, GLMM) の備忘録。
一般化線形モデル (Generalized Linear Model, GLM) との違いだったり、SAS や R での使い方だったり。


GLM と GLMM


GLMでは固定効果しか扱えなかったが、GLMMではランダム効果を扱えるようになった。
ランダム効果は次のようなもの。ランダム効果を投入することで、測定データ間の関係を細かく検討できるようになった。

  • 被験者によるズレ(同じ痛みでも人により反応が違う)
  • 評価者によるズレ(甘い採点/辛い採点)
  • 施設・地域によるズレ(組織・地域的な性格・性質)

GLM


■ モデル、期待値・分散

 \quad \bf{Y} = \bf{X' \beta} + \bf{\varepsilon}

 \quad E \left( \bf{Y} \right) = \bf{X' \beta}

 \quad V \left( \bf{Y} \right) = V \left( \bf{\varepsilon} \right) = \bf{R}

GLMM


■ モデル、期待値・分散

 \quad \bf{Y} = \bf{X' \beta} + \bf{Z' \gamma} + \bf{\varepsilon}

 \quad E \left( \bf{Y} \right) = \bf{X' \beta}

 \quad V \left( \bf{Y} \right) = \bf{Z \ V \left( \bf{\gamma} \right) \ Z' } + V \left( \bf{\varepsilon} \right) = \bf{Z \ G \ Z' } + \bf{R}


式をみると、GLMの分散にランダム効果に関わる項(Z' G Z)が加わった形となっている。

SAS (proc mixed)


SASの mixed プロシジャでの指定の仕方。固定効果、ランダム効果、反復効果をそれぞれ Model、Random、Repeatedステートメントで指定する。

f:id:cochineal19:20211010162011p:plain:w350

f:id:cochineal19:20211010164603p:plain:w600

モデルにランダム効果や反復効果が要らなければ、そのステートメント書かなくて良い。

ランダム効果と反復効果の分散共分散構造を指定したい場合、各ステートメントのオプションで指定できる。詳細はプロシジャガイドや参考資料を参照。
オプションに 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 反復測定データの分析1 分割法タイプのデータ - 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'>テーブル1</a></li>
          <li><a href='#Table2'>テーブル2</a></li>
          <li><a href='#Table3'>テーブル3</a></li>
        </ol>
      </div> 
    </div> 
    
    <div class='content'>
      <hr/>
      <div id='Table1'>■テーブル1</div>
      <iframe width='90%' height='500px' src='./SampleHTMLtable.html' frameborder='0'></iframe>
      <br/>
      
      <hr/>
      <div id='Table2'>■テーブル2</div>
      <iframe width='40%' height='200px' src='./SampleHTMLtable.html' frameborder='0'></iframe>
      <br/>
      
      <hr/>
      <div id='Table3'>■テーブル3</div>
      <iframe width='70%' height='500px' src='./SampleHTMLtable.html' frameborder='0'></iframe>
      <br/>
      
    </div>
    
  </article>
  
</body>
</html>


うまくHTMLテーブルを表示できた。 f:id:cochineal19:20210829160959p:plain

サイドメニューの設定は以下サイト様を参考にさせていただきました。
直帰率改善に!スクロールしてもサイドバーを固定して記事一覧を表示する方法 | イズクル

【VBA & HTML】ExcelのテーブルをHTML形式に変換する

Excelを開いて表を見るのは面倒だったりする。
Excelの起動を待ったり、別作業でExcelを使っていたりするとき)

そんなときブラウザで表形式データが見れたら便利だ、というモチベーションでExcelで作成したテーブルをHTMLに変換するプログラムを作ってみた。

Excel の準備


こんな感じでExcelを準備。セル色:グレーに設定情報、白に実データが入ります。
f:id:cochineal19:20210829152628p:plain
* 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, "&", "&amp;")
    EscapeTxt = Replace(EscapeTxt, "<", "&lt;")
    EscapeTxt = Replace(EscapeTxt, ">", "&gt;")
    EscapeTxt = Replace(EscapeTxt, "'", "&#39;")
    EscapeTxt = Replace(EscapeTxt, """", "&quot;")
    EscapeTxt = Replace(EscapeTxt, " ", "&nbsp;")
    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ファイルを生成してみるとこんな感じ。
f:id:cochineal19:20210829154816p:plain

PDFのリンクも開け、画像も表示され、最低限の機能は揃えられたかな?

一元管理ができそうで、何より何より。


今回作ったHTMLテーブルを別のHTMLファイルに埋め込む方法は以下の記事で。
【HTML】別HTMLを埋め込む(iframe) - こちにぃるの日記

【VBA】パラメタを変えてスクリプトを大量生成

久々にVBA

R、PythonSAS ... なんでも良いが、
共有ロジックを外部ファンクションや外部マクロにしておいて、 パラメタを変えてぐるぐる回すことがある。
例えば、疾患名を変えて患者数を推移図にするなど。

そんなとき、スクリプトをいちいち作っていると面倒だし、ミスにもつながる。

だからVBAで一括生成してしまおうというお話。


こんな感じでエクセルに情報を埋めて、
f:id:cochineal19:20210828145754p:plain:w400

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スクリプト)。
f:id:cochineal19:20210828145957p:plain:w350

スクリプトの中身はこんな感じ。外部ファンクションへの指定パラメタをスクリプト毎に変えて生成している。
f:id:cochineal19:20210828150239p:plain:w550
f:id:cochineal19:20210828150319p:plain:w550

これで少し楽ができる。めでたし、めでたし。

【Python】Kerasでロジスティック回帰

Kerasの勉強用にロジスティック回帰と同じ構造を作って実行してみる。

【目次】


基本の整理


ニューラルネットワークは入力層ー中間層ー出力層からなり、図示するとこんな感じ。

f:id:cochineal19:20210725104049p:plain:w350

このうち、ロジスティック回帰は中間層がない単純な形。

f:id:cochineal19:20210725104305p:plain:w220

活性化関数にはシグモイド関数が使われる。

\quad y= \dfrac{1}{1+exp \left( - w^{T} X \right)}


また、損失関数には交差エントロピー誤差が使われる。

\quad E=-\sum \{ t_{i} \cdot log\ y_{i} + \left(1-t_{i} \right) \cdot log \left( 1-y_{i} \right) \}


使用データ


Kaggleチュートリアルタイタニックデータを使う。

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))

f:id:cochineal19:20210725114727p:plain:w400

アルゴリズムが違うだろうが、だいたい同じ?

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クラスロジスティック回帰 - 人工知能に関する断創録

本ブログで参考になりそうなもの
【統計】ロジスティック回帰分析 - こちにぃるの日記

本ブログは個人メモです。 本ブログの内容によって生じた損害等の一切の責任を負いかねますのでご了承ください。