R グラフ描画

データの特徴を知るために，グラフ描画を行なってみた際のメモです．この内容は， 統計ソフト「R」超入門〈最新版〉  の記載を参考に処理してみた結果です． 

インデックスプロットは，各ケースのインデックスを横軸に，その観測値を縦軸に取ってプロットしたグラフです．インデックスプロットによって，データを全体的にみて，分布の特徴を大まかに捉えたり，外れ値を発見したり，ケースの並びに傾向があるかどうかを知ることができます．

Rでインデックスプロットを作成するには，R Commander のメニューからGraph > Index plot... を選択します．

すると，Index Plot ウィンドウが開きます．Data タブで Variable (pick one or more) の枠からプロットする変数を選択します（ここでは StressResponseScoreを選択）．

Option タブに切り替えて，Plot Options（プロットオプション） の Style of plot （プロットのスタイル）から"Spike（スパイク）"を選びます．必要に応じて Plot Labels（ラベルを表示）に y-axis label（y軸のラベル），Graph title（グラフのタイトル）を指定します．

完成したら，OKボタンをクリックするとグラフが描画されます．

グラフのタイプをpoints（点）にすると，以下のようになります．

ドットプロットは，横軸の変数の該当する値の位置にデータの個数分ドットを上に積み上げて分布をみるためのグラフです．因子の水準別にドットプロットを描くことで，水準間の分布の比較が可能になります．

R Commanderのメニューから Graph > Dot Plot を選ぶとドットプロットウィンドウが現れます．

Variable (pick one) （変数（1つ選別））の枠から，今回はLocusOfControl（健康統制感）を選んで Plot by groups...（層別のプロット）をクリックします．すると，以下のようなウィンドウが開くので，Sex（性別）を選択してOKボタンをクリックします．

このドットプロットでは，by = 1が男性の健康統制感のプロット図で，by=2が女性のプロット図です．この二つを比較すると，健康統制感には男女で大きな違いはないように感じます．

Rの付属データを読み込んでグラフを描いてみる

 今回は，Rのパッケージ"datasets"に付属するデータセット"airquality"を読み込みます．この内容は， 統計ソフト「R」超入門〈最新版〉  の記載を参考に処理してみた結果のメモです．

手順は以下のとおりです．

1. R Commander のメニューから Data > Data In Package > Read data set from an attached package... を選択すると，読み込むパッケージとデータセットを選択するための"Read Data From Package"ウィンドウが表示される．
2. Package (Double-click to select) の中からパッケージを選択してダブルクリックする（今回は"datasets"）．Data set (Double-click to select) の枠内に"datasets"パッケージに含まれるデータセットの一覧が表示される．
3. 読み込みたいデータセット（今回は"airquality"）をダブルクリックして選択する．この時，"Help on selected data set"をクリックするとデータセット（今回は"airquality"）の説明が表示される．ちなみに，このデータはニューヨークの1973年5月1日〜9月30日における大気（オゾン濃度，日射量，気温，風速，月，日）を測定したデータ．
4. "OK"ボタンをクリックすると，データセットが読み込まれます．

ここまでの操作で，R Script には以下のような表示がされています．

    data()

    help("airquality", package="datasets")

    data(airquality, package="datasets")

Output には以下のような表示がされています．

    > data()

    > help("airquality", package="datasets")

    > data(airquality, package="datasets")

また，R Scriptのメニュー2行目にある" Dataset : "の右にあるボタンには読み込まれているデータセットとして"airquality"と表示されているはずです．

R Commander の”View data set"ボタンをクリックすると，読み込んだ"airquality"のデータを表示でき，"airquality"のデータの変数は6つであることがわかります．

ここまでが，データセットの読み込みです．続いてグラフ（散布図）を描いてみます．

R Commander のメニューから Graph > Scatterplot を選択すると，Scatterplot... ウィンドウが表示されます．

右下の OK ボタンをクリックすると，散布図が表示されます．

----------

