Gradle で Jetty9 を使用

Java

今のところ Gradle 標準の jetty プラグインでは Jetty6 しか使えないようなので、Servlet 3.0 を使いたいケースでは不便です。

この場合、build.gradle で Jetty9 の起動処理を自前で実装する手も考えられますが、Gradle Jetty Plugin for Eclipse Jetty （gradle-jetty-eclipse-plugin） を使った方が簡単だと思います。 （他にも gradle-jetty9-plugin というのもあります ）

gradle-jetty-eclipse-plugin を使うと gradle jettyEclipseRun で Jetty 9.0.6 を起動できます。

gradle-jetty-eclipse-plugin の使い方

使用するには、build.gradle へ下記の設定を追加するだけです。

  • gradle-jetty-eclipse-plugin を使用するための buildscript を定義
  • jettyEclipse を apply plugin する
gradle-jetty-eclipse-plugin を使用するための build.gradle 設定
apply plugin: 'jettyEclipse'

buildscript {
    repositories {
        jcenter()
        maven {
            url 'http://dl.bintray.com/khoulaiz/gradle-plugins'
        }
    }
    dependencies {
        classpath (group: 'com.sahlbach.gradle', name: 'gradle-jetty-eclipse-plugin', version: '1.9.+')
    }
}

なお、jcenter() の部分は Jetty9 のモジュールを取得できるリポジトリであればよいようなので mavenCentral() でも問題ないようです。

Jetty9 の実行

それでは @WebServlet アノテーションを使った単純な Servlet を実行してみます。

サンプルソースは http://github.com/fits/try_samples/tree/master/blog/20140203/ に配置しています。

まずは build.gradle です。

gradle-jetty-eclipse-plugin の利用設定を追加し、providedCompile で Servlet 3.0 の API を指定しました。

build.gradle
apply plugin: 'jettyEclipse'

buildscript {
    repositories {
        jcenter()
        maven {
            url 'http://dl.bintray.com/khoulaiz/gradle-plugins'
        }
    }
    dependencies {
        classpath (group: 'com.sahlbach.gradle', name: 'gradle-jetty-eclipse-plugin', version: '1.9.+')
    }
}

repositories {
    mavenCentral()
}

dependencies {
    providedCompile 'javax.servlet:javax.servlet-api:3.0.1'
}

Servlet は "sample" という文字を出力するだけの単純な処理を実装しました。

SampleServlet.java
package fits.sample;

import java.io.IOException;
import javax.servlet.ServletException;
import javax.servlet.annotation.WebServlet;
import javax.servlet.http.HttpServlet;
import javax.servlet.http.HttpServletRequest;
import javax.servlet.http.HttpServletResponse;

@WebServlet(urlPatterns = {"/sample"})
public class SampleServlet extends HttpServlet {
    @Override
    protected void doGet(HttpServletRequest req, HttpServletResponse res) throws ServletException, IOException {
        res.getWriter().print("sample");
    }
}

それでは gradle jettyEclipseRun で実行してみます。

ちなみに、標準 jetty プラグインのような jettyRun と jettyRunWar のような区別はありません。

実行
> gradle jettyEclipseRun

:compileJava
:processResources UP-TO-DATE
:classes
:war
:jettyEclipseRun
Empty contextPath
!RequestLog

Hit <ENTER> to reload the webapp.
Hit r + <ENTER> to rebuild and reload the webapp.
Hit R + <ENTER> to rebuild the webapp without reload

> Building 80% > :jettyEclipseRun > Running at http://localhost:8080/

これで Jetty9 が起動し、http://localhost:8080/sample へ接続すると "sample" という文字が返ってくるはずです。

Commons OGNL でマッピング・フィルタリング・畳み込み

Java Groovy

以前、Java でマッピング・フィルタリング・畳み込みを試しましたが、 今回は Commons OGNL を使って OGNL 式によるマッピング・フィルタリング・畳み込みを Groovy で試してみました。

サンプルソースは http://github.com/fits/try_samples/tree/master/blog/20140118/

はじめに

Commons OGNL ではマッピングとフィルタリング用の式は用意されてますが、今のところ畳み込みは用意されていないようです。 （一応、畳み込みはラムダ式 :[ e ] で実装できました）

処理 機能名 OGNL式
マッピング Projection e1.{ e2 }
フィルタリング Selection e1.{? e2 }
畳み込み - -

Selection には、マッチした最初の要素を返す e1.{^ e2 } や最後の要素を返す e1.{$ e2 } のようなバリエーションも用意されています。

サンプルで使用するモデルクラス

今回のサンプルでは下記のようなモデルクラスを使いました。

Order.groovy
import groovy.transform.*

@CompileStatic
class Order {
    List<OrderLine> lines = []
}

@CompileStatic
@Immutable
class OrderLine {
    String code
    BigDecimal price = 0
}

groovyc でコンパイルしておきます。

コンパイル
> groovyc Order.groovy

マッピング・フィルタリング

OGNL 式を使ったマッピング・フィルタリング処理を Groovy スクリプトで試してみます。

まず、@Grab を使って普通に Ognl.getValue(<OGNL式>, <オブジェクト>) を実行しようとすると下記のようなエラーが発生します。

エラー内容
java.lang.IllegalArgumentException: Javassist library is missing in classpath! Please add missed dependency!
        at org.apache.commons.ognl.OgnlRuntime.getCompiler(OgnlRuntime.java:210)        ・・・
Caused by: java.lang.ClassNotFoundException: Unable to resolve class: javassist.ClassPool
        at org.apache.commons.ognl.OgnlRuntime.classForName(OgnlRuntime.java:665)
        ・・・

