Julia と R のコラボレーションを jupyter notebook で記録

jupyter notebook というのは，いろいろな言語に対応している。

R の RStudio や，これまた多くの言語に対応している Atom と少し違って，スクリプトと実行結果をまとめて記録でき，成果物を LaTeX, PDF をはじめとして様々な形式で保存・閲覧できるという点で気に入っている。

たとえば，RStudio は，プログラム編集とコンソールとPlotなどのペインが分かれているので，作業しているときに煩わしい（別のやり方があるのかも知れないが）。その点，jupyter notebook はやったとおり，できたとおりのものが順に表示されていくので心地よい。あとで，不必要な部分や，追加処理を挿入するなどしてまとめて全部再度実行するのも簡単にできる。

 

それはさておき，今回は，Julia と R の環境を行き来して，いいとこ取りの分析環境を構築できるかどうか，ちょっとだけつまみ食いしてみよう。

以下では，jupiter notebook の結果を text ファイルに出力したまま（都合で html ファイルの一部分を切り取って）表示しているので不細工だが，本当は遙かに綺麗（そのまま出版できるレベル）なので，誤解なきよう。

 

まずは，Julia で，独立二標本の平均値の差の検定を行ってみる。
使うのは，ちょっと，長い名前だが UnequalVarianceTTest() だ。

 

In [1]:

using HypothesisTests

In [2]:

x = [2,1,3,2,4,3,2,3,4,5];

In [3]:

y = [3,2,3,4,5,6,3,2,3,3,3,4,3,2,6];

In [4]:

UnequalVarianceTTest(x, y)

Out[4]:

Two sample t-test (unequal variance)
------------------------------------
Population details:
    parameter of interest:   Mean difference
    value under h_0:         0
    point estimate:          -0.566667
    95% confidence interval: (-1.6209, 0.4876)

Test summary:
    outcome with 95% confidence: fail to reject h_0
    two-sided p-value:           0.2760

Details:
    number of observations:   [10,15]
    t-statistic:              -1.1192102478745312
    degrees of freedom:       20.56776717709893
    empirical standard error: 0.506309397848002

次に，Julia から RCall を介して，R の t.test() で同じデータを分析してみる

 

In [5]:

using RCall

長いスクリプトは，R""" と """ で挟んで，何行でも書ける。
$x, $y というのは Julia の x, y ということを表す。
結果を明示的に print() で書かないと，最後の出力しか表示されないので注意。

In [6]:

R"""
options(digits=15)
a = t.test($x, $y)
print(a)
print(a$statistic)
print(a$p.value)
"""

 

	Welch Two Sample t-test

data:  `#JL`$x and `#JL`$y
t = -1.119210247875, df = 20.5677671771, p-value = 0.275951369175
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.620943659272958  0.487610325939624
sample estimates:
       mean of x        mean of y 
2.90000000000000 3.46666666666667 

                t 
-1.11921024787453 
[1] 0.275951369174595

Out[6]:

RObject{RealSxp}
[1] 0.275951369174595

In [7]:

R"a"

Out[7]:

RObject{VecSxp}

	Welch Two Sample t-test

data:  `#JL`$x and `#JL`$y
t = -1.119210247875, df = 20.5677671771, p-value = 0.275951369175
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.620943659272958  0.487610325939624
sample estimates:
       mean of x        mean of y 
2.90000000000000 3.46666666666667 

図を描くと，図もその場に出力されるので快適。
ただ，希望しない出力もしてくれるのは大きなお世話ともいえるが，普通は気にしない分析結果なので，統計解析もしっかりやらなければならないよという戒めにもなるのかな？

In [8]:

R"""
group = rep(factor(c("x", "y")), c(length($x),length($y)))
value = c($x, $y)
a = boxplot(value ~ group)
"""

Out[8]:

RObject{VecSxp}
$stats
     [,1] [,2]
[1,]    1    2
[2,]    2    3
[3,]    3    3
[4,]    4    4
[5,]    5    5
attr(,"class")
        x 
"integer" 

$n
[1] 10 15

$conf
                 [,1]             [,2]
[1,] 2.00072025938679 2.59204575419949
[2,] 3.99927974061321 3.40795424580051

$out
[1] 6 6

$group
[1] 2 2

$names
[1] "x" "y"