上記の処理は，MacBookPro (2017, Intel) で行っていますが，プルダウンメニューの反応に時間がかかることがあります． 

macOSにRのRコマンダーをインストールする．

CRANからRのインストーラーをダウンロードします．

macOSへのインストールなので，Download R for macOS を選択します．

Intel Macの場合はFor older Intel Macs: （2023/08/14時点で R-4.3.1-x86_64.pkg），M1/M2 Macの場合は For Apple silicon (M1/M2) Macs: （2023/08/14時点で R-4.3.1-arm64.pkg）をクリックして，インストーラーをダウンロードします．

ダウンロードしたインストーラーを起動して，Rをインストールします．続いて，XQuartzをインストールします．

Rを立ち上げると，R Consoleが立ち上がります．試しに，以下のように平均を求めるコマンド"mean()"と数値をベクトル形式で入力するコマンド"c()"を使用してみます．

    > mean(c(4.8, 5.2, 6, 8))

すると，"6"という答えが返ってきました．これで，Rが上手く動いていることが確認できました．

R Commander のインストールには直接関係ありませんが，

    Package & Data > Package Installer

を立ち上げて，Update All をクリックしてパッケージを更新しておきます．

続けて，R Consoleで以下のコマンドを実行します．

    > install.packages("Rcmdr")

    > library(Rcmdr)

すると，下図のように R Commander が立ち上がります（理由はよくわかりませんが，インストール直後にR Commanderが立ち上がらないことがあります．このような時は，Macを再起動して，再度Rの立ち上げからを繰り返すと上手くいくことがあります）．

次からは，Rを立ち上げた後に

    > library(Rcmdr)

を実行するとR commanderが立ち上がるようになりました．

R Commanderを使ってみます．これ以降は 統計ソフト「R」超入門〈最新版〉  の記載を参考に進めてみた結果のメモです．

R Commander のR Script欄に

    4/3*pi*4.8^3

と入力（半径4.8の球の体積を求める計算式）してSubmitをクリックするとOutput欄に以下のように返ってきます．

    > 4/3*pi*4.8^3

    [1] 463.2467

続いて，異なる半径 4.8, 5.2, 6, 8 の4つの球の体積を同時に計算してみます．R Scriptに

    r <- c(4.8, 5.2, 6, 8)

と入力します．ここで，<- は左辺の変数（r）に右辺で記述された値を代入する記号です．

そして r と入力して Submit をクリックすると，Outputnに以下のように返ってきます．

    > r

    [1] 4.8 5.2 6.0 8.0

次に半径4.8の球の体積を求めた式の半径4.8を r に置き換えて，半径 r の球の体積を求めた結果を変数 v に代入する式に修正してみます．

    4/3*pi*r^3 -> v

ここで， -> は左辺の計算を右辺の変数（v）に代入する記号です．

そして，R Script にv と入力してSubmit をクリックすると，Outputに以下のように返ってきます．

    > v

    [1]  463.2467  588.9774  904.7787 2144.6606

上記の4つの体積の平均をとるには，関数 mean() を用います．R Script欄に mean(v) と入力してSubmit をクリックすると，Output欄に以下のように返ってきます．

    > mean(v)

    [1] 1025.416

Pythonのモジュールを追加する

Pythonのモジュールを追加する際のメモです．

Pythonでモジュールがインストールされていないというエラーが出た際には，インストールする必要があります．

例えば，Jupyter-Notebookで以下のようなエラーが出たとします．

ModuleNotFoundError: No module named 'requests'

そんな時には，ターミナルから，以下のコマンドで requests というモジュールをインストールします．

% python3 -m pip install --upgrade requests

すると，以下のように流れてインストールが完了します．

Collecting requests

  Downloading requests-2.28.2-py3-none-any.whl (62 kB)

     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 62.8/62.8 kB 2.7 MB/s eta 0:00:00