これは、デフォルトで使用される DefaultClassResolver が SystemClassLoader から javassist.ClassPool をロードしようとする事に起因するようで、Servlet 等の Java EE Web アプリケーションで実行する際にも同様のエラーが発生します。 （war ファイル内の JAR からロードするような場合）

エラー内容から CLASSPATH に javassist.jar が含まれていないように思ってしまいますが、そうではありません。

DefaultClassResolver に原因があるので、自前の ClassResolver を使って作成した Context を使うようにすればエラーを回避できます。

map_filter_sample.groovy
@GrabResolver('http://repository.apache.org/snapshots/')
@Grab('org.apache.commons:commons-ognl:4.0-SNAPSHOT')
import org.apache.commons.ognl.Ognl
import org.apache.commons.ognl.ClassResolver

def data = new Order()
data.lines << new OrderLine('1', 100)
data.lines << new OrderLine('2', 200)
data.lines << new OrderLine('3', 300)

try {
    // (1) DefaultClassResolver が原因でエラーが発生
    println Ognl.getValue('lines.{? #this.code == "2" }', data)
} catch (e) {
    println '(1) ' + e
    //e.printStackTrace()
}
println '-----------------------'

// エラー回避のために自前の ClassResolver を使って Context 作成
def ctx = Ognl.createDefaultContext(null, { String className, Map<String, Object> context ->
    Class.forName(className)
} as ClassResolver)

// (2) マッピング
println '(2) ' + Ognl.getValue('lines.{ #this.price > 100 }', ctx, data)

// (3) フィルタリング
println '(3) ' + Ognl.getValue('lines.{? #this.price > 100 }', ctx, data)

// (4) マッチした最初の要素
println '(4) ' + Ognl.getValue('lines.{^ #this.price > 100 }', ctx, data)

// (5) マッチした最後の要素
println '(5) ' + Ognl.getValue('lines.{$ #this.price > 100 }', ctx, data)

// (6) OGNL式をコンパイルして使用
def exprNode = Ognl.compileExpression(ctx, null, 'lines.{? #this.code not in {"2", "4"} }')
println '(6) ' + Ognl.getValue(exprNode, data)

実行結果は下記の通りです。

実行結果
> groovy map_filter_sample.groovy
(1) java.lang.IllegalArgumentException: Javassist library is missing in classpath! Please add missed dependency!
-----------------------
(2) [false, true, true]
(3) [OrderLine(2, 200), OrderLine(3, 300)]
(4) [OrderLine(2, 200)]
(5) [OrderLine(3, 300)]
(6) [OrderLine(1, 100), OrderLine(3, 300)]

畳み込み

OGNL のラムダ記法 :[ 処理 ] を使って畳み込みを試しに実装してみました。

OGNL のラムダは引数を一つしか取れないようなので（ラムダ内の #this にバインドされます）、リストを使ってラムダ #fold へ引数 {<処理>, <値>, <要素リスト>} を渡すようにしています。

