H2OとRでグリッドサーチする

H2Oはグリッドサーチを使うことでハイパーパラメータの探索をある程度、効率的にできます。Rを使った場合のグリッドサーチの書き方を説明します。

grid <- h2o.grid(
  algorithm = "<アルゴリズム名>",
  x = <説明変数リスト>,
  y = <目的変数>,
  training_frame = <訓練集合>,
  validation_frame = <評価集合>,

  hyper_params = list(
     <ハイパーパラメータのリスト>
  )
)

GBMなどはそうハイパーパラメータの数が多くないので学習をただ繰り返してもいいのですが、Deep Learningなどはハイパーパラメータも多いのでグリッドサーチがないとチューニングが現実的ではなくなります。

Cox回帰の結果を使ってMCMCしてみようとした

まず、Cox回帰は以下のように定義されます。すなわち、ある瞬間での瞬間死亡率を積分して、累積死亡率を出します。このとき、h_0(t)がベースラインハザード関数と呼ばれ、時間依存の死亡率を表します。そして、それに共変量の影響を加えますが影響は対数線形モデルで書けること、さらに時間によらず一定なこと(比例ハザード性)が仮定されています。

h(t)=\lim_{\Delta t \to 0} \frac{S(t)-S(t - \Delta t)}{\Delta t} \cdot \frac{1}{S(t)}
H(t)=\int_0^{t} h(t)du=-\log S(t)
h(t)=h_0(t) \cdot \mathrm{e}^{(\alpha + \beta_{1} X_{1} + \beta_{2} X_{2} + \ldots + \beta_{p} X_{p})} = h_0(t) \cdot RISK

たとえば、仮にCox回帰でゲームをやめてしまうモデルを作ったとすれば、お客さんがゲームをやめてしまう可能性はたとえば性別だとかそういったものを説明変数とした場合、その影響は時間によらず一定なことが求められます。というか、そうであることが仮定されているモデルだからですが。

裏を返せば、説明変数の影響は時間によらず一定なので、共変量の影響部分だけを切り離してスコア化できるモデルともとらえられると思います。つまり、先の式で言うと\mathrm{e}^{(\alpha + \beta_{1} X_{1} + \beta_{2} X_{2} + \ldots + \beta_{p} X_{p})}の部分を切り離して顧客のリスクとしてとらえられます。

もともと、Cox回帰は生存時間分析の一角で医学で用いられてきたため、用語として瞬間死亡確率、累積死亡確率というなかなかヘビーな単語が頻出します。今日では、離反予測などでも使われているのでイベントヒストリー分析なとども言います。

累積死亡率ですから、その区間の前までの累積死亡率からその前までの累積死亡率を引けばその区間での死亡率が求められるはずです。そうすると、その死亡率でさいころを振ってみてその区間で死亡するかどうかが求められます。その結果を集計すればその結果はほぼほぼ集計レベルでは元と近似するはずです。

ベースラインハザード関数の部分はノンパラメトリックなの推定ですから推定部分は定数としておくなり、データベースに格納するなりが考えられます。パラメトリックな共変量からの計算部分は関数としておけます。たとえば、Rのサンプルデータのkidneyのモデルはこのようにおけるでしょうか。

let KidneyModel sex diseaseGN diseaseAN diseasePKD = 
    let x00 = sex - 1.73684210526316  in
        let x01 = diseaseGN - 0.236842105263158 in
            let x02 = diseaseAN - 0.315789473684211  in
                let x03 = diseasePKD - 0.105263157894737  in
                    System.Math.Exp(-1.47739255480755 * x00 + 0.139220775300325  * x01 + 0.413157096502011  * x02 + (-1.36707122522468 ) * x03)

このコードはベースラインハザード関数の部分、先ほどの式ではh_0(t)の部分を含んでいないので式としては不完全です。最初で共変量から引いているのはRのcoxphは内部で共変量から平均を引き算しているため同様な処理を入れています。当然ながら、書き下した関数とRの結果がほぼほぼ誤差込みで一致するのは確認済みです。

