Triad sou.

Log-gamma distribution

A log-gamma random variable, $Y$, is defined as,
\[
Y=\exp(-X)
\] where $X \in [0, \infty)$ follows a gamma distribution, and $Y \in (0, 1]$.


The probability density function of the log-gamma distribution is
\[
f_Y(y)=\frac{b^a\left\{\log(1/y) \right\}^{a-1}y^{b-1}}{\Gamma(a)},
\] where $b > 0$, and $a > 0$.


The expected value of $Y$ is
\[
\mathrm{E}[Y]=\int_0^1 yf_Y(y) \mathrm{d}y=\left( \frac{b}{1+b} \right )^a.
\] The variance of $Y$ is
\[
\mathrm{Var}[Y]=\mathrm{E}[Y^2]-(\mathrm{E}[Y])^2= \left( \frac{b}{2+b} \right)^a - \left( \frac{b}{1+b} \right)^{2a}.
\]


R code:

# shape = a, rate = b
a <- 4
b <- 7
r <- -rgamma(100000, shape = a, rate = b)
mean(exp(r))
(b/(1 + b))^a

mean(exp(r - log((b/(1 + b))^a)))

var(exp(r))
(b/(2 + b))^a - (b/(1 + b))^(2*a)


Result:

> mean(exp(r))
[1] 0.5855565
> (b/(1 + b))^a
[1] 0.5861816
> 
> mean(exp(r - log((b/(1 + b))^a)))
[1] 0.9989336
> 
> var(exp(r))
[1] 0.02245031
> (b/(2 + b))^a - (b/(1 + b))^(2*a)
[1] 0.0223414

Note

\[
\int_0^1 yf_Y(y) \mathrm{d}y= \frac{b^a}{\Gamma(a)} \int_0^1 \left\{\log(1/y) \right\}^{a-1}y^b \mathrm{d}y
\] \[
\begin{align*} z &= \log(1/y) \\
\mathrm{d}z &= -\frac{1}{y} \mathrm{d}y \\
\Re(z) &> 0, \log(1/0)=\infty, \log(1/1)=0.
\end{align*}
\] \[
\begin{align*}
\int_0^1 yf_Y(y) \mathrm{d}y &= \frac{b^a}{\Gamma(a)} \int_{\infty}^0 -z^{a-1} e^{-z(1+b)} \mathrm{d}z \\
&= \frac{b^a}{\Gamma(a)} \frac{\Gamma(a)}{(1+b)^a} \\
&= \left( \frac{b}{1+b} \right)^a \end{align*}.
\]

Java Compiler API を使って何かしてみよう

Java のバイナリから Java ソースコードコンパイルするための API を使って、String 型の文字列をコンパイルし、クラスから static なメソッドを呼んで値を返す方法を模索していました。


