Homebrew経由でRをインストール

Homebrew経由でRをインストールした際のメモです.

ターミナルから,以下のコマンドを入力します.

$ brew tap brewsci/science
以前は,外部リポジトリの `homebrew/science` に `tap` すれば良かったようですが,外部リポジトリ `brewsci/science` に `tap` しなければならなくなっているようです.
brew tapでscienceリポジトリを追加します*.

その後に `R` をインストールします.

$ brew install r
これで,無事 `R` がインストールされます.

以下のようにして,`R`を起動します.

$ R
 
無事,インストールが成功していれば,以下のような画面が表示されます.
R version 3.6.0 (2019-04-26) -- "Planting of a Tree"
Copyright (C) 2019 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin17.7.0 (64-bit)
 
R は、自由なソフトウェアであり、「完全に無保証」です。 
一定の条件に従えば、自由にこれを再配布することができます。 
配布条件の詳細に関しては、'license()' あるいは 'licence()' と入力してください。 
 
  Natural language support but running in an English locale
 
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
 
'demo()' と入力すればデモをみることができます。 
'help()' とすればオンラインヘルプが出ます。 
'help.start()' で HTML ブラウザによるヘルプがみられます。 
'q()' と入力すれば R を終了します。 
 
 起動準備中です -  警告メッセージ: 
1: Setting LC_COLLATE failed, using "C" 
2: Setting LC_TIME failed, using "C" 
3: Setting LC_MESSAGES failed, using "C" 
4: Setting LC_MONETARY failed, using "C" 

さらに,RStanをインストールするには,Rが立ち上がっている状態で,以下のコマンドを入力します.
> install.packages("rstan")
 
すると,以下のようにインストールが始まります.
下記ののように,ミラーサイトの選択を要求されるので,Japan(42)を選択します.
 パッケージを ‘/usr/local/lib/R/3.6/site-library’ 中にインストールします 
 (‘lib’ が指定されていないため) 
 --- このセッションで使うために、CRAN のミラーサイトを選んでください --- 
Secure CRAN mirrors 
 
 1: 0-Cloud [https]                   2: Algeria [https]                
 3: Australia (Canberra) [https]      4: Australia (Melbourne 1) [https]
 5: Australia (Melbourne 2) [https]   6: Australia (Perth) [https]      
 7: Austria [https]                   8: Belgium (Ghent) [https]        
 9: Brazil (PR) [https]              10: Brazil (RJ) [https]            
11: Brazil (SP 1) [https]            12: Brazil (SP 2) [https]          
13: Bulgaria [https]                 14: Chile 1 [https]                
15: Chile 2 [https]                  16: China (Hong Kong) [https]      
17: China (Guangzhou) [https]        18: China (Lanzhou) [https]        
19: China (Shanghai 1) [https]       20: China (Shanghai 2) [https]     
21: Colombia (Cali) [https]          22: Czech Republic [https]         
23: Denmark [https]                  24: Ecuador (Cuenca) [https]       
25: Ecuador (Quito) [https]          26: Estonia [https]                
27: France (Lyon 1) [https]          28: France (Lyon 2) [https]        
29: France (Marseille) [https]       30: France (Montpellier) [https]   
31: France (Paris 2) [https]         32: Germany (Erlangen) [https]     
33: Germany (Göttingen) [https]      34: Germany (Münster) [https]      
35: Greece [https]                   36: Hungary [https]                
37: Iceland [https]                  38: Indonesia (Jakarta) [https]    
39: Ireland [https]                  40: Italy (Padua) [https]          
41: Japan (Tokyo) [https]            42: Japan (Yonezawa) [https]       
43: Korea (Busan) [https]            44: Korea (Gyeongsan-si) [https]   
45: Korea (Seoul 1) [https]          46: Korea (Ulsan) [https]          
47: Malaysia [https]                 48: Mexico (Mexico City) [https]   
49: Norway [https]                   50: Philippines [https]            
51: Serbia [https]                   52: Spain (A Coruña) [https]       
53: Spain (Madrid) [https]           54: Sweden [https]                 
55: Switzerland [https]              56: Turkey (Denizli) [https]       
57: Turkey (Mersin) [https]          58: UK (Bristol) [https]           
59: UK (London 1) [https]            60: USA (CA 1) [https]             
61: USA (IA) [https]                 62: USA (KS) [https]               
63: USA (MI 1) [https]               64: USA (OR) [https]               
65: USA (TN) [https]                 66: USA (TX 1) [https]             
67: Uruguay [https]                  68: (other mirrors)                
 
 
Selection: 42
すると,インストールが始まります(ズラズラと画面が流れて,完了までに結構時間がかかります).
 