更に、OrderLine の price を合計するための処理もラムダで定義して （:[ #this[0] + #this[1].price ]）、#fold へ渡しています。

fold_sample.groovy
・・・
// 畳み込みの OGNL 式
def foldOgnl = '''
#fold = :[
  #this[2].size() > 0 ? #fold({ #this[0], #this[0]({ #this[1], #this[2][0] }), #this[2].subList(1, #this[2].size()) }) : #this[1]
],
#fold({ :[ #this[0] + #this[1].price ], 0, lines })
'''
// 0 + 100 + 200 + 300
println Ognl.getValue(foldOgnl, ctx, data)
実行結果
> groovy fold_sample.groovy
600

なお、合計を計算するだけなら以下のようにした方がシンプルです。

println Ognl.getValue('#v = 0, lines.{ #v = #v + #this.price }, #v', ctx, data)

並列実行時の注意点

最後に、並列処理で実行する際の注意点です。

問題確認のために、下記 3種類の処理を 50回並列でそれぞれ行ってみて OGNL の処理結果が正しかったかどうかを true・false で返すようにしてみます。

  • (1) 文字列の OGNL 式で Context を再利用して getValue を実行
  • (2) 文字列の OGNL 式で Context を都度作成して getValue を実行
  • (3) compileExpression した結果で getValue を実行
parallel_sample.groovy
@GrabResolver('http://repository.apache.org/snapshots/')
@Grab('org.apache.commons:commons-ognl:4.0-SNAPSHOT')
import org.apache.commons.ognl.Ognl
import org.apache.commons.ognl.ClassResolver

import groovyx.gpars.*

def createData = { i ->
    def data = new Order()
    data.lines << new OrderLine('1', i)
    data.lines << new OrderLine('2', i + 1)
    data.lines << new OrderLine('3', i + 2)
    data
}

def sum = { data ->
    data.lines.inject(0){ acc, val ->
        acc + val.price
    }
}

def printResult = {
    it.groupBy().each { k, v -> println "${k} ${v.size()}" }
}

def createContext = {
    Ognl.createDefaultContext(null, { String className, Map<String, Object> context ->
        Class.forName(className)
    } as ClassResolver)
}

def ctx = createContext()

// OGNL式をコンパイル
def exprNode = Ognl.compileExpression(ctx, null, '#v = 0, lines.{ #v = #v + #this.price }, #v')

def count = 50

GParsPool.withPool(20) {
    // (1) 文字列の OGNL 式で並列処理（Context再利用）
    def res1 = (0..<count).collectParallel {
        def d = createData(it)
        try {
            // 稀に NoSuchPropertyException: Order.price が発生するため try-catch
            sum(d) == Ognl.getValue('#v = 0, lines.{ #v = #v + #this.price }, #v', ctx, d)
        } catch (e) {
            println e
            false
        }
    }

    // (2) 文字列の OGNL 式で並列処理（Context毎回作成）
    def res2 = (0..<count).collectParallel {
        def d = createData(it)
        sum(d) == Ognl.getValue('#v = 0, lines.{ #v = #v + #this.price }, #v', createContext(), d)
    }

    // (3) compileExpression 結果で並列処理
    def res3 = (0..<count).collectParallel {
        def d = createData(it)
        sum(d) == Ognl.getValue(exprNode, d)
    }

    println '----- (1) 文字列の OGNL 式で並列処理（Context再利用）------'
    printResult res1

    println '----- (2) 文字列の OGNL 式で並列処理（Context毎回作成） ---'
    printResult res2

    println '----- (3) compileExpression 結果で並列処理 ----------------'
    printResult res3
}

処理結果は以下のようになり、(1) のケースで問題が発生しています。

実行結果
> groovy parallel_sample.groovy
org.apache.commons.ognl.NoSuchPropertyException: Order.price
----- (1) 文字列の OGNL 式で並列処理（Context再利用）------
true 36
false 14
----- (2) 文字列の OGNL 式で並列処理（Context毎回作成） ---
true 50
----- (3) compileExpression 結果で並列処理 ----------------
true 50

これは (1) のように Context を再利用すると #v や #this を並列処理間で共有してしまう事が原因と考えられます。

そのため、並列処理で実行する際は Context をその都度作成するか、compileExpression した結果を使う事になると思います。 （他の方法もあるかもしれません）

MathJax + AngularJS で数式の編集と表示

JavaScript

LaTeX の数式表示に JavaScript・HTML・CSS といった Web 技術で数式を表示できる MathJax が便利です。

簡易な数式編集ツールなら簡単に作れそうだったので、AngularJS を試すついでに以下のようなサンプルを作ってみました。

  • TextArea に入力した数式を動的にプレビュー表示

サンプルソースは http://github.com/fits/try_samples/tree/master/blog/20140111/

数式編集のサンプル

MathJax で特定の要素の数式表示を更新するには下記のように処理します。

  • (1) MathJax が適用された要素を取得
  • (2) MathJax.Hub.Queue() を使って (1) の要素へ数式を表示

なお、数式を表示する要素には MathJax の適用対象となるように ${}$ を設定しておきます。

下記サンプルでは、TextArea の内容に変化が生じるとプレビュー表示を更新するよう AngularJS の ng-change を使いました。

index.html
<!doctype html>
<html ng-app>
<head>
    <title>MathJax + AngularJS sample</title>
    <script type="text/x-mathjax-config">
        MathJax.Hub.Config({
            tex2jax: {inlineMath: [["$","$"],["\\(","\\)"]]}
        });
    </script>
    <script src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML"></script>
    <script src="https://ajax.googleapis.com/ajax/libs/angularjs/1.2.8/angular.min.js"></script>
    <script>
        function mathjax($scope) {
            $scope.formula = '';

            $scope.changeFormula = function(formula) {
                // (1) 要素の取得
                var el = MathJax.Hub.getAllJax('result')[0];
                // (2) 数式の表示
                MathJax.Hub.Queue(['Text', el, formula])();
            };
        }
    </script>
</head>
<body>
    <div ng-controller="mathjax">
        <textarea ng-model="formula" ng-change="changeFormula(formula)"></textarea>
    </div>
    <br />
    <!-- プレビュー表示 -->
    <div id="result">${}$</div>
</body>
</html>
画面例

f:id:fits:20140111195353p:plain

R でロジスティック回帰 - glm, MCMCpack

R GLM Bayes

前回 に続き、今回も書籍 「 データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学) 」のサンプルを使って GLM とベイズ統計を試してみます。

題材は、6章「GLM の応用範囲をひろげる -ロジスティック回帰など-」のサンプルデータを使ったロジスティック回帰です。

  • (1) GLM によるロジスティック回帰 （glm 関数）
  • (2) MCMCpack を使ったベイズ統計によるロジスティック回帰（MCMCmetrop1R 関数）

ここで、MCMCpack のロジスティック回帰に MCMClogit 関数を使えると思ったのですが、使い方がよく分からなかったので MCMCmetrop1R 関数のみ利用しました。

サンプルソースは http://github.com/fits/try_samples/tree/master/blog/20131229/

はじめに

MCMCpack パッケージを install.packages("MCMCpack") で R の実行環境へインストールしておきます。

ロジスティック回帰を試すためのデータは書籍のサポート web サイトから取得できます。

データ data4a.csv
N,y,x,f
8,1,9.76,C
8,6,10.48,C
8,5,10.83,C
・・・

データ内容は下記のようになっており、個体 i それぞれにおいて 「 { N_i } 個の観察種子のうち生きていて発芽能力があるものは { y_i } 個、死んだ種子は { N_i - y_i } 個」 となっています。

項目 内容
N 観察種子数
y 生存種子数
x 植物の体サイズ
f 施肥処理 （C: 肥料なし, T: 肥料あり）

体サイズ x と肥料による施肥処理 f が種子の生存する確率（ある個体 i から得られた 1個の種子が生存している確率）にどのように影響しているかをロジスティック回帰で解析します。

(1) GLM によるロジスティック回帰 （glm 関数）

まず、glm 関数を使ったロジスティック回帰を行って予測値の曲線を描画します。

特徴は以下のようになります。

  • family に binomial （二項分布）を指定 （デフォルトのリンク関数が logit となる）
  • family が binomial で応答変数 （cbind(y, N - y) の箇所） が 2値（0, 1）以外なら cbind する必要あり （2値なら y だけでよい）