Rのモデルの情報はこんな感じです。

Call:
coxph(formula = Surv(time, status) ~ sex + disease, data = kidney, 
    model = T)


             coef exp(coef) se(coef)     z       p
sex        -1.477     0.228    0.357 -4.14 3.5e-05
diseaseGN   0.139     1.149    0.364  0.38    0.70
diseaseAN   0.413     1.512    0.336  1.23    0.22
diseasePKD -1.367     0.255    0.589 -2.32    0.02

Likelihood ratio test=17.6  on 4 df, p=0.0015
n= 76, number of events= 58 

このとき、KidneyBaseline t なる時刻tを引数としてもち、時刻tのベースラインハザード曲線を返す関数があったならば、時刻tにおける累積死亡確率は

h(t)=h_0(t) \cdot \mathrm{e}^{(\alpha + \beta_{1} X_{1} + \beta_{2} X_{2} + \ldots + \beta_{p} X_{p})}

に基づき、(KidneyModel sex diseaseGN diseaseAN diseasePKD) * (KidneyBaseline t)になります。従って。累積死亡確率を求める関数 KidneyCumsumは以下のように定義されます。

let KidneyCumsum t sex diseaseGN diseaseAN diseasePKD = (KidneyModel sex diseaseGN diseaseAN diseasePKD) * (KidneyBaseline t)

その上で、区間t1からt2における死亡確率は(KidneyCumsum t2) – (KidneyCumsum t1)になります。従って、KidneyDeathProbabilityなる関数を定義するとすれば

let KidneyDeathProbability t1 t2 sex diseaseGN diseaseAN diseasePKD = (KidneyCumsum t2) - (KidneyCumsum t1)

となります。

これをサンプルデータ一つ一つについて、順次計算し、同じ区間について集計すれば平均的なその区間における死者を推定することができます。Cox回帰は比例ハザード性を仮定しているので共変量の効き具合は時間にかかわらずという強力な仮定があるのでこの、仮定が崩れていると使えませんので注意がいります。

これらを踏まえてMCMCを実装するとこのようになります。

open Kidney;

type Row = {sex: float; diseaseGN: float; diseaseAN: float; diseasePKD: float; mutable dead: bool}

[]
let main argv = 
    let sr = new System.IO.StreamReader("kidneyp.csv", System.Text.Encoding.ASCII)
    let mutable buffer = sr.ReadLine()
    let getGN (buffer:string) =
        match buffer.Replace("\"", "") with
        | "GN" -> 1.0
        | _ -> 0.0
    let getAN (buffer:string) =
        match buffer.Replace("\"", "") with
        | "AN" -> 1.0
        | _ -> 0.0
    let getPKD (buffer:string) =
        match buffer.Replace("\"", "") with
        | "PKD" -> 1.0
        | _ -> 0.0

    let parse (buffer:string) =
        let items = buffer.Split [|','|]
        {sex = System.Double.Parse items.[5];diseaseGN = getGN items.[6];diseaseAN = getAN items.[6];diseasePKD = getPKD items.[6];dead = false}

    let patients = [|
        while not sr.EndOfStream do
            buffer <- sr.ReadLine()
            yield parse buffer
    |]

    let rg = new System.Random()

    let results = [|
        for i in 1..(KidneyModel.KidneyBaselineSize - 1) do
            let ps = [|
                for patient in patients do
                    if patient.dead = false then
                        let p = KidneyModel.KidneyDeathProbability (i - 1) i patient.sex patient.diseaseGN patient.diseaseAN patient.diseasePKD
                        let r = rg.NextDouble()
                        let c = match (r <= p) with
                        | true -> 1.0
                        | _ -> 0.0
                        if c > 0.0 then
                            patient.dead <- true
                        yield c
                    else
                        yield 0.0
                done
            |]
            yield Array.sum(ps)
        done
    |]

    let sw = new System.IO.StreamWriter("dead.csv", false, System.Text.Encoding.ASCII);
    for result in results do
        sw.WriteLine(result)
    done
    sw.Close()

    printfn "%A" argv
    0 // 整数の終了コードを返します