Collecting charset-normalizer<4,>=2

  Downloading charset_normalizer-3.0.1-cp310-cp310-macosx_11_0_arm64.whl (122 kB)

     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 122.5/122.5 kB 4.6 MB/s eta 0:00:00

Collecting certifi>=2017.4.17

  Downloading certifi-2022.12.7-py3-none-any.whl (155 kB)

     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 155.3/155.3 kB 3.1 MB/s eta 0:00:00

Collecting idna<4,>=2.5

  Downloading idna-3.4-py3-none-any.whl (61 kB)

     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 61.5/61.5 kB 4.8 MB/s eta 0:00:00

Collecting urllib3<1.27,>=1.21.1

  Downloading urllib3-1.26.14-py2.py3-none-any.whl (140 kB)

     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 140.6/140.6 kB 5.0 MB/s eta 0:00:00

Installing collected packages: charset-normalizer, urllib3, idna, certifi, requests

Successfully installed certifi-2022.12.7 charset-normalizer-3.0.1 idna-3.4 requests-2.28.2 urllib3-1.26.14

Juliaで作って学ぶベイズ統計学のコードを実行してみました

Juliaで作って学ぶベイズ統計学のJulia コードを実行してみた結果を，GItHubで公開してみています．

1.  Juliaの基礎
2. 数値計算の基礎
3. 確率計算の基礎
4. 確率分布の基礎
5. 統計モデリングと推論
6. 勾配を利用した近似推論手法
7. 発展的な統計モデル

Jupyter Notebookが完成し次第，随時アップロードしていく予定です．

JuliaをアップデートしてJupiter Notebookでも使えるようにする

 アップデートなので，すでにJuliaがインストールされている前提です．

Homebrewを使います．まずはターミナルから，以下のようしてJuliaをアップデートします．

% brew install julia --cask

以前は，以下のように cask install juliaだったのですが，エラーが出てしまいました．

% brew cask install julia

Error: `brew cask` is no longer a `brew` command. Use `brew <command> --cask` instead.

% brew cask install julia

アップデートが終わったら，Applicationフォルダに新しいJulia（以下の例では，最初に1.0がインストールされていて，新たに1.8をインストールした）のアイコンがあるので，アイコンをダブルクリックして立ち上げます．

その上で，"]"キーを押してパッケージモードに切り替えて，以下のように（add IJulia）IJuliaパッケージを追加します．

/Applications/Julia-1.8.app/Contents/Resources/julia/bin/julia ; exit;