下記では、生存する確率 { q_i } の関数が { q_i = \frac{1}{1 + \exp(-z_i)} } 、
線形予測子 { z_i = \beta_1 + \beta_2 x_i + \beta_3 f_i } で解析する事になります。

logisticGlm.R
d <- read.csv('data4a.csv')

d.res <- glm(cbind(y, N - y) ~ x + f, data = d, family = binomial)

summary(d.res)

png("logisticGlm.png")
# C: 肥料なしを 赤、T: 肥料ありを 青 で描画
plot(d$x, d$y, col = c("red", "blue")[d$f])

xx <- seq(min(d$x), max(d$x), length = 1000)
ft <- factor("T", levels = levels(d$f))
fc <- factor("C", levels = levels(d$f))
# 下記でも可
#ft <- factor("T", levels = c("C", "T"))
#fc <- factor("C", levels = c("C", "T"))

# T: 肥料ありの予測値を算出
qq.t <- predict(d.res, data.frame(x = xx, f = ft), type="response")
# C: 肥料なしの予測値を算出
qq.c <- predict(d.res, data.frame(x = xx, f = fc), type="response")

# T: 肥料ありの予測値の曲線を 緑 で描画
lines(xx, max(d$N) * qq.t, col = "green")
# C: 肥料なしの予測値の曲線を 黄 で描画
lines(xx, max(d$N) * qq.c, col = "yellow")

dev.off()

predict 関数を使えば glm 結果の回帰モデルから新しいデータによる予測値を求めることができますので、"T:肥料あり" と "C:肥料なし" のそれぞれの予測値を求めて曲線を描画しています。

predict の結果は、生存する確率 { q_i } の予測値ですので、縦軸が生存種子数 { y_i } のグラフへ描画するには観察種子数 N （ここでは 8）を乗算する事になります。

ここで、"T:肥料あり" と "C:肥料なし" は factor 関数を使ってファクタ（因子）として作成しています。 d$f と同じ水準を持たせる必要がありますので levels 指定しています。

また、predict には glm で使った説明変数と同じ名前（上記の x と f）を使ったデータを渡す点に注意が必要です。

実行

実行すると下記のような結果になりました。

> R CMD BATCH logisticGlm.R
実行結果 logisticGlm.Rout
・・・
> d.res <- glm(cbind(y, N - y) ~ x + f, data = d, family = binomial)
> 
> summary(d.res)

Call:
glm(formula = cbind(y, N - y) ~ x + f, family = binomial, data = d)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.2584  -0.7948   0.1753   0.8757   3.1589  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -19.5361     1.4138  -13.82   <2e-16 ***
x             1.9524     0.1389   14.06   <2e-16 ***
fT            2.0215     0.2313    8.74   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 499.23  on 99  degrees of freedom
Residual deviance: 123.03  on 97  degrees of freedom
AIC: 272.21

Number of Fisher Scoring iterations: 5
・・・

この結果より、線形予測子 { z_i = \beta_1 + \beta_2 x_i + \beta_3 f_i } のパラメータ推定値 { { \hat\beta_1, \hat\beta_2, \hat\beta_3 } } が -19.5361, 1.9524, 2.0215 となります。

f:id:fits:20131229184250p:plain

(2) MCMCpack を使ったベイズ統計によるロジスティック回帰（MCMCmetrop1R 関数）

次は、MCMCpack の MCMCmetrop1R 関数を使ったロジスティック回帰です。

今回のモデルの対数尤度は { \log L ( \beta_j ) = \sum_i \log( \small N_i \large C \small y_i \normalsize ) + y_i \log(q_i) + (N_i - y_i) \log(1 - q_i) } で、
{ q_i = \frac{1}{1 + \exp(-z_i)} } 、
{ z_i = \beta_1 + \beta_2 x_i + \beta_3 f_i } なので、
これらを関数化 （下記の func） して MCMCmetrop1R 関数へ渡しています。

ちなみに、{ \small N_i \large C \small y_i \normalsize = \frac{N_i !}{y_i ! (N_i - y_i) ! } } です。（combination）

前回は、項目毎のデータ（d$x や d$y）を func 関数へ渡すようにしていましたが、今回はデータ d を丸ごと渡すようにしました。 MCMCmetrop1R(・・・, data = d, ・・・)

ファクタ（因子）のデータ （d$f） を直接計算に使えないようなので as.numeric で数値化して使っています。

なお、{ \log( \small N_i \large C \small y_i \normalsize ) + y_i \log(q_i) + (N_i - y_i) \log(1 - q_i) } の部分は dbinom と log 関数を使うとシンプルになります。

パラメータ { \beta_1 } { \beta_2 } { \beta_3 } の初期値は theta.init で c(0, 0, 0) と指定しました。

logisticMcmcMetrop.R
library(MCMCpack)

d <- read.csv('data4a.csv')

func <- function(beta, data) {
  z <- beta[1] + beta[2] * data$x + beta[3] * as.numeric(data$f)
  q <- 1.0 / (1.0 + exp(-z))
  sum(log(choose(data$N, data$y)) + data$y * log(q) + (data$N - data$y) * log(1 - q))

  # 下記でも可
  # sum(log(dbinom(data$y, data$N, q)))
}

d.res <- MCMCmetrop1R(func, theta.init = c(0, 0, 0), data = d, logfun = TRUE)

summary(d.res)

png("logisticMcmcMetrop.png")

plot(d$x, d$y, col = c("red", "blue")[d$f])

xx <- seq(min(d$x), max(d$x), length = 1000)
ft <- factor("T", levels = levels(d$f))
fc <- factor("C", levels = levels(d$f))

