[R] Windows10にRStudioをインストール

WIndows10にRStudioをインストールした際のメモです．

まずは，CRANからRStudioのインストーラーをダウンロードします．

"Download R for Windows" を選択して，その後，base > Download R 3.6.0 for Windows （3.6.0は現時点でのバージョン）をクリックするとインストール実行ファイルのダウンロードが始まります．

その後，ダウンロードした"R-3.6.0-win.exe （3.6.0は現時点でのバージョン）"を実行するとインストールが始まります．
基本的にデフォルトのまま進めば良いのですが，Windows10は64bit OSなので，32-bit Filesのチェックは外します．

インストールが完了したら，R（Rx64 3.6.0）を起動します．

外部パッケージをインストールするには，あらかじめ，外部パッケージをダウンロードしておきます（今回はTSSSをインストールしたいので，こちら（TSSS_1.2.4.zip）からダウンロード）．

メニューの "パッケージ > Install package(s) from local files..." を選択し，"Select files"でダウンロードした TSSS_1.2.4.zip を選択します．

メニューのパッケージから "パッケージの読み込み"，"Select one"で"TSSS"を選択します．

これで，パッケージの読み込みは完了です．

試しに以下のように入力してみます．

うまくいっているようなので，とりあえずこれで良しとします．

[R] Jupyter NotebookでRを使う

以前に，Homebrew経由でRをインストールする記事を投稿しました．
今回は，RをJupyter Notebookで使えるようにするための手順のメモです．

まずは，ターミナルからRを起動します．

$ r

続いて，IRkernelのページにある手順に従って必要になるパッケージをインストールします．

> options(repos="https://cran.ism.ac.jp")

> install.packages(c('rzmq','repr','IRkernel','IRdisplay'),repos = c('http://irkernel.github.io/', getOption('repos')))

パッケージのインストールが完了したら，Rのコマンドラインから，以下を実行します．

> IRkernel::installspec()

その結果，以下のように表示されれば完了です．

[InstallKernelSpec] Installed kernelspec ir in /Users/hide/Library/Jupyter/kernels/ir

Jupyter Notebookを起動して，"New"を選択（クリック）した際に，以下のように"R"が表示されれば成功です．

ただ，現時点では外部のパッケージをJupyter Notebookで使う設定ができていません...

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

以下に，円周率のモンテカルロ計算を実行した時にかかる時間の計測例を示します．
（オリジナルは，黒木玄さんの「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.13768

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

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ファイルは，GitHubのMechanism of calculation.ipynbというファイルで見ることができます．

[macOS] HomeBrew list コマンド

HomeBrewでインストールしたソフトを確認するには，ターミナルで以下のように入力します．

$ brew list

すると，HomeBrew経由でインストールしたものの一覧が表示されます．

一覧の中でアンインストールしたいものがあるときは，以下のように入力します．

$ brew uninstall xxx

HomeBrew自体のアンインストールを行いたい場合は，以下のように入力します．

$ ruby -e “$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/uninstall)”

[Statistics] 標準偏差

Julia, Python, Rで標準偏差の計算を行う際のメモです．

標準偏差に関しては，下式

$\sigma = \sqrt{\sigma^2} = \sqrt{\cfrac{1}{N} ( x_i - \mu )^2}$

の標本分散の平方根とする場合と，下式

$\sigma^2 = \cfrac{1}{N-1} ( x_i - \mu )^2$

の不偏分散の平方根と説明してある場合があります．

名称に関してはどう呼ぶのが正しいのかはわかりませんが，以下の説明では最初の式を標本の標準偏差，2番目の式を不偏標準偏差と呼ぶことにします．

Julia, Python, Rの各言語の標準偏差を求める関数を用いたときに標本の標準偏差と不偏標準偏差のどちらが返ってくるのかを調べた結果です．結論としては，3種類の言語ともに（デフォルトでは）不偏標準偏差の値を返してきます．

Julia（JuliaのStatisticsパッケージのドキュメントはこちら）
Statistics パッケージの std 関数で標準偏差を求めると，不偏標準偏差が計算されます．これを標本の標準偏差だと思っていると当然，計算結果が合わないので悩むことになります．

以下に，その例を示します．
まずは，ターミナルからJuliaを立ち上げて以下のように1〜10の平均を求めます．

julia> using Statistics

julia> x = collect(1:1:10)         #[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

10-element Array{Int64,1}:

  1

  2

  3

  4

  5

  6

  7

  8

  9

 10

julia> μ = Statistics.mean(x)     # \mu[tab]

5.5

julia> n = length(x)

10

ここで，定義に従って分散を求めると以下の値になります．

julia> sum((x .- μ).^2) / n

8.25

不偏分散を求めると以下の値となります．