Rへの読み込みは,以下のようにコマンドを入力します.
> library(rstan)
 
以下の設定を行うことが推奨されているようです.
> rstan_options(auto_write=TRUE)
> options(mc.cores=parallel::detectCores())
最初のコマンドは,コンパイルした結果を保存するもので,次のコマンドは複数のコアを使用するものです.
`R` を終了させるには,以下のコマンドを入力します.
> quit()
___
* `tap`を実行すると,以下のような画面が流れます.
Updating Homebrew...
==> Auto-updated Homebrew!
Updated 2 taps (homebrew/core and homebrew/cask).
==> New Formulae
...
==> Tapping brewsci/science
Cloning into '/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science'...
remote: Enumerating objects: 435, done.
remote: Counting objects: 100% (435/435), done.
remote: Compressing objects: 100% (431/431), done.
remote: Total 435 (delta 2), reused 181 (delta 2), pack-reused 0
Receiving objects: 100% (435/435), 398.22 KiB | 152.00 KiB/s, done.
Resolving deltas: 100% (2/2), done.
Tapped 416 formulae (455 files, 1.2MB).
`...` は実際に表示されず,略していることを表します.
 
** `R`のインストールが始まると,以下のような画面が流れます.
==> Installing dependencies for r: isl, mpfr, gcc, gettext, libpng, openblas, pcre and readline
...
==> readline
readline is keg-only, which means it was not symlinked into /usr/local,
because macOS provides the BSD libedit library, which shadows libreadline.
In order to prevent conflicts when programs look for libreadline we are
defaulting this GNU Readline installation to keg-only.
 
For compilers to find readline you may need to set:
  export LDFLAGS="-L/usr/local/opt/readline/lib"
  export CPPFLAGS="-I/usr/local/opt/readline/include"
 
For pkg-config to find readline you may need to set:
  export PKG_CONFIG_PATH="/usr/local/opt/readline/lib/pkgconfig"

`...` は実際に表示されず,略していることを表します.

Juliaの計算の仕組み

Juliaは,
「関数を実行したときに,与えられた引数の型情報を使い,その関数をネイティブコードにコンパイルしてから実行する仕組み」
なので,高速に計算したい場合は関数化して計算する必要があります.

以下に,円周率のモンテカルロ計算を実行した時にかかる時間の計測例を示します.
(オリジナルは,黒木玄さんの「JuliaとJupyterのすすめ」にあるコードを転用してREPLで確認)

円周率のモンテカルロ計算を関数化せずに実行する場合

 

julia> @time begin

       L = 10^7

       c = 0

       for i in 1:L

       global c

       c += ifelse(rand()^2+rand()^2 ≤ 1, 1, 0)

       end

       4c/L

       end

  1.940780 seconds (40.00 M allocations: 762.914 MiB, 1.05% gc time)

3.1415724


関数化して実行場合する場合

 

julia> function simpi(L=10^7)

       c = 0

       for i in 1:L

       c += ifelse(rand()^2+rand()^2 ≤ 1, 1, 0)

       end

       4c/L

       end

simpi (generic function with 2 methods)

 

julia> @time simpi(10^5)

  0.020499 seconds (25.11 k allocations: 1.288 MiB)

3.12924


実行環境にもよりますが,上記の例では処理速度が,数十倍も違うことが確認できルはずです.
なお,オリジナルの,黒木玄さんの「JuliaとJupyterのすすめ」には,gccと比較しても断然速いことが示されていいます.

ただし,gccで適切なライブラリを選定すると,Juliaより速く計算できる例も示されていますが,最適なライブラリ設定は面倒な(というよりわからない)事も多々あるので,Juliaで計算することのメリットを実感することができます.

上記の例のJupyter Notebookファイルは,GitHubMechanism of calculation.ipynbというファイルで見ることができます.

-Julia- 平均の計算

Juliaで平均の計算を行ってみた際のメモです.実行例はREPLで行なった結果です.