% /Applications/Julia-1.8.app/Contents/Resources/julia/bin/julia ; exit;

               _

   _       _ _(_)_     |  Documentation: https://docs.julialang.org

  (_)     | (_) (_)    |

   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.

  | | | | | | |/ _` |  |

  | | |_| | | | (_| |  |  Version 1.8.1 (2022-09-06)

 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release

|__/                   |

(@v1.8) pkg> add IJulia

この↑例では省略していますが，add IJuliaを実行すると，ずらずらとインストールが始まります．

IJuliaパッケージの追加が終わった段階で，Available kernelsを確認すると，julia-1.8 が現れます．

% jupyter kernelspec list

Available kernels:

  julia-1.0    /Users/hide/Library/Jupyter/kernels/julia-1.0

  julia-1.8    /Users/hide/Library/Jupyter/kernels/julia-1.8

  python3      /usr/local/share/jupyter/kernels/python3

Jupyter Notebookを立ち上げて，インストールしたJuliaを使えるかどうか確認しておきます．

Julia（今回は1.8.1）が使えるようになっていれば成功です．

使わないカーネルや，存在しないカーネル（今回の例では，Julia 1.0.5）を削除するにはターミナルから以下のように入力します．なお，"KERNEL_NAME"には使わなくなった，もしくは存在しないカーネルの名前を書きます．

% jupyter kernelspec uninstall KERNEL_NAME

今回の例で，julia-1.0を削除すると以下のようになります．

% jupyter kernelspec uninstall julia-1.0  

Kernel specs to remove:

  julia-1.0           /Users/xxx/Library/Jupyter/kernels/julia-1.0

Remove 1 kernel specs [y/N]: y

[RemoveKernelSpec] Removed /Users/xxx/Library/Jupyter/kernels/julia-1.0

今回の例では，Juliaを1.0から1.8にアップデートしました．1.0でインストールしたパッケージは継承されないようなので，改めて追加する必要があります．以下の例は，とりあえず追加しておいたパッケージ達です．

(@v1.8) pkg> add CSV

(@v1.8) pkg> add Combinatorics

(@v1.8) pkg> add DataFrames

(@v1.8) pkg> add Distributions

(@v1.8) pkg> add JuMP

(@v1.8) pkg> add Plots

(@v1.8) pkg> add PyCall

(@v1.8) pkg> add PyPlot

(@v1.8) pkg> add RDatasets

(@v1.8) pkg> add StatsBase

(@v1.8) pkg> add StatsPlots

(@v1.8) pkg> add SumOfSquares

(@v1.8) pkg> add Statistics

(@v1.8) pkg> add Gadfly

(@v1.8) pkg> add LinearAlgebra

(@v1.8) pkg> add Optim

インストールしたパッケージを確認するには，以下のように"status"コマンドを実行します．

(@v1.8) pkg> status

Status `~/.julia/environments/v1.8/Project.toml`

  [336ed68f] CSV v0.10.4

  [861a8166] Combinatorics v1.0.2

  [a93c6f00] DataFrames v1.3.5

⌅ [31c24e10] Distributions v0.23.4

  [c91e804a] Gadfly v1.3.4

  [7073ff75] IJulia v1.23.3

  [4076af6c] JuMP v1.3.0

  [5424a776] Mamba v0.12.5

  [429524aa] Optim v1.7.2

  [91a5bcdd] Plots v1.32.1

  [438e738f] PyCall v1.94.1

  [d330b81b] PyPlot v2.11.0

  [ce6b1742] RDatasets v0.7.7

  [2913bbd2] StatsBase v0.33.21

⌃ [f3b207a7] StatsPlots v0.15.1

  [4b9e565b] SumOfSquares v0.6.2

  [37e2e46d] LinearAlgebra

  [10745b16] Statistics

Info Packages marked with ⌃ and ⌅ have new versions available, but those with ⌅ cannot be upgraded. To see why use `status --outdated`

ベータ分布

   Pythonによるベイズ統計学入門に掲載されているコードを実行してみた際のメモです．

ベータ分布を描画するPythonのコード．

# -*- coding: utf-8 -*-

#%% NumPyの読み込み

import numpy as np

# SciPyのstatsモジュールの読み込み

import scipy.stats as st

# MatplotlibのPyplotモジュールの読み込み

import matplotlib.pyplot as plt

#%% ベータ分布の確率密度関数

q = np.linspace(0, 1, 250)

value_a = np.array([0.5, 1.0, 2.0, 4.0])

value_b = np.array([0.5, 1.0, 2.0, 4.0])

rows = value_a.shape[0]

cols = value_b.shape[0]

fig, ax = plt.subplots(rows, cols, sharex='all', sharey='all',

num=1, facecolor='w')

ax[0, 0].set_xlim(0.0, 1.0)

ax[0, 0].set_ylim(0.0, 4.5)

for row_index in range(rows):

a = value_a[row_index]

ax[row_index, 0].set_ylabel('$\\alpha$ = {0:3.1f}'.format(a),

fontsize=12)

for column_index in range(cols):

b = value_b[column_index]

ax[row_index, column_index].plot(q, st.beta.pdf(q, a, b), 'k-')

if row_index == 0:

ax[0, column_index].set_title('$\\beta$ = {0:3.1f}'.format(b),

fontsize=12)

plt.tight_layout()

plt.savefig('pybayes_fig_beta_distribution.png', dpi=300)

plt.show()

このコードを実行すると，以下の図が描画されます．この図はα，βを変えた時のベータ分布の形状です．