以下のプログラムと解説を参考にしています

  1. Compiling from Memory : Java Compiler tools « JDK 6 « Java
  2. Javaリフレクションメモ(Hishidama's Java Reflection Memo)
import java.io.IOException;
import java.lang.reflect.InvocationTargetException;
import java.lang.reflect.Method;
import java.net.URI;
import java.util.Arrays;

import javax.tools.Diagnostic;
import javax.tools.DiagnosticCollector;
import javax.tools.JavaCompiler;
import javax.tools.JavaFileObject;
import javax.tools.SimpleJavaFileObject;
import javax.tools.ToolProvider;
import javax.tools.JavaCompiler.CompilationTask;

public class CompileSourceInMemory {

  public static void main(String args[]) throws IOException {

    JavaCompiler compiler = ToolProvider.getSystemJavaCompiler();
    DiagnosticCollector<JavaFileObject> diagnostics =
      new DiagnosticCollector<JavaFileObject>();

    String src =
      "public class HelloWorld {" +
      "public static int main(String args[]) {" +
      "System.out.println(\"This is in another java file\");" +
      "return(1);" +
      "}" +
      "}";

    JavaFileObject file = new JavaSourceFromString("HelloWorld", src);

    String[] compileOptions = new String[] {"-d", "bin"} ;
    Iterable<String> compilationOptionss = Arrays.asList(compileOptions);

    Iterable <? extends JavaFileObject > compilationUnits = Arrays.asList(file);
    CompilationTask task = compiler.getTask(
                             null,
                             null,
                             diagnostics,
                             compilationOptionss,
                             null,
                             compilationUnits
                           );

    boolean success = task.call();
    for (Diagnostic<?> diagnostic : diagnostics.getDiagnostics()) {
      System.out.println(diagnostic.getCode());
      System.out.println(diagnostic.getKind());
      System.out.println(diagnostic.getPosition());
      System.out.println(diagnostic.getStartPosition());
      System.out.println(diagnostic.getEndPosition());
      System.out.println(diagnostic.getSource());
      System.out.println(diagnostic.getMessage(null));

    }
    System.out.println("Success: " + success);

    if (success) {
      Object ret;
      try {
        Class<?> clazz = Class.forName("HelloWorld");
        Method method = clazz.getMethod("main", new Class[] { String[].class });
        ret = method.invoke(null, new Object[] { null });
        System.out.println(ret);
      } catch (ClassNotFoundException e) {
        System.err.println("Class not found: " + e);
      } catch (NoSuchMethodException e) {
        System.err.println("No such method: " + e);
      } catch (IllegalAccessException e) {
        System.err.println("Illegal access: " + e);
      } catch (InvocationTargetException e) {
        System.err.println("Invocation target: " + e);
      }
    }
  }
}

class JavaSourceFromString extends SimpleJavaFileObject {

  final String code;

  JavaSourceFromString(String name, String code) {
    super(URI.create("string:///" + name.replace('.', '/') +
                     Kind.SOURCE.extension), Kind.SOURCE);
    this.code = code;
  }

  @Override
  public CharSequence getCharContent(boolean ignoreEncodingErrors) {
    return code;
  }

}

ポイント

私の環境では Eclipse 上で実行しているので、

String[] compileOptions = new String[] {"-d", "bin"} ;

を入れています。


もう一つは、

Object ret;
try {
  Class<?> clazz = Class.forName("HelloWorld");
  Method method = clazz.getMethod("main", new Class[] { String[].class });
  ret = method.invoke(null, new Object[] { null });
  System.out.println(ret);
}

という感じにリフレクションを使って、返り値がある場合に拡張してみました。
とりあえず Eclipse 上では実行できた。ので、jar ファイルなどにしても動けばいいのですが・・・

2012/06/25 追記

実行可能 jar ファイル化して動かなかった原因は、jdk ではないと動かない (tools.jar 必須) のに、jre から動いていたため・・・
最新 jdk 以外を全てアンインストール / jar ファイルの関連づけを変更したら普通に動きました。
とりあえず、

JavaCompiler compiler = ToolProvider.getSystemJavaCompiler();
if (compiler == null) {
  System.out.println("null dayo!");
}

で、compiler を取ってこれているかどうかはチェックはしておいた方が良いでしょう。
null の場合は tools.jar が読み込めていないので、jdk から実行できていそうかなどをちゃんとチェックしておくとよいみたいです。

2012/06/26 リフレクションメモ

複数引数のリフレクションサンプルを作ってみました。
可変長引数がないので、以下の様に列挙するだけで簡単に使えますね。

import java.lang.reflect.InvocationTargetException;
import java.lang.reflect.Method;

public class a {
  public static void main(String[] args) {
    Class<?> clazz;
    try {
      clazz = Class.forName("test");
      Object obj = clazz.newInstance();
      Method method =
        obj.getClass().getMethod("fx", double.class, double.class, double[].class);
      System.out.println((Double) method.invoke(null, 2.3, 1.0, new double[] {1.0, 2.0, 1.0}));
      // 7.3 が出力される
    } catch (ClassNotFoundException e) {
      e.printStackTrace();
    } catch (InstantiationException e) {
      e.printStackTrace();
    } catch (IllegalAccessException e) {
      e.printStackTrace();
    } catch (SecurityException e) {
      e.printStackTrace();
    } catch (NoSuchMethodException e) {
      e.printStackTrace();
    } catch (IllegalArgumentException e) {
      e.printStackTrace();
    } catch (InvocationTargetException e) {
      e.printStackTrace();
    }
  }
}

class test {
  public static double fx(double a, double b, double[] c) {
    double cc = 0;
    for (int i = 0; i < c.length; i++) {
      cc += c[i];
    }
    return a + b + cc;
  }
}

R Commander による Kaplan-Meier プロットの作図手順

  1. R のインストール
  2. バージョン 2.15.0 以上の R をインストールします。ソフトウェアのインストールをした経験がある方なら特に困ることはないでしょう。


  3. R Commander のインストール

  4. R を起動し、以下のコマンドを実行するか、パッケージのインストール画面等から "Rcmdr" を選択します。

    install.packages("Rcmdr")
    

    • パッケージインストール画面の表示


    • ミラーサーバーを選択 (なるべく近そうなサーバを選ぶと良いでしょう)


    • インストールするパッケージの選択


    インストールエラーが表示された場合は、もう一度同じ操作を行ってみてください。


  5. R Commander プラグインのインストール

  6. 同様にプラグインをインストールします。
    以下のコマンドを実行するか、パッケージのインストール画面等から "RcmdrPlugin.KMggplot2" を選択します。上の画像のように、複数のパッケージを同時にインストールすることもできます。

    install.packages("RcmdrPlugin.KMggplot2")
    

  7. R Commander の起動

  8. 以下のコマンドを実行するか、パッケージのロード画面から "Rcmdr" を選択します。

    library("Rcmdr")
    


    RcmdrPlugin.KMggplot2 0.2-0 以降では、R Commander の起動とプラグインのロードを同時に行うことができます。
    以下のコマンドを実行するか、パッケージのロード画面から "RcmdrPlugin.KMggplot2" を選択します。

    library("RcmdrPlugin.KMggplot2")
    

    ロードが終了したら、5 番目のステップ (プラグインの読み込み) を実行する必要はありません。



    • パッケージのロード


    • 初回起動時は、以下の様な画面が表示されるはずなので、両方共そのまま [はい] を選択します。選択すると必要なパッケージが自動的にインストールされます




    • Rcmdr の起動画面



  9. R Commander プラグインの読み込み

  10. グラフ作成に必要なプラグインを読み込みます。


    • プラグインのロード




    • 再起動を促されるので再起動します



  11. データの読み込み

    • イベント変数をイベントありの場合が 1、打ち切りの場合が 0 でコーディングした形式のデータを準備し、クリップボードにコピーしておきます




    • データ、データのインポート、テキストファイル・・・を選択します


    • データセット名を指定し、変数名はあり、クリップボード、タブをチェックし、OK をクリックします



  12. 作図

    • KMggplot2、Kaplan-Meier プロット... を選択します


    • 時間変数、イベント変数、必要に応じてその他のオプションを設定し、OK をクリックします。グラフをファイルにセーブしたい場合はグラフの保存にチェックを入れておきます


    • 実行結果です


記号微分

deriv を利用して、R で記号微分してみます。

正規分布の対数尤度関数を記号微分してみよう

正規分布の対数尤度関数は、
\[
l(\mu, \s; x)=-\frac{\log(2\pi\s^2)}{2}-\frac{(x - \mu)^2}{2\s^2}
\] ですので、これを偏微分すると、
\[
\begin{align*}
\frac{\partial l(\mu, \s; x)}{\partial \mu} &= \frac{x-\mu}{\s^2} \\
\frac{\partial l(\mu, \s; x)}{\partial \s} &= -\frac{1}{\s}+\frac{(x-\mu)^2}{\s^3}
\end{align*}
\] となる予定です。


deriv を用いた数値微分では、まず被微分関数を expressioncall、または formula の形で定義します。

fx <- expression(
  -log(2*pi*sd^2)/2 - (x - mean)^2/(2*sd^2)
)


あとは、

  • 定義した被微分関数
  • 微分したい変数名
  • 出力される微分された関数に対する引数名
  • 必要に応じてオプション

を指定して実行するだけです。

dfx <- deriv(
  fx,
  c("mean", "sd"),
  function.arg = c("x", "mean", "sd")
)
dfx


結果は以下のようになりました。

> dfx
function (x, mean, sd) 
{
    .expr1 <- 2 * pi
    .expr2 <- sd^2
    .expr3 <- .expr1 * .expr2
    .expr7 <- x - mean
    .expr8 <- .expr7^2
    .expr9 <- 2 * .expr2
    .expr14 <- 2 * sd
    .value <- -log(.expr3)/2 - .expr8/.expr9
    .grad <- array(0, c(length(.value), 2L), list(NULL, c("mean", 
        "sd")))
    .grad[, "mean"] <- 2 * .expr7/.expr9
    .grad[, "sd"] <- -(.expr1 * .expr14/.expr3/2 - .expr8 * (2 * 
        .expr14)/.expr9^2)
    attr(.value, "gradient") <- .grad
    .value
}


一応確認しておくと、

.expr2 <- sd^2
.expr7 <- x - mean
.expr9 <- 2 * .expr2
.grad[, "mean"] <- 2 * .expr7/.expr9

# 2 * .expr7/.expr9 =
#   2 * (x - mean) / 2 * sd^2 =
#   (x - mean) / sd^2

および、

.expr1 <- 2 * pi
.expr2 <- sd^2
.expr3 <- .expr1 * .expr2
.expr7 <- x - mean
.expr8 <- .expr7^2
.expr9 <- 2 * .expr2
.expr14 <- 2 * sd
.grad[, "sd"] <- -(.expr1 * .expr14/.expr3/2 - .expr8 * (2 * 
        .expr14)/.expr9^2)

# -.expr1 * .expr14/.expr3/2 =
#   -2 * pi * 2 * sd / (2 * pi * sd^2) / 2 =
#   -1/sd
# .expr8 * (2 * .expr14)/.expr9^2 =
#   (x - mean)^2 * (2 * 2 * sd) / (2 * sd^2)^2 =
#   (x - mean)^2 / sd^3

となります。
このぐらいなら手計算の方が早いですが、いろいろな場面で利用できるのではないでしょうか。


また、Hessian を自動で出してくれる deriv3 もあります。

注意

deriv のヘルプのタイトルは、

Symbolic and Algorithmic Derivatives of Simple Expressions

となっていて、以下の説明があります。

The internal code knows about the arithmetic operators +, -, *, / and ^, and the single-variable functions exp, log, sin, cos, tan, sinh, cosh, sqrt, pnorm, dnorm, asin, acos, atan, gamma, lgamma, digamma and trigamma, as well as psigamma for one or two arguments (but derivative only with respect to the first). (Note that only the standard normal distribution is considered.)

deriv では基本的な演算子と、いくつかの単変数の関数ぐらいしか使えないため、被微分関数を作るときに注意が必要です。

RcmdrPlugin.KMggplot2_0.1-0 is on CRAN now

I posted a new version of the "RcmdrPlugin.KMggplot2" package, which is an Rcmdr plug-in for a "ggplot2" GUI front-end.
This package assists you to make "ggplot2" graphics.

NEWS

Changes in version 0.1-0 (2012-05-18)
theme_simple()

The theme_simple() is a simple black and white theme.
Top and right borders and both grids are removed from the theme_bw().

A list box for facet_xx() functions

A list box for facet_xx() functions was added.

Confidence intervals for Kaplan-Meier plots

Violin plots

Stacked bar charts

Univariate plots at diagonal positions (ggplot2::plotmatrix)

"..count.. + facet_grid" の挙動について

追記

ggplot-0.9.1 で修正されていました!
ggplot-0.9.1 での実行結果

メモ

スクリプトの書き方が良くないのかもしれませんが、..count..facet_gridを併用したときの挙動をメモっておきます。
バージョンは R-2.15.0 + ggplot-0.9.0 です。


こう書くと全ての度数が 1 になってしまう。

require("ggplot2")

mydf1 <- data.frame(
  x = sample(letters[1:2], 500, replace = T),
  y = sample(letters[10:11], 500, replace = T)
)

ggplot(mydf1, aes(x = x, y = ..count..)) +
  geom_bar() + facet_grid( ~ y) + opts(title = "..count.. + facet_grid")


facet_wrap の場合は意図したグラフが得られました。

ggplot(mydf1, aes(x = x, y = ..count..)) +
  geom_bar() + facet_wrap( ~ y) + opts(title = "..count.. + facet_wrap")


また、データフレームに id 変数を加えてみると facet_wrap の結果と一致しました。

mydf2 <- data.frame(
  mydf1,
  id = 1:500
)
ggplot(mydf2, aes(x = x, y = ..count..)) +
  geom_bar() + facet_grid( ~ y) + opts(title = "..count.. + facet_grid + id"))