zt <- mean(d.res[,1]) + mean(d.res[,2]) * xx + mean(d.res[,3]) * as.numeric(ft)
zc <- mean(d.res[,1]) + mean(d.res[,2]) * xx + mean(d.res[,3]) * as.numeric(fc)

lines(xx, max(d$N) * 1.0 / (1.0 + exp(-zt)), col="green")
lines(xx, max(d$N) * 1.0 / (1.0 + exp(-zc)), col="yellow")

dev.off()

MCMCmetrop1R 結果の各パラメータの mean の値を使って、
肥料の有無（T と C）の生存確率の予測値をそれぞれ算出して、
{ \frac{1}{1 + \exp(-(\hat\beta_1 + \hat\beta_2 x_i + \hat\beta_3 f_i))} }
観察種子数 N （ここでは 8）を乗算し曲線を描画しています。

実行

実行すると下記のような結果になります。

> R CMD BATCH logisticMcmcMetrop.R
実行結果 logisticMcmcMetrop.Rout
・・・
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
The Metropolis acceptance rate was 0.44590
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
> 
> summary(d.res)

Iterations = 501:20500
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 20000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

        Mean     SD  Naive SE Time-series SE
[1,] -21.694 1.5583 0.0110185       0.037186
[2,]   1.966 0.1386 0.0009798       0.003322
[3,]   2.027 0.2268 0.0016035       0.005383

2. Quantiles for each variable:

        2.5%     25%     50%     75%   97.5%
var1 -24.836 -22.733 -21.667 -20.596 -18.766
var2   1.703   1.870   1.963   2.059   2.248
var3   1.594   1.873   2.022   2.172   2.494
・・・

Mean の値が glm の結果と似たような値になりました。 グラフにすると違いが分からないくらいです。

f:id:fits:20131229184326p:plain

{ \beta_1 } の分布

f:id:fits:20131229184348p:plain

{ \beta_2 } の分布

f:id:fits:20131229184713p:plain

{ \beta_3 } の分布

f:id:fits:20131229184417p:plain

R でポアソン回帰 - glm, MCMCpack

R GLM Bayes

書籍 「 データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学) 」の 3章 「一般化線形モデル（GLM）」 と 9章 「GLMのベイズモデル化と事後分布の推定」 で説明されていたポアソン回帰を下記のような 3通りで試してみました。

書籍では、R から WinBUGS を呼び出して MCMC サンプリングを行っていましたが、今回は R 上でベイズ統計解析を実施する MCMCpack パッケージを試してみました。

サンプルソースは http://github.com/fits/try_samples/tree/master/blog/20131215/

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

はじめに

まず、MCMCpack パッケージを使用するためにパッケージを R の実行環境へインストールしておきます。

MCMCpack のインストール
> install.packages("MCMCpack")

次に、ポアソン回帰を試すデータを書籍のサポート web サイトから取得してきます。

データ data3a.csv
y,x,f
6,8.31,C
6,9.44,C
6,9.5,C
12,9.07,C
10,10.16,C
・・・

ここで、データの内容は下記のようになっており、体サイズ x と肥料による施肥処理 f が種子数 y にどのように影響しているかをポアソン回帰で解析します。

項目 内容
y 種子数
x 植物の体サイズ
f 施肥処理

(1) GLM によるポアソン回帰 （glm 関数）

それでは glm 関数を使ったポアソン回帰を試してみます。 ここでは stepAIC() を使った AIC によるモデル選択を試してみました。

poissonGlm.R
d <- read.csv('data3a.csv')

# 説明変数を全て投入したモデル （y ~ x + f と同じ）
d.all <- glm(y ~ ., data = d, family = poisson)

library(MASS)
# AIC によるモデル選択
d.res <- stepAIC(d.all)

summary(d.res)

png("poissonGlm.png")
plot(d$x, d$y, col = c("red", "blue")[d$f])
xx <- seq(min(d$x), max(d$x), length = 1000)
# y ~ x のモデルを使った平均種子数の予測値の曲線を描画
lines(xx, exp(d.res$coefficients["(Intercept)"] + d.res$coefficients["x"] * xx), col="green")
# 以下でも可
#lines(xx, predict(d.res, newdata = data.frame(x = xx), type = "response"), col = "green")
dev.off()

なお、今回は y ~ x のモデルが選択される事 （施肥処理は効果が無い） を予め分かっているので、y ~ x のモデルを使って平均種子数の予測値を描画 （緑の線） しています。

ここで、 y ~ x モデルのリンク関数 （対数リンク関数） は { \log \lambda_i = \beta_1 + \beta_2 x_i } で、平均種子数 { \lambda_i } は { \lambda_i = \exp (\beta_1 + \beta_2 x_i) } となります。

また、{ \beta_1 } の最尤推定値が d.res$coefficients["(Intercept)"]、{ \beta_2 } の最尤推定値が d.res$coefficients["x"] に該当します。

実行

実行すると下記のような結果になります。

> R CMD BATCH poissonGlm.R
実行結果 poissonGlm.Rout
・・・
> d.res <- stepAIC(d.all)
Start:  AIC=476.59
y ~ x + f

       Df Deviance    AIC
- f     1   84.993 474.77
<none>      84.808 476.59
- x     1   89.475 479.25

Step:  AIC=474.77
y ~ x

       Df Deviance    AIC
<none>      84.993 474.77
- x     1   89.507 477.29
> 
> summary(d.res)

Call:
glm(formula = y ~ x, family = poisson, data = d)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.3679  -0.7348  -0.1775   0.6987   2.3760  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  1.29172    0.36369   3.552 0.000383 ***
x            0.07566    0.03560   2.125 0.033580 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 89.507  on 99  degrees of freedom
Residual deviance: 84.993  on 98  degrees of freedom
AIC: 474.77