シミュレーションの結果は以下の通りです。

集計結果

集計結果

ソードアート・オンラインをword2vecする

ソードアート・オンライン (Web版)をword2vecしてみました。これをmecabを用いて分かち書きします。分かち書きしたものをgensimに与えて、word2vecでモデルを作成してみました。

モデリングのコードは以下のとおりです。基本的にあんちべさんの『自然言語処理の最新手法”word2vec”で艦これ加賀さんから乳を引いてみる』と同じです。

import codecs
from gensim.models import word2vec
import logging
logging.basicConfig(format='%(asctime)s : %(levelname)s : %(message)s', level=logging.INFO)
sentences = word2vec.Text8Corpus('SAO1.wakati')
model = word2vec.Word2Vec(sentences, size=200)

原作でキリトの対するヒースクリフについてアウトプットを見てみました。

ランク 類似度
1 自分 0.797701954842
2 茅場 0.772888898849
3 姿 0.766229510307
4 意思 0.732955038548
5 0.720785975456

茅場はヒースクリフの正体でありセットで出てくるので類似度が高いのは妥当です。自分は恐らく、1巻がキリト目線だからでしょう。キリトの前に立ちはだかるのがヒースクリフであり、キリト自身意識しているため高い類似度を示したものと思います。紅は恐らくヒースクリフがKoBの団長であるためでしょう。

この結果を見る限り、word2vecのアウトプットはかなり妥当性の高いものではないかと思います。

AzureMLでfitdistrplusを使った確率分布の推定

ここでは確率分布の種類は既知のものとして、データから確率分布のパラメータの推定を行う過程を紹介する。Azure MLでは当該機能は機能セットとしては提供していないため、Rで実装されているfitdistrplusを用いて行う。
サンプルデータとしてパラメータが既知の分布を基づいて乱数を生成し、その乱数をデータとしてパラメータの推定を行う。この場合のML StudioでのExperimentを下記図に示す。

Experiment for Estimation Parameters

Experiment for Estimation Parameters

Figure 1は2つのRスクリプトで構成されている。このスクリプトを下記Source 1およびSource 2に示す。Source 1はサンプルデータの生成を行うスクリプトで、Source 2がデータに基づき分布の推定を行うスクリプトである。

maml.mapOutputPortおよび、maml.mapInputPortはAzure MLのRにおける固有の関数でデータの受け渡しを行う関数である。アウトプットされたグラフィックスは既定で特定の割り当てがなされるので特に実装を行わなくてもPNGフォーマットで取得することができる。出力結果を以下のEstimate Resultに示す。

# Create Data Set
data.set
Source 1: Data Generate
# Load required packages
library(fitdistrplus)

# load data frame
data.set
Source 2: Parameter Estimation
Estimate Result

Estimate Result


データサイエンティスト養成読本 R活用編

いささか、旧聞にはなるが上記、書籍を書いておきます。基本的には前に刊行されたムックの続きでRなどにフォーカスしたムックです。R活用編となってはいますが、商業書籍では珍しいJuliaの入門記事やRを.NET環境から利用するR.NETの記事など割と面白いものになっています。特にR.NETを含めた.NETのライブラリ群の配布は現在はNuGetであることは.NETの開発者の諸兄には自明なのですが、意外と気づかない向きもありますので書いておきます。

さて、私の目的はJuliaの入門記事であったわけですが、文法の解説に関しては割りとしっかり書かれているように思われます。とはいえ、ムックの一セクションなので食い足りないのも確かです。書籍ベースだと後は洋書の”Seven More Languages in Seven Weeks: Languages That Are Shaping the Future”と先に書評を書いた、同人誌の実例で学ぶJulia言語入門 v0.3.3/v0.4.0-dev 対応版あたり。