まずは相加平均(算術平均)についてです.

\mu = \frac{1}{n} \sum\limits_{i=1}^n x_i
1〜10までの整数の和の平均を求めてみます.
総和は`sum()`関数を用いて求めます.
julia> x1 = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
10-element Array{Float64,1}:
  1.0
  2.0
  3.0
  4.0
  5.0
  6.0
  7.0
  8.0
  9.0
 10.0
 
julia> n = length(x1)
10
 
julia> mu_a = sum(x1) / n
5.5
 
続いて幾何平均についてです.
( \prod\limits_{i=1}^{n} x_i ) ^ \frac {1} {n}
1〜10の幾何平均を求めてみます.
総乗は`prod()`関数を用いて求めます.
julia> x2 = [1.0, 2.0, 3.0, 4.0, 5.0]
5-element Array{Float64,1}:
 1.0
 2.0
 3.0
 4.0
 5.0
 
julia> n = length(x2)
5
 
julia> mu_g = prod(x2) ^ (1 / n)
2.605171084697352
 
Python3で確認のための計算を行なってみると以下のようになります.
>>> import scipy.stats.mstats as mstats
>>> mstats.gmean(x2)
2.6051710846973517

幾何平均(対数)- log-average -についてです.
\exp \left( \cfrac{1}{n} \sum_{i=1}^n \ln {x_i} \right)
この平均は,元のデータの数値 $a_{i}$ を対数*に変換して相加(算術)平均を求め、指数関数*を適用して元の数値の幾何平均を得ます.
log-averageを対数平均(logarithmic mean )と混同しないように注意が必要です.

julia> x3 = [1.0, 2.0, 3.0, 4.0, 5.0]
5-element Array{Float64,1}:
 1.0
 2.0
 3.0
 4.0
 5.0
 
julia> x4 = map(a -> log(a), x3)
5-element Array{Float64,1}:
 0.0               
 0.6931471805599453
 1.0986122886681098
 1.3862943611198906
 1.6094379124341003
 
julia> n = length(x3)
5
 
julia> mu_g_ln = exp*1
2.605171084697352
 
Python3で確認のための計算を行なってみると以下のようになります.
>>> import numpy as np
>>> from statistics import mean
>>> x3 = [1.0, 2.0, 3.0, 4.0, 5.0]
>>> np.exp(mean (np.log(x3)))
2.6051710846973517
 

なお,上記の例のJupyter Notebookファイルは,GitHubMean.ipynb というファイルで見ることができます.

* `ln`は自然対数で,`e`を底とする対数です.`exp`は`e`を底とする指数関数です.`exp`は`log`と逆の意味になるので,`log(x)`を`exp`すると元の値に戻し,`exp(x)`を`log`しても元に戻ります.

julia> exp(log(pi))
3.141592653589793
 
julia> log(exp(pi))
3.141592653589793
 