Number of Fisher Scoring iterations: 4

・・・

この結果より、平均種子数の予測値を求める関数は { \lambda = \exp(1.29172 + 0.07566 x) } になります。

f:id:fits:20131215205033p:plain

(2) MCMCpack を使ったベイズ統計によるポアソン回帰1 （MCMCpoisson 関数）

MCMCpack でポアソン回帰を行うには MCMCpoisson 関数を使うのが簡単だと思います。 基本的に glm 関数と同様のモデル式とデータを指定するだけです。

ここでは y ~ x のモデルだけをポアソン回帰してみました。

なお、glm 関数では線形予測子のパラメータ { \beta_1 } と { \beta_2 } の最尤推定値を取得しますが、MCMCpoisson 関数では { \beta_1 } と { \beta_2 } のそれぞれの分布を取得する点が大きく異なります。

poissonMcmcPoisson.R
library(MCMCpack)

d <- read.csv('data3a.csv')

d.res <- MCMCpoisson(y ~ x, data = d)

summary(d.res)

png("poissonMcmcPoisson.png")
plot(d$x, d$y, col = c("red", "blue")[d$f])
xx <- seq(min(d$x), max(d$x), length = 10000)
lines(xx, exp(mean(d.res[,1]) + mean(d.res[,2]) * xx), col="green")
dev.off()

glm の時とは異なり、{ \beta_1 } と { \beta_2 } のそれぞれの算術平均値 （mean の結果） を使って、平均種子数の予測値を描画 （緑の線） しています。

実行

実行すると下記のような結果になります。

> R CMD BATCH poissonMcmcPoisson.R
実行結果 poissonMcmcPoisson.Rout
・・・
> summary(d.res)

Iterations = 1001:11000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 10000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

               Mean     SD Naive SE Time-series SE
(Intercept) 1.29315 0.3602 0.003602       0.010757
x           0.07547 0.0352 0.000352       0.001047

2. Quantiles for each variable:

                2.5%     25%    50%     75%  97.5%
(Intercept) 0.586606 1.05208 1.2956 1.54153 1.9917
x           0.006025 0.05125 0.0756 0.09951 0.1438

・・・

Mean の値が glm の結果とほぼ同じ値になっています。

f:id:fits:20131215205111p:plain

ちなみに、d.res には MCMC サンプリング結果として 1万個の { \beta_1 } と { \beta_2 } の値が格納されています。

{ \beta_1 } の分布

f:id:fits:20131215205124p:plain

{ \beta_2 } の分布

f:id:fits:20131215210052p:plain

(3) MCMCpack を使ったベイズ統計によるポアソン回帰2 （MCMCmetrop1R 関数）

最後に MCMCmetrop1R 関数を使ってみます。

MCMCmetrop1R 関数では自前で定義した尤度関数を使って MCMC サンプリングを実施できるので汎用的に使えます。

今回のモデルの対数尤度は { \log L ( \beta_1, \beta_2 ) = \sum_i \log \frac{\lambda_i ^ y_i \exp(- \lambda_i)}{y_i!} } で、{ \lambda_i = \exp (\beta_1 + \beta_2 x_i) } なので、これらを関数化 （下記の func） しています。

尤度関数が対数尤度か否かは logfun で指定するようになっており （TRUE の場合が対数尤度）、デフォルト値は TRUE となっています。 （今回は対数尤度なので logfun = TRUE です）

また、theta.init で { \beta_1 } と { \beta_2 } の初期値を c(0, 0) と指定しています。

poissonMcmcMetrop.R
library(MCMCpack)

d <- read.csv('data3a.csv')

# 尤度関数（対数尤度関数）
func <- function(beta, x, y) {
  lambda <- exp(beta[1] + beta[2] * x)
  sum(log(dpois(y, lambda)))
}

d.res <- MCMCmetrop1R(func, theta.init = c(0, 0), x = d$x, y = d$y, logfun = TRUE)
# 下記でも同じ
#d.res <- MCMCmetrop1R(func, theta.init = c(0, 0), x = d$x, y = d$y)

summary(d.res)

png("poissonMcmcMetrop.png")
plot(d$x, d$y, col = c("red", "blue")[d$f])
xx <- seq(min(d$x), max(d$x), length = 10000)
lines(xx, exp(mean(d.res[,1]) + mean(d.res[,2]) * xx), col="green")
dev.off()

実行

実行すると下記のような結果になります。

> R CMD BATCH poissonMcmcMetrop.R
実行結果 poissonMcmcMetrop.Rout
・・・
> summary(d.res)

Iterations = 501:20500
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 20000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

        Mean      SD  Naive SE Time-series SE
[1,] 1.29545 0.35797 0.0025313      0.0077477
[2,] 0.07528 0.03498 0.0002473      0.0007561

2. Quantiles for each variable:

         2.5%     25%     50%     75%  97.5%
var1 0.584646 1.05368 1.29884 1.53850 1.9830
var2 0.007798 0.05159 0.07518 0.09907 0.1448

・・・

glm や MCMCpoisson と似たような結果となりました。

f:id:fits:20131215205205p:plain

なお、MCMCpoisson と MCMCmetrop1R ではバーンイン burnin とサンプリング数 mcmc のデフォルト値が異なっています。

関数 burnin のデフォルト値 mcmc のデフォルト値
MCMCpoisson 1000 10000
MCMCmetrop1R 500 20000

Groovy で Dempsy を分散実行

Groovy CEP

前回 と 前々回に続き、今回は Dempsy を Groovy で分散実行してみます。

サンプルソースは http://github.com/fits/try_samples/tree/master/blog/20131208/

はじめに

Dempsy を分散実行するには ZooKeeper が必要となりますので 「GroovyでZooKeeperを組み込み実行」 で作成した Groovy スクリプトを使う事にします。