PowerQueryでRemember the milkからタスク一覧を抽出する (拡張)

先ほどのPowerQueryのクエリを拡張して他のカラムも取り込んで見ました。タグや期日なども取り込めるようにしてのでサンプルとしてより実用的になったと思います。

let
    Source = Web.Contents("< source uri >"),
    Document = Xml.Document(Source),
    Entries = Table.SelectRows(Document{0}[Value], each ([Name] = "entry")),
    Rows = Table.ToRows(Entries),
    ProjectedEntries = List.Transform(Rows, (x) => x{2}),
    P = List.Transform(ProjectedEntries, (x) => Table.SelectRows(x, each (([Name] = "id") or ([Name] = "updated") or ([Name] = "title") or ([Name] = "content")))),
    R = Table.FromList(List.Transform(P, (x) => Table.Range(Table.PromoteHeaders(Table.Transpose(x)),1,1)), Splitter.SplitByNothing(), null, null, ExtraValues.Error),
    S = Table.ExpandTableColumn(R, "Column1", {"updated", "id", "title", "content"}, {"updated", "id", "title", "content"}),
    SX = Table.AddColumn(S, "additional", each Table.PromoteHeaders(Table.Transpose(Table.FromRecords(List.Transform([content]{0}[Value][Value],(x) => [Name = Text.Trim(Text.Replace(x[Value]{0},":","")),Value = x[Value]{1}]))))),
    SY = Table.ExpandTableColumn(SX, "additional", {"Due", "Priority", "Time estimate", "Tags", "Location", "Postponed"}, {"additional.Due", "additional.Priority", "additional.Time estimate", "additional.Tags", "additional.Location", "additional.Postponed"}),
    S1 = Table.RenameColumns(SY,{{"updated", "更新日時"},{"title", "タスク"},{"id", "ID"},{"additional.Due", "期日"},{"additional.Priority", "優先度"},{"additional.Tags", "タグ"}}),
    S2 = Table.RemoveColumns(Table.ReorderColumns(S1,{"更新日時", "タスク", "期日", "優先度", "タグ", "ID"}),{"content","additional.Time estimate","additional.Location","additional.Postponed"})
    
in
    S2

PowerQueryでRemember the milkのRSSからタスク一覧を抽出する

PowerQueryでRemember the milkのフィードを加工してタスク一覧を抽出するクエリを書いてみました。データ加工のサンプルとして最適なので公開します。フィードのURIのみ変更してあります。

let
    Source = Web.Contents(< feed address >),
    Document = Xml.Document(Source),
    Entries = Table.SelectRows(Document{0}[Value], each ([Name] = "entry")),
    Rows = Table.ToRows(Entries),
    ProjectedEntries = List.Transform(Rows, (x) => x{2}),
    P = List.Transform(ProjectedEntries, (x) => Table.SelectRows(x, each (([Name] = "id") or ([Name] = "updated") or ([Name] = "title") or ([Name] = "content")))),
    R = Table.FromList(List.Transform(P, (x) => Table.Range(Table.PromoteHeaders(Table.Transpose(x)),1,1)), Splitter.SplitByNothing(), null, null, ExtraValues.Error),
    S = Table.ExpandTableColumn(R, "Column1", {"updated", "id", "title", "content"}, {"updated", "id", "title", "content"}),
    S0 = Table.AddColumn(S, "filtered content", each [content]),
    S1 = Table.RenameColumns(S,{{"updated", "更新日時"},{"title", "タスク"},{"id", "ID"}}),
    S2 = Table.RemoveColumns(Table.ReorderColumns(S1,{"更新日時", "タスク", "content", "ID"}),{"content"})
    
in
    S2