何でだろう。
多分 ..count.. の挙動を良く理解していないからだろう。

行列のスカラー微分のメモ

$\mathbf{A}$, $\mathbf{B}$ について、スカラー $\theta$ で微分する場合を考える。
だたし、
\[
\mathbf{A}= \left\{ a_{ij} \right\}, \frac{\partial \mathbf{A}}{\partial \theta}= \left\{ \frac{\partial a_{ij}}{\partial \theta} \right\}
\] とし、$\mathbf{B}$ も同様に定義する。

行列積の微分

行列積 $\mathbf{A} \mathbf{B}$ が定義できる場合、
\[
\frac{\partial \mathbf{A}\mathbf{B}}{\partial \theta}= \frac{\partial \mathbf{A}}{\partial \theta}\mathbf{B} + \mathbf{A}\frac{\partial \mathbf{B}}{\partial \theta}
\]

逆行列微分

$\mathbf{A}$ が正方行列で逆行列が存在するならば、
\[
\frac{\partial \mathbf{A}^{-1}}{\partial \theta}= -\mathbf{A}^{-1} \frac{\partial \mathbf{A}}{\partial \theta} \mathbf{A}^{-1}
\]

逆行列微分の導出

\[
\begin{align*}
\frac{\partial \mathbf{A}\mathbf{A}^{-1}}{\partial \theta} &=
\frac{\partial \mathbf{I}}{\partial \theta} \\
\frac{\partial \mathbf{A}}{\partial \theta}\mathbf{A}^{-1} + \mathbf{A}\frac{\partial \mathbf{A}^{-1}}{\partial \theta} &=
\mathbf{0}\\
\mathbf{A}\frac{\partial \mathbf{A}^{-1}}{\partial \theta} &= -\frac{\partial \mathbf{A}}{\partial \theta}\mathbf{A}^{-1} \\
\frac{\partial \mathbf{A}^{-1}}{\partial \theta} &= -\mathbf{A}^{-1}\frac{\partial \mathbf{A}}{\partial \theta}\mathbf{A}^{-1}
\end{align*}
\]