また、分散実行するには Dempsy オブジェクトの構成を下記のようにします。 （※ の箇所が単独実行と異なる部分）

Dempsy 分散実行時の構成
プロパティ名 クラス
clusterSessionFactory ZookeeperSessionFactory ※
clusterCheck SpecificClusterCheck ※
defaultRoutingStrategy DecentralizedRoutingStrategy
defaultSerializer KryoSerializer
defaultStatsCollectorFactory StatsCollectorFactoryCoda
defaultTransport TcpTransport ※

Dempsy 分散実行アプリケーション

それでは、前回の KeySource 指定無しのサンプルを元に分散実行への対応を行ってみます。

基本的には Dempsy オブジェクトの構成を変えるだけです。

単独実行と異なり、分散実行では ClusterDefinition 毎にプロセスを起動しますので、実行時引数でクラスター名を指定するようにしています。

なお、dempsy.waitToBeStopped() はクラスター adaptor が終了しないようにするための措置です。 （dempsy.waitToBeStopped() が無いと、すぐに終了してしまいます）

money_count_vertx_cluster.groovy
・・・前回と同じ実装・・・

// 実行する ClusterDefinition の名称
def cluster = args[0]

def dempsy = new Dempsy()

dempsy.applicationDefinitions = [app]
dempsy.clusterSessionFactory = new ZookeeperSessionFactory('localhost:2181', 5000)
// 実行する ClusterDefinition を選定
dempsy.clusterCheck = new SpecificClusterCheck(new ClusterId('money-count', cluster))
// mp を 3ノード構成で処理する設定
dempsy.defaultRoutingStrategy = new DecentralizedRoutingStrategy(5, 3)
dempsy.defaultSerializer = new KryoSerializer()
dempsy.defaultStatsCollectorFactory = new StatsCollectorFactoryCoda()
dempsy.defaultTransport = new TcpTransport()

dempsy.start()

Runtime.runtime.addShutdownHook { ->
    println 'shutdown ...'
    dempsy.stop()
}

// Adaptor が終了するのを防止（MessageProcessor の場合は不要）
dempsy.waitToBeStopped()

MessageProcessor （クラスター mp） の実行ノード数などは DecentralizedRoutingStrategy で設定します。

DecentralizedRoutingStrategy のコンストラクタは DecentralizedRoutingStrategy(int defaultTotalSlots, int defaultNumNodes) となっており、第1引数でスロット（シャード）数、第2引数でノード数を指定するようになっています。

キー毎の担当スロットが messageKey.hashCode() % defaultTotalSlots で算出され、ノード毎にスロットが割り当てられて処理されます。

今回のサンプルのキー構成でスロットが適切にバラけるように defaultTotalSlots を 5 にしました。

スロット数 5 の場合のキー毎の messageKey.hashCode() % 5 結果
キー名 スロット
1 4
5 3
10 2
50 1
100 0
500 4
1000 3
2000 4
5000 2
10000 1

分散実行

それでは分散実行してみます。

まずは ZooKeeper を実行しておきます。

ZooKeeper 実行
> groovy zk_run.groovy zoo.cfg

次に Adaptor を実行します。

Adaptor 実行
> groovy money_count_vertx_cluster.groovy adaptor
・・・
MoneyAdaptor.start ...
・・・

そして MessageProcessor を実行します。

DecentralizedRoutingStrategy で指定したのは 3ノードですが、ノードが停止した場合の挙動を確認するため今回は 4ノード分のプロセスを実行しました。

(1) mp 実行 （1つ目）
> groovy money_count_vertx_cluster.groovy mp
・・・
2013-12-08 19:42:11,187 [myid:] - INFO  [main-SendThread(localhost:2181):ClientCnxn$SendThread@738] - Session establishment complete on server localhost/127.0.0.1:2181, sessionid = 0x142d19fda340012, negotiated timeout = 6000
2013-12-08 19:42:15,664 [myid:] - INFO  [main:SimpleThreadPool@270] - Job execution threads will use class loader of thread: main
・・・
(2) mp 実行 （2つ目）
> groovy money_count_vertx_cluster.groovy mp
・・・
(3) mp 実行 （3つ目）
> groovy money_count_vertx_cluster.groovy mp
・・・
(4) mp 実行 （4つ目）
> groovy money_count_vertx_cluster.groovy mp
・・・

結果1 - 全キーへのアクセス

下記のように全キーの URL へアクセスしてカウントアップしてみます。

$ curl http://localhost:8080/1
$ curl http://localhost:8080/10
$ curl http://localhost:8080/100
$ curl http://localhost:8080/1000
$ curl http://localhost:8080/10000
$ curl http://localhost:8080/5000
$ curl http://localhost:8080/500
$ curl http://localhost:8080/50
$ curl http://localhost:8080/5
$ curl http://localhost:8080/2000

(1) ～ (4) の mp の状態と割り当てられたスロットは下記のようになりました。

なお、(4) のノードにはスロットが割り当てられておらずキーも送信されていない事を確認できました。