*1:1 / n) * sum(x4

-Python v Julia- 二項分布

以前に,Pythonで二項分布のグラフをプロットする投稿をしました.今回は,同様の内容をJuliaで行ってみます.
二項分布の確率質量関数は

 \mathrm{Bin}(m | M, \mu) = {}_M C_m \mu^{m} (1 - \mu)^{M - m}

です.
実装は以下のようになります.実装例はREPLで行った結果です.

julia> using Statistics
 
julia> using Distributions
 
julia> using StatsPlots
 
julia> M = 50
50
 
julia> mu = 0.5
0.5
 
julia> m = Binomial(M, mu)
Binomial{Float64}(n=50, p=0.5)
 
 
julia> scatter(m, leg=false, xlims = (0,50))
 
以下のようなグラフが描かれます.
この図を保存するには,以下のコマンドを実行します(この例ではファイル名をplot1として.png形式のファイルで保存)
julia> savefig("plot1.png")
 
また,参考までに平均と分散を求めると以下のようになります.
定義に基づいて計算した平均
julia> println("mean = ", M*mu)            
mean = 25.0
 
Statics.mean()関数を用いて計算した平均
julia> println("mean = ", Statistics.mean(m))
mean = 25.0
 
定義に基づいて計算した分散
julia> println("variance = ", M*mu*(1-mu))         
variance = 12.5
 
Statics.var()関数を用いて計算した平均
julia> println("variance = ", Statistics.var(m))   # Calculate using function
variance = 12.5
 
なお,上記の例のJupyter Notebookファイルは,GitHubBinomialDistribution.ipynb というファイルで見ることができます.

-Python vs Julia- 順列・組合せ

以前に,順列・組合せに関する計算をPythonで行う投稿 をしました.
今回は,それをJuliaで行ってみます.計算例はREPLで行なったものです.

順列
1, 2, 3の3つの数字の並べ方は 3! = 3 \times 2 \times 1 = 6 通りです.
階乗を求めるには factorial() 関数を使用します.

julia> factorial(3)
6
 
順列
{}_n P _r = \cfrac{n!}{(n - r)!}
を表示するには,以下のように permutations()  関数を使用します.
 

1, 2, 3, 4, 5 から3つを選ぶ並べ方は,{}_5 P _3=5 \times 4 \times 3 = 60 通りとなります.
何通りあるかだけを知りたい場合には,以下のようにして求めます.

julia> binomial(5, 3) * factorial(3)
60
 

順列の並びを表示させるには,permutations() 関数を用いて,第一引数に対象とする配列,第二引数に選ぶ数を指定します.

julia> using Combinatorics
 

julia> seq = (1, 2, 3)
(1, 2, 3)

julia> collect(permutations(seq))
6-element Array{Array{Int64,1},1}:
 [1, 2, 3]
 [1, 3, 2]
 [2, 1, 3]
 [2, 3, 1]
 [3, 1, 2]
 [3, 2, 1]
 
順列,組合せを扱うには,Combinatorics パッケージ*を使用します.
なお,collect は,配列を展開して表示するためのコマンドです.
 
上記のように特に配列(seq)を定義しなくても実行可能です.
julia> collect(permutations([1, 2, 3]))
6-element Array{Array{Int64,1},1}:
 [1, 2, 3]
 [1, 3, 2]
 [2, 1, 3]
 [2, 3, 1]
 [3, 1, 2]
 [3, 2, 1]
 
julia> collect(permutations([1, 2, 3, 4, 5], 3))
60-element Array{Array{Int64,1},1}:
 [1, 2, 3]
 [1, 2, 4]
 [1, 2, 5]
 [1, 3, 2]
 [1, 3, 4]
 [1, 3, 5]
 [1, 4, 2]
 [1, 4, 3]
 [1, 4, 5]
 [1, 5, 2]
 [1, 5, 3]
 [1, 5, 4]
 [2, 1, 3]
 [2, 1, 4]
 [2, 1, 5]
 [2, 3, 1]
 [2, 3, 4]
 [2, 3, 5]
 [2, 4, 1]
 [2, 4, 3]
 [2, 4, 5]
 ⋮        
 [4, 2, 1]
 [4, 2, 3]
 [4, 2, 5]
 [4, 3, 1]
 [4, 3, 2]
 [4, 3, 5]
 [4, 5, 1]
 [4, 5, 2]
 [4, 5, 3]
 [5, 1, 2]
 [5, 1, 3]
 [5, 1, 4]
 [5, 2, 1]
 [5, 2, 3]
 [5, 2, 4]
 [5, 3, 1]
 [5, 3, 2]
 [5, 3, 4]
 [5, 4, 1]
 [5, 4, 2]
 [5, 4, 3]


組合せ
1, 2, 3, 4, 5 から3つを選ぶ組合せは,{}_5 C _3= \frac{{}_5 P_3}{3!} = 10 通りとなります.
何通りあるかを知りたい場合には以下のようにします.

julia> binomial(5, 3)
10

組合せ
[tex:{}_n C _r= \cfrac{{}_n P_r}{r!} = \cfrac{n!}{r! (n - r)!}]

を表示するには,以下のように combinations() 関数を用いて,第一引数に対象とする配列,第二引数に選ぶ数を指定します.
julia> seq2 = (1, 2, 3, 4, 5)
(1, 2, 3, 4, 5)
 
julia> collect(combinations(seq2, 3))
10-element Array{Array{Int64,1},1}:
 [1, 2, 3]
 [1, 2, 4]
 [1, 2, 5]
 [1, 3, 4]
 [1, 3, 5]
 [1, 4, 5]
 [2, 3, 4]
 [2, 3, 5]
 [2, 4, 5]
 [3, 4, 5]
 
順列の例と同様に,配列を事前に定義せずに,第一引数に記載しても同様な結果が得られます.
 
なお,上記の例のJupyter Notebookファイルは,GitHubDiscreteProbabilityCombinationというファイルで見ることができます.


* Combinatorics パッケージを追加する際には,パッケージモードで " add Combinatorics "と入力します.すると以下のようにインストールが行われます.

(v1.1) pkg> add Combinatorics
 Resolving package versions...
 Installed Combinatorics ─ v0.7.0
 Installed Polynomials ─── v0.5.2
  Updating `~/.julia/environments/v1.1/Project.toml`
  [861a8166] + Combinatorics v0.7.0
  Updating `~/.julia/environments/v1.1/Manifest.toml`
  [861a8166] + Combinatorics v0.7.0
  [f27b6e38] + Polynomials v0.5.2
 

 

-Python vs Julia- 整数の乱数を発生させる

過去にPythonを用いて整数の乱数を発生させるコード例を投稿しました.

その中では,20回サイコロを振ることを想定して,整数の乱数を発生させるに,Pythonで以下のように書きました.

>>> import numpy as np
>>> np.random.choice(np.arange(1, 7), 20)
array([4, 2, 5, 1, 6, 5, 5, 3, 1, 5, 4, 4, 3, 4, 1, 5, 2, 1, 2, 2])


これをJuliaで行う場合には,以下のように書けば良いようです.

julia> rand(1:6, 20) # rand(start:stop, times)
20-element Array{Int64,1}:
 4
 3
 5
 5
 2
 3
 2
 2
 5
 3
 3
 1
 4
 4
 6
 5
 4
 4
 2
 5

上記のREPLの例にもコメントしてありますが,関数rand()の使い方は以下のようになります.

   rand(start:stop, times)


上記の例は復元抽出なので,非復元抽出を行う場合は以下のように書きました.

>>> np.random.choice(np.arange(1, 7), 6, replace = False)
array([5, 3, 4, 6, 1, 2])


これをJuliaで行う場合には,StatsBaseパッケージを用いて,以下のように書けば良いようです.

julia> using StatsBase
 
julia> sample(1:6, 6, replace=false)
6-element Array{Int64,1}:
 6
 4
 1
 3
 2
 5
関数sample()の使い方は以下のようになります.
    sample(start:stop, times, replace = "true" or "false")
 
なので,最初のコードと同じことを実行するためには,以下のように書いても良いようです.

julia> sample(1:6, 20, replace=true)

20-element Array{Int64,1}:

 4

 1

 6

 4

 1

 5

 5

 6

 4

 2

 6

 3

 2

 1

 1

 6

 2

 1

 5

 1

なお,上記の例のJupyter Notebookファイルは,GitHubPython_vs_Julia_Julia_03というファイルで見ることができます.

-Julia- パッケージのインストール

Juliaにおいてパッケージをインストールした際のメモです.

現時点ではバージョン1.1ですが,Julia 1.0以降では"pkgモード"が作られています.
REPLを起動して"]"を入力するとプロンプトが "(v1.1)pkg>"と変わります.
この状態で
    add "Package Name"
とするとインストールが始まります.以下は,"StatsBase"パッケージをインストールした際の例です.

(v1.1) pkg> add StatsBase
  Updating registry at `~/.julia/registries/General`
 Resolving package versions...
 Installed Missings ─────────── v0.4.0
 Installed SortingAlgorithms ── v0.3.1
 Installed DataStructures ───── v0.15.0
 Installed StatsBase ────────── v0.30.0
 Installed OrderedCollections ─ v1.1.0
  Updating `~/.julia/environments/v1.1/Project.toml`
  [2913bbd2] + StatsBase v0.30.0
  Updating `~/.julia/environments/v1.1/Manifest.toml`
  [864edb3b] + DataStructures v0.15.0
  [e1d29d7a] + Missings v0.4.0
  [bac558e1] + OrderedCollections v1.1.0
  [a2af1166] + SortingAlgorithms v0.3.1
  [2913bbd2] + StatsBase v0.30.0
(v1.1) pkg> 

pkgモードからは,"delete"を押すか,Ctr+Cで抜けることができます.