julia> sum((x .- μ).^2) / (n - 1)

9.166666666666666

標本の標準偏差を求めると以下の値になります．

julia> sqrt(sum((x .- μ).^2) / n)

2.8722813232690143

不偏標準偏差を求めると以下のようになります．

julia> sqrt(sum((x .- μ).^2) / (n - 1))

3.0276503540974917

続いて，JuliaのStatisticsモジュールのvar（分散）関数を使うと以下のようになります．

julia> s²₂ = Statistics.var(x)

9.166666666666666

この値は不偏分散です（この例では，Statistics.var としていますが，var()だけでも良い）．

std（標準偏差）関数を使うと以下のようになります．

julia> s₂ = Statistics.std(x)

3.0276503540974917

この値は不偏標準偏差の値です（この例では，Statistics.std としていますが，std()だけでも良い）．

以上から，JuliaのStatisticsモジュールのstd関数を用いると不偏標準偏差の値が返ってくることが確認できました．Jupyter Notebookでの実行例はGiHubで見ることができます．コードもGitHubで見ることができます．

Python（Numpyのstd関数に関するドキュメントはこちら）
Juliaと同様にPythonでも確認してみます．

>>> import numpy as np

標準偏差を求める配列を設定します．

>>> x = [i for i in range(1,11)]

デフォルトでのstd関数を用いた場合の値を求めます．

>>> s2 = np.std(x) # np.std(x, ddof=0)

>>> print ("s2 = ", s2)

s2 =  2.8722813232690143

上記のs2の計算式をしている行でもコメントアウトしていますが，デフォルトではddof = 0として計算を行っており，不偏標準偏差の値が返ってきています．

標本の標準偏差を求める場合には，以下のように引数の設定でddof = 1とします．

>>> s1 = np.std(x, ddof=1)

>>> print ("s1 = ", s1)

s1 =  3.0276503540974917

以上から，PythonのNumpyモジュールのstd関数を用いると不偏標準偏差の値が帰ってくることが確認できました．Jupyter Notebookでの実行例はGitHubで見ることができます．コードもGitHubで見ることができます．

R
Rでも確認してみます．
まずは，標準偏差を求める配列を設定します．

> x = (1:10)

標準偏差を求めてみます．

> s2 = sd(x)

> s2

[1] 3.02765

不偏標準偏差の値が返ってきています．

標本の標準偏差を求めるには以下のようにします．

> v = var(x)*(length(x)-1)/length(x)

> s1 = sqrt(v)

> s1

[1] 2.872281

以上から，Rのsd関数を用いると不偏標準偏差の値が帰ってくることが確認できました．Jupyter Notebookでの実行例はGitHubで見ることができます．コードもGitHubで見ることができます．

[Julia] バイオリンプロットを描く

前回はCSVファイルを読み込む際のメモを投稿しましたが，CSVファイルで読み込んだデータからバイオリンプロットを描いてみた際のメモです．

まずはCSVパッケージを読み込みます．

julia> using CSV

続いてCSVデータを読み込みます．読み込むデータは 3-3-2-fish_multi.csv*というcsvファイルです．

julia> fish_multi = CSV.read("3-3-2-fish_multi_2.csv", header = true, delim = ',')

20×2 DataFrames.DataFrame

│ Row │ species │ length │

│     │ String  │ Int64  │

├─────┼─────────┼────────┤

│ 1   │ A       │ 2      │

│ 2   │ A       │ 3      │

│ 3   │ A       │ 3      │

│ 4   │ A       │ 4      │

│ 5   │ A       │ 4      │

│ 6   │ A       │ 4      │

│ 7   │ A       │ 4      │

│ 8   │ A       │ 5      │

│ 9   │ A       │ 5      │

│ 10  │ A       │ 6      │

│ 11  │ B       │ 5      │

│ 12  │ B       │ 6      │

│ 13  │ B       │ 6      │

│ 14  │ B       │ 7      │

│ 15  │ B       │ 7      │

│ 16  │ B       │ 7      │

│ 17  │ B       │ 7      │

│ 18  │ B       │ 8      │

│ 19  │ B       │ 8      │

│ 20  │ B       │ 9      │

続いて，描画のためのライブラリを読み込みます．

julia> using PyCall, PyPlot, DataFrames

julia> sns = pyimport("seaborn")

julia> pd  = pyimport("pandas")

データをPandasのDataFrameに渡します．

julia> fish_multi_Pd = pd.DataFrame( 

           Dict(

               "species" => fish_multi.species, 

               "length"  => fish_multi.length 

               )

           );

Seabornを使って描画します．

julia> sns.violinplot(x = "species", y = "length",

                             data = fish_multi_Pd, color = "gray")