(1) mp の状態 （スロット 3）
> groovy money_count_vertx_cluster.groovy mp
・・・
MoneyCount.activation : 1000
key: 1000, count: 1
MoneyCount.activation : 5
key: 1000, count: 1
key: 5, count: 1
・・・
(2) mp の状態 （スロット 0, 1）
> groovy money_count_vertx_cluster.groovy mp
・・・
MoneyCount.activation : 100
MoneyCount.activation : 10000
key: 100, count: 1
key: 10000, count: 1
MoneyCount.activation : 50
key: 100, count: 1
key: 10000, count: 1
key: 50, count: 1
・・・
(3) mp の状態 （スロット 2, 4）
> groovy money_count_vertx_cluster.groovy mp
・・・
MoneyCount.activation : 1
MoneyCount.activation : 10
key: 1, count: 1
key: 10, count: 1
MoneyCount.activation : 5000
MoneyCount.activation : 500
key: 1, count: 1
key: 500, count: 1
key: 5000, count: 1
key: 10, count: 1
MoneyCount.activation : 2000
key: 1, count: 1
key: 500, count: 1
key: 2000, count: 1
key: 5000, count: 1
key: 10, count: 1
・・・
(4) mp の状態 （変化なし）
> groovy money_count_vertx_cluster.groovy mp
・・・

結果2 - (1) の mp を停止

(1) の mp を Ctrl + c でプロセス終了して、下記 URL へアクセスしてみます。

$ curl http://localhost:8080/1000
$ curl http://localhost:8080/5

(1) の担当スロットが、代わりに (4) へ割り当てられてキーが送信されてくるようになります。

(4) mp の状態 （(1) の停止後）
> groovy money_count_vertx_cluster.groovy mp
・・・
MoneyCount.activation : 1000
MoneyCount.activation : 5
key: 1000, count: 1
key: 5, count: 1
・・・

Groovy で Dempsy を単独実行2 - KeySource

Groovy CEP

前回に続き、Dempsy を Groovy で単独実行してみます。 今回は KeySource の設定有無でどのように挙動が変わるかを簡単に調べます。

サンプルソースは http://github.com/fits/try_samples/tree/master/blog/20131201/

KeySource 無しの場合

まず、前回サンプルの Adaptor クラスを変更し、http://localhost:8080/<金額> へのアクセスがあった際に Dispatcher へメッセージを dispatch() するように変更しました。 （HTTP サーバー処理は Vert.x を組み込み実行）

また、outputExecuter に RelativeOutputSchedule を設定し 10秒毎にカウント値を出力するように変更しています。

money_count_vertx.groovy
・・・
class Constants {
    final static def MONEYS = [
        '1', '5', '10', '50', '100', '500', '1000', '2000', '5000', '10000'
    ]
}

class MoneyAdaptor implements Adaptor {
    Dispatcher dispatcher

    void start() {
        println 'MoneyAdaptor.start ...'

        def vertx = Vertx.newVertx()

        def rm = new RouteMatcher()
        rm.get '/:money', { req ->
            def money = req.params['money']

            if (Constants.MONEYS.contains(money)) {
                dispatcher.dispatch(new Money(req.params['money']))
            }
            req.response.end()
        }

        vertx.createHttpServer().requestHandler(rm.asClosure()).listen 8080
    }

    void stop() {
        println 'MoneyAdaptor.stop'
    }
}

def mp = new ClusterDefinition('mp', new MoneyCount())
mp.outputExecuter = new com.nokia.dempsy.output.RelativeOutputSchedule(10, java.util.concurrent.TimeUnit.SECONDS)

def app = new ApplicationDefinition('money-count').add(
    new ClusterDefinition('adaptor', new MoneyAdaptor()),
    mp
)
・・・

実行

それでは実行します。

出力1
> groovy money_count_vertx.groovy
・・・
MoneyAdaptor.start ...

MoneyAdaptor が start しただけで、MoneyCount は 1つもアクティブ化されていない状態です。

ここで http://localhost:8080/5 へアクセスすると、キー値 5 を処理する MoneyCount がアクティブ化され、10秒毎に現在のカウント値が出力されるようになります。

出力2
MoneyCount.activation : 5
key: 5, count: 1

この後、http://localhost:8080/5 へアクセスするたびにカウントアップされます。

次に、http://localhost:8080/10000 へアクセスするとキー値 10000 を処理する新しい MoneyCount がアクティブ化されます。

出力3
MoneyCount.activation : 10000
key: 10000, count: 1
key: 5, count: 4

このように、KeySource を設定しなかった場合は新しくキーが発生する度に MessageProcessor の新しいインスタンスがアクティブ化されます。

KeySource 有りの場合

次に、KeySource を設定してみます。 KeySource を設定するには下記のようにします。

  • KeySource インターフェースの実装オブジェクトを MessageProcessor を設定した ClusterDefinition へ setKeySource() で設定

なお、KeySource インターフェースは getAllPossibleKeys() で全てのキーを返すように実装します。

変更点は下記のようになります。

money_count_vertx_ks.groovy
・・・
def mp = new ClusterDefinition('mp', new MoneyCount()).setKeySource({
    Constants.MONEYS
} as KeySource<String>)
・・・

実行

それでは実行します。

出力1
> groovy money_count_vertx_ks.groovy
・・・
MoneyAdaptor.start ...
MoneyCount.activation : 1
MoneyCount.activation : 5
MoneyCount.activation : 10
MoneyCount.activation : 50
MoneyCount.activation : 100
MoneyCount.activation : 500
MoneyCount.activation : 1000
MoneyCount.activation : 2000
MoneyCount.activation : 5000
MoneyCount.activation : 10000

KeySource を設定しなかった場合とは異なり、起動時に全キーの MoneyCount がアクティブ化されます。

http://localhost:8080/5 等へアクセスすると該当キーがカウントアップされます。

出力2
key: 1, count: 0
key: 100, count: 0
key: 10000, count: 1
key: 500, count: 0
key: 1000, count: 0
key: 2000, count: 0
key: 5, count: 4
key: 5000, count: 0
key: 10, count: 0
key: 50, count: 0

まとめ

簡単にまとめると下記のようになります。

KeySource の設定 MessageProcessor のアクティブ化タイミング
無 新しくキーが発生する度
有 起動時