正方行列 $\mathbf{A}$, $\mathbf{B}$ について、
\[
\begin{align*}
\frac{\partial \mathbf{A}^{-1}\mathbf{B}\mathbf{A}^{-1}}{\partial \theta} &=
\frac{\partial \mathbf{A}^{-1}}{\partial \theta}\mathbf{B}\mathbf{A}^{-1} + \mathbf{A}^{-1}\frac{\partial \mathbf{B}\mathbf{A}^{-1}}{\partial \theta} \\ &= -\mathbf{A}^{-1}\frac{\partial \mathbf{A}}{\partial \theta}\mathbf{A}^{-1}\mathbf{B}\mathbf{A}^{-1} + \mathbf{A}^{-1}\frac{\partial \mathbf{B}}{\partial \theta}\mathbf{A}^{-1} - \mathbf{A}^{-1}\mathbf{B}\mathbf{A}^{-1}\frac{\partial \mathbf{A}}{\partial \theta}\mathbf{A}^{-1} \\ &=
\mathbf{A}^{-1} \left( \frac{\partial \mathbf{B}}{\partial \theta} -\frac{\partial \mathbf{A}}{\partial \theta}\mathbf{A}^{-1}\mathbf{B} -\mathbf{B}\mathbf{A}^{-1}\frac{\partial \mathbf{A}}{\partial \theta} \right) \mathbf{A}^{-1}
\end{align*}
\]

$\mathbf{A}$ の要素が全て定数の場合、$\partial \mathbf{A} / \partial \theta = \mathbf{0}$ より $\mathbf{A}^{-1} (\partial \mathbf{B} / \partial \theta) \mathbf{A}^{-1}$。