上記の事柄をJupyter Notebookで実行した結果はこちら（GiHub）で見ることができます．また，同じことをPythonで実行した結果はこちら（Github）で見ることができます．

* 3-3-2-fish_multi.csv の中身

species,length

A,2

A,3

A,3

A,4

A,4

A,4

A,4

A,5

A,5

A,6

B,5

B,6

B,6

B,7

B,7

B,7

B,7

B,8

B,8

B,9

[Julia] CSVデータを読み込む

JuliaでCSVデータを読み込む際のメモです．

まずはdownload関数でサンプルデータファイルをインターネットからダウンロードして，データフレームに読み込みます．以下の例では，インターネット上のCSVデータを読み込んでいます．

julia> download("https://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data", "iris.csv")

  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current

                                 Dload  Upload   Total   Spent    Left  Speed

100  4551  100  4551    0     0   3003      0  0:00:01  0:00:01 --:--:--  3003

"iris.csv"

isfile関数でファイルがダウンロードできたことを確認します．

julia> isfile("iris.csv")

true

readline関数で1行目を表示します．

julia> readline("iris.csv")

"5.1,3.5,1.4,0.2,Iris-setosa"

readline関数による表示で，1行目がカンマで区切られていて，ピリオドは小数点でヘッダ行がないことがわかります．

CSV.read関数でデータを読み込み，データフレームを確認する．CSV.read関数はCSV形式のファイルを読み込んで，内容をDataFrameに書き込む関数．

julia> using CSV, DataFrames

julia> df = CSV.read("iris.csv", 

                     header = ["PetalLength", "PetalWidth", "SepalLength", "SepalWidth", "Class"]);

thread = 1 warning: only found 1 / 5 columns on data row: 151. Filling remaining columns with `missing`

julia> describe(df)

5×8 DataFrame. Omitted printing of 1 columns

│ Row │ variable    │ mean    │ min         │ median │ max            │ nunique │ nmissing │

│     │ Symbol      │ Union…  │ Any         │ Union… │ Any            │ Union…  │ Int64    │

├─────┼─────────────┼─────────┼─────────────┼────────┼────────────────┼─────────┼──────────┤

│ 1   │ PetalLength │ 5.84333 │ 4.3         │ 5.8    │ 7.9            │         │ 1        │

│ 2   │ PetalWidth  │ 3.054   │ 2.0         │ 3.0    │ 4.4            │         │ 1        │

│ 3   │ SepalLength │ 3.75867 │ 1.0         │ 4.35   │ 6.9            │         │ 1        │

│ 4   │ SepalWidth  │ 1.19867 │ 0.1         │ 1.3    │ 2.5            │         │ 1        │

│ 5   │ Class       │         │ Iris-setosa │        │ Iris-virginica │ 3       │ 1        │

データフレームの最後の6行を見てみます．

julia> last(df, 6)

6×5 DataFrame

│ Row │ PetalLength │ PetalWidth │ SepalLength │ SepalWidth │ Class          │

│     │ Float64⍰    │ Float64⍰   │ Float64⍰    │ Float64⍰   │ String⍰        │

├─────┼─────────────┼────────────┼─────────────┼────────────┼────────────────┤

│ 1   │ 6.7         │ 3.0        │ 5.2         │ 2.3        │ Iris-virginica │

│ 2   │ 6.3         │ 2.5        │ 5.0         │ 1.9        │ Iris-virginica │

│ 3   │ 6.5         │ 3.0        │ 5.2         │ 2.0        │ Iris-virginica │

│ 4   │ 6.2         │ 3.4        │ 5.4         │ 2.3        │ Iris-virginica │

│ 5   │ 5.9         │ 3.0        │ 5.1         │ 1.8        │ Iris-virginica │

│ 6   │ missing     │ missing    │ missing     │ missing    │ missing        │

julia> eltypes(df)

┌ Warning: `eltypes(df::AbstractDataFrame)` is deprecated, use `eltype.(eachcol(df))` instead.

│   caller = top-level scope at none:0

└ @ Core none:0

5-element Array{Union,1}:

 Union{Missing, Float64}

 Union{Missing, Float64}

 Union{Missing, Float64}

 Union{Missing, Float64}

 Union{Missing, String} 

データ中の最後の行には，missingしかありません（ファイルの末尾に余計な改行が入っていたため）．この問題を回避するために，最後の行を切り捨てて，行の型をmissing値を許さない方に変換し，データフレーム中にmissing値がないようにします．

julia> df = disallowmissing!(df[1:end-1, :]);

julia> eltype(df)

Any

データを整理できたので，CSVファイルに書き出します．

julia> CSV.write("iris2.csv", df);

上記の事柄をJupyter Notebook上で実行してみた結果は，GitHubで見ることができます．