LinuxでRAID構築

前にラズパイで構築したサーバーを使って来たけどHDD故障でデータが飛ぶのが怖くなってきた。。。
なのでもう一枚HDDを用意してRAID1を構築した。

今回はその手順について。

目次

はじめに

RAIDのイメージ

Linux上のRAIDのイメージは以下の図のようになる。
複数のHDDをまとめて1つのデバイス(Multiple Device)として扱う。

f:id:ricrowl:20201130045630p:plain
RAID1のイメージ図

HDD確認

まずはRAIDに使いたいHDD確認(今回は実験のためUSBメモリ使用)。
以下で確認できる。

$ sudo fdisk -l
...
Disk /dev/sda: 14.5 GiB, 15504900096 bytes, 30283008 sectors
...
Disk /dev/sdb: 14.9 GiB, 15938355200 bytes, 31129600 sectors
...

今回はsdaとsdbのHDDでRAIDを構築する。

RAIDパーティション作成

新規パーティション作成

sdaとsdbのHDDにRAIDパーティションを作成する。
※このときHDDのデータは消去されるのでどこかにコピーしておく!

$ sudo fdisk /dev/sda

Welcome to fdisk (util-linux 2.33.1).
Changes will remain in memory only, until you decide to write them.
Be careful before using the write command.


Command (m for help): n
Partition type
   p   primary (0 primary, 0 extended, 4 free)
   e   extended (container for logical partitions)
Select (default p): p
Partition number (1-4, default 1): 
First sector (2048-30283007, default 2048): 
Last sector, +/-sectors or +/-size{K,M,G,T,P} (2048-30283007, default 30283007): 

Created a new partition 1 of type 'Linux' and of size 14.4 GiB.

特にこだわりがなければデフォルトのままEnterを押していけば良い。

パーティションが作成できないとき

以下のようにパーティション作成ができない場合は一回パーティションを消してから実行する。

$ sudo fdisk /dev/sda

Welcome to fdisk (util-linux 2.33.1).
Changes will remain in memory only, until you decide to write them.
Be careful before using the write command.


Command (m for help):  n
All space for primary partitions is in use. # パーティション作成できない

Command (m for help): d # パーティション削除
Selected partition 1
Partition 1 has been deleted.

パーティションタイプの変更

パーティションタイプをRAID用のLinux raid autoに変更する。

$ sudo fdisk /dev/sda

Welcome to fdisk (util-linux 2.33.1).
Changes will remain in memory only, until you decide to write them.
Be careful before using the write command.


Command (m for help): t
Selected partition 1
Hex code (type L to list all codes): fd
Changed type of partition 'Linux' to 'Linux raid autodetect'.

HDDの変更を保存

下記で変更を保存する。

Command (m for help): w
The partition table has been altered.
Calling ioctl() to re-read partition table.
Syncing disks.

一連の作業をsdbでも同様にして行う。

RAID設定

まずはRAID用パッケージmdadmをインストール。

$ sudo apt update
$ sudo apt install mdadm

mdadmによりRAIDを構築する。

$ sudo mdadm --create /dev/md0 --raid-devices=2 --level=raid1 /dev/sda1 /dev/sdb1
  • 基本構成: mdadm [モード] <RAIDデバイス名> [オプション] <RAID構成デバイス>
  • --create: RAID作成モード
  • /dev/md0: RAIDバイス名を/dev/md0に
  • --raid-devices=2: RAID構成デバイス数2
  • --level=raid1: RAID種類をRAID1に
  • /dev/sda1 /dev/sdb1: RAID構成デバイスを/dev/sda1, /dev/sdb1に

RAID構成デバイスの容量に差が1%以上あると以下の警告がでる。
特に問題がなければyで続行。
RAIDバイスの容量は小さい方のデバイスの容量となる。

mdadm: largest drive (/dev/sdb1) exceeds size (15131264K) by more than 1%
Continue creating array? y
mdadm: Defaulting to version 1.2 metadata
mdadm: array /dev/md0 started.

上記のコマンドを実行するとしばらくRAID構成の処理がバックグラウンドで実行される(結構時間かかる)。
進捗状況は以下で確認できる。

$ sudo watch cat /proc/mdstat
Personalities : [raid1] 
md0 : active raid1 sdb1[1] sda1[0]
      15131264 blocks super 1.2 [2/2] [UU]
      [===========>.........]  resync = 58.8% (8901056/15131264) finish=12.0min speed=8620K/sec
      
unused devices: <none>

処理完了後、mkfsRAIDバイスをフォーマットする。

$ sudo mkfs -t ext4 /dev/md0

RAIDバイスの使用

構築したRAIDバイスは通常のHDDと同様にマウントして使用することができる。

$ mkdir /home/user00/Share/md0 #マウント先ディレクトリ作成
$ sudo mount /dev/md0 /home/user00/Share/md0 #マウント

マウント後、rootでないと書き込みが出来なくなっていて煩わしかったら、chownで所有者を変更する。

$ sudo chown -R user00:user00 /home/user00/Share/md0

設定の保存

再起動後も使えるように設定を保存する。
設定情報は以下で取得できる。

$ sudo mdadm --detail --scan
ARRAY /dev/md/0 metadata=1.2 name=raspberrypi:0 UUID=a6a07a86:3ea1907b:727d6c92:0bfc2e73

これを/etc/mdadm.confに書き込む。

$ sudo vim.tiny /etc/mdadm.conf
ARRAY /dev/md0 metadata=1.2 name=raspberrypi:0 UUID=a6a07a86:3ea1907b:727d6c92:0bfc2e73
~                                                                                                      
~ 
~

さらに起動時にマウントするよう/etc/fstabにも以下を追記。

/dev/md0      /home/user00/Share/md0     ext4    defaults,nofail        0       0
  • nofail: デバイスがマウントできなくても起動できるようにする(通常はエマージェンシーモードになってsshができない)

有事の際

構成HDDの故障検知

RAIDなので構成HDDに異常があってもmd0上のデータはいつもどおりアクセスできる。
しかし、気づかないうちに残りのHDDが故障してしまったらせっかくRAIDにした意味がなくなるので知る必要がある。

構成HDDの異常は以下で調べることができる。

$ sudo cat /proc/mdstat

実際に一方のHDDを抜いて挙動を見てみる。
sdbのHDDを抜いてみた。

# 正常時の出力
Personalities : [raid1] 
md0 : active raid1 sdb1[1] sda1[0]
      15131264 blocks super 1.2 [2/2] [UU]
      
unused devices: <none>

# 異常時の出力
Personalities : [raid1] 
md0 : active raid1 sda1[0]
      15131264 blocks super 1.2 [2/1] [U_]
      
unused devices: <none>

下記で更に詳細に調べることができる。

$ sudo mdadm --detail /dev/md0
/dev/md0:
           Version : 1.2
     Creation Time : Sun Nov 29 17:41:25 2020
        Raid Level : raid1
        Array Size : 15131264 (14.43 GiB 15.49 GB)
     Used Dev Size : 15131264 (14.43 GiB 15.49 GB)
      Raid Devices : 2
     Total Devices : 1
       Persistence : Superblock is persistent

       Update Time : Sun Nov 29 19:33:34 2020
             State : clean, degraded 
    Active Devices : 1
   Working Devices : 1
    Failed Devices : 0
     Spare Devices : 0

Consistency Policy : resync

              Name : raspberrypi:0  (local to host raspberrypi)
              UUID : a6a07a86:3ea1907b:727d6c92:0bfc2e73
            Events : 28

    Number   Major   Minor   RaidDevice State
       0       8        1        0      active sync   /dev/sda1
       -       0        0        1      removed

故障HDDの交換

故障したHDDを新しいのに交換する手順。

まず、故障したHDDをRAID構成から除外する。

$ sudo mdadm /dev/md0 --remove /dev/sdb1

次に新しいHDDを取り付け、RAIDパーティション作成(前述)を行う。

そして、新しいHDDをRAID構成に加える。

$ sudo mdadm --add /dev/md0 /dev/sdb1

はじめてRAIDを構成したときのように、しばらくRAID構成の処理がバックグラウンドで実行されるので待つ。

処理が完了したら交換完了。以下のように正常にRAIDが動いていることが確認できる。

$ sudo cat /proc/mdstat
Personalities : [raid1] 
md0 : active raid1 sdb1[2] sda1[0]
      15131264 blocks super 1.2 [2/2] [UU]
      
unused devices: <none>

その他

RAIDの解除/再構築

RAIDの解除。

$ sudo umount /dev/md0

$ sudo mdadm --stop /dev/md0
mdadm: stopped /dev/md0

RAIDの再構築。

$ sudo mdadm --assemble /dev/md0 /dev/sda1 /dev/sdb1
mdadm: /dev/md0 has been started with 2 drives.

$ sudo mount /dev/md0 /home/user00/Share/md0

XGboostとりあえず動かす

「とりあえずこれやっとけ」で有名なXGboost。とりあえず動かす。

使用するモジュールはこちら

xgboost.readthedocs.io

目次

回帰編

import

回帰では以下のモジュールを使用。

import xgboost as xgb
from sklearn import datasets
import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
  • datasets, pandas: サンプルデータセット
  • train_test_split: 訓練/テストデータ分割用
  • mean_squared_error: 平均二乗誤差計算用

データセット取得

今回はBoston house pricesを使用。

boston = datasets.load_boston()
X = pd.DataFrame(boston.data, columns=boston.feature_names)
y = boston.target

こんなデータが得られる。 f:id:ricrowl:20201020022107p:plainf:id:ricrowl:20201020021800p:plain

訓練/テストデータに分割

X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.3, shuffle=True)
  • test_size: テストデータの比率
  • shuffle: データをシャッフルして分割するかどうか

データをxgboost用のフォーマットに変換

dtrain = xgb.DMatrix(X_tr, label=y_tr)
dtest = xgb.DMatrix(X_te, label=y_te)

xgboostのパラメータ設定

params = {'objective':'reg:squarederror'}
num_round = 50
  • objective: 目的関数
  • num_round: 学習回数

普通の回帰なので目的関数に二乗誤差(reg:squarederror)を指定。
タスクによっては二乗対数誤差(reg:squaredlogerror)やロジスティック回帰誤差(reg:logistic)などを選択。

他の目的関数: https://xgboost.readthedocs.io/en/latest/parameter.html?highlight=multi%20softmax#learning-task-parameters

さらに学習係数や木の深さなどのパラメータを設定したいときにはparamsに加える。

# 学習係数:0.3, 木の深さ:6
params = {'objective':'reg:squarederror', 'eta':0.3, 'max_depth':6}

他のパラメータとか: https://xgboost.readthedocs.io/en/latest/parameter.html?highlight=multi%20softmax#parameters-for-tree-booster

学習

watchlist = [(dtrain, 'train'), (dtest, 'eval')]
bst = xgb.train(params, dtrain, num_round, evals=watchlist)

テストデータ評価

pred = bst.predict(dtest)
score = mean_squared_error(y_te, pred)  
print('score:', format(score))

重要度変数可視化

xgb.plot_importance(bst)

分類編

import

分類では以下のモジュールを使用。

import xgboost as xgb
from sklearn import datasets
import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
  • datasets, pandas: サンプルデータセット
  • train_test_split: 訓練/テストデータ分割用
  • accuracy_score: クラス分類の正答率計算用

データセット取得

今回はIrisを使用。

iris = datasets.load_iris()
X = pd.DataFrame(iris.data, columns=iris.feature_names)
y = iris.target

こんなデータが得られる。 f:id:ricrowl:20201020024452p:plain

訓練/テストデータに分割

X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.3, shuffle=True)
  • test_size: テストデータの比率
  • shuffle: データをシャッフルして分割するかどうか

データをxgboost用のフォーマットに変換

dtrain = xgb.DMatrix(X_tr, label=y_tr)
dtest = xgb.DMatrix(X_te, label=y_te)

xgboostのパラメータ設定

params = {'objective':'multi:softmax', 'num_class':3}
num_round = 50
  • objective: 目的関数
  • num_class: 分類のクラス数
  • num_round: 学習回数

普通の分類なので目的関数に交差エントロピー(multi:softmax)を指定。
さらに分類のクラス数も指定する必要がある。
二値分類の場合はロジスティック回帰誤差(binary:logistic)を選択。

他の目的関数: https://xgboost.readthedocs.io/en/latest/parameter.html?highlight=multi%20softmax#learning-task-parameters

さらに学習係数や木の深さなどのパラメータを設定したいときにはparamsに加える。

# 学習係数:0.3, 木の深さ:6
params = {'objective':'multi:softmax', 'num_class':3, 'eta':0.3, 'max_depth':6}

他のパラメータとか: https://xgboost.readthedocs.io/en/latest/parameter.html?highlight=multi%20softmax#parameters-for-tree-booster

学習

watchlist = [(dtrain, 'train'), (dtest, 'eval')]
bst = xgb.train(params, dtrain, num_round, evals=watchlist)

テストデータ評価

pred = bst.predict(dtest)
score = accuracy_score(y_te, pred)  
print('score:{0:.4f}'.format(score))

重要度変数可視化

xgb.plot_importance(bst)

matplotlibの日本語化(フォント変更)

matplotlibで日本語を表示しようとすると以下のように文字化けしてしまう。

import matplotlib.pyplot as plt
plt.plot(list(range(10)), list(range(10)))
plt.title('日本語タイトル')

f:id:ricrowl:20200914013415p:plain

これはmatplotlibで使っているフォントが日本語に対応していないからなので対応しているフォントを使うようにしてやる。

ちなみに今使っているフォントはplt.rcParams['font.family']を出力して確認できる。

フォントの取得

日本語対応したフォントを外部サイトから取得する。
情報処理推進機構がフォントを公開してくれている(DLリンク)。

ライセンスも商用・非商用が可なので使用しやすい。以下IPAexフォント Ver.004.01のライセンスから抜粋。

6.受領者は、3条2項の定めに従い、商用・非商用を問わず、許諾プログラムをそのままの状態で改変することなく複製して第三者への譲渡し、公衆送信し、その他の方法で再配布することができます(以下、「再配布」といいます。)。

ダンロードしたzip内のipaexg.ttfを使用する。

フォントをmatplotlibに追加

matplotlibのパッケージデータにipaexg.ttfを加える。

まず、plt.rcParams['datapath']を出力してパッケージデータが保存されているディレクトリを確認する。
おそらく.../matplotlib/mpl-dataのような構造となっている。

次に.../matplotlib/mpl-data/fonts/ttfにipaexg.ttfを追加する。自分の環境では以下のようにした。

cp ipaexg.ttf /opt/conda/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/

フォントをデフォルトとして登録

追加したフォントをmatplotlibのデフォルトとして登録する。

.../matplotlib/mpl-data/matplotlibrcの以下の部分を編集する。
編集前

# font.family         : sans-serif

編集後

font.family         : IPAexGothic

キャッシュを削除

matplotlibのフォントのキャッシュデータがあると追加したフォントが認識されないので削除する。

まず、以下のようにしてキャッシュデータのディレクトリを確認する。

import matplotlib 
print(matplotlib.get_cachedir())

おそらく.../.cache/matplotlibのような構造となっている。

次に.../.cache/matplotlib/下にfontがつくファイルがある(例えばfontlist-v300.jsonみたいな)ので全て削除する。自分の環境では以下のようにした。

rm /root/.cache/matplotlib/font*

プロット再確認

pythonやjupyterを再起動して冒頭と同じものをプロットしてみる。

import matplotlib.pyplot as plt
plt.plot(list(range(10)), list(range(10)))
plt.title('日本語タイトル')

f:id:ricrowl:20200914024846p:plain

日本語で表示できるようになった。

参考:一時的にフォントを変えたいだけのとき

matplotlibのデフォルトのフォントはそのままで一時的にだけ変えたい場合は以下2通りの方法がある。

特定のスクリプトだけフォントを変えたいとき

このスクリプトだけ日本語化したいなーというときはこっち。

スクリプトのはじめにplt.rcParams['font.family']を変更する。

plt.rcParams['font.family'] = 'IPAexGothic'

以下、ずっと変更したフォントで表示される。

特定の図だけフォントを変えたいとき

この図だけ日本語化したいなーというときはこっち。

変更したいフォントのFontPropertiesを作成し、pltの引数に加える。

import matplotlib.pyplot as plt
import matplotlib.font_manager as fm
font_path = '/opt/conda/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/ipaexg.ttf'
fp = fm.FontProperties(fname=font_path)
plt.plot(list(range(10)), list(range(10)))
plt.title('日本語タイトル', fontproperties=fp)
plt.text(5, 5, 'デフォルトの日本語タイトル')

f:id:ricrowl:20200914030101p:plain

plt.text(5, 5, 'デフォルトの日本語タイトル')はデフォルトのフォントのままで文字化けしているけど、plt.title('日本語タイトル', fontproperties=fp)は日本語化できている。

さらにこの方法ではフォントのディレクトリはどこでもOK。

参考:フォント名を確認する方法

今回ipaexg.ttfはIPAexGothicというフォント名だったけどファイル名からは名前がわからない。

そんなときはmatplotlibのフォントのキャッシュデータを見ると名前が書いてある。
自分の環境では/root/.cache/matplotlib/fontlist-v300.jsonを確認すると以下の部分が確認がある。

    {
      "fname": "fonts/ttf/ipaexg.ttf",
      "name": "IPAexGothic",
      "style": "normal",
      "variant": "normal",
      "weight": 400,
      "stretch": "normal",
      "size": "scalable",
      "__class__": "FontEntry"
    },

ここの"name"からIPAexGothicとわかる。

※キャッシュデータはフォント追加後に削除したあとのもの

PythonでMCMCを実装してみた

MCMCの理論 - ricrowlのブログの続き。 マルコフ連鎖モンテカルロ法(MCMC)のアルゴリズムであるMetropolis-Hastings(MH)法を実装してみた。

MH法アルゴリズム

確率変数$\boldsymbol{w}$のデータ$\boldsymbol{X}$観測後の事後分布$P(\boldsymbol{w}|\boldsymbol{X})$を求めたいとする。

MH法では、以下の遷移確率で定義されたマルコフモデルからサンプリングする。

$P(\boldsymbol{w}|\boldsymbol{X})=\frac{P(\boldsymbol{X}|\boldsymbol{w})P(\boldsymbol{w})}{\int P(\boldsymbol{X}, \boldsymbol{w})d\boldsymbol{w}}=\frac{f(\boldsymbol{w})}{Z}$としたとき、 $$ \begin{align} T(\boldsymbol{w}^{t+1}|\boldsymbol{w}^{t}) &= q(\boldsymbol{w}^{t+1}|\boldsymbol{w}^{t}) \times \min(1, \frac{ f(\boldsymbol{w^{t+1}}) q(\boldsymbol{w}^{t}|\boldsymbol{w}^{t+1}) }{ f(\boldsymbol{w^{t}}) q(\boldsymbol{w}^{t+1}|\boldsymbol{w}^{t}) }) \\ T(\boldsymbol{w}^{t}|\boldsymbol{w}^{t}) &= 1 - T(\boldsymbol{w}^{t+1}|\boldsymbol{w}^{t}) \end{align} $$ ここで、$q(\boldsymbol{w}'|\boldsymbol{w})$は提案分布と呼ばれ、$\boldsymbol{w}$が遷移する候補の値を取ってくるための分布。

サンプリング手順

  1. 提案分布$q(\boldsymbol{w}'|\boldsymbol{w})$から遷移先の候補点$\boldsymbol{w}'$を取得する
  2. 受理確率$a = \frac{ f(\boldsymbol{w^{t+1}}) q(\boldsymbol{w}^{t}|\boldsymbol{w}^{t+1}) }{ f(\boldsymbol{w^{t}}) q(\boldsymbol{w}^{t+1}|\boldsymbol{w}^{t}) }$を計算する
  3. $a$が1以上だったら候補点$\boldsymbol{w}'$をサンプルとして受理し、1より小さかったら$a$の確率で$\boldsymbol{w}'$をサンプルとして受理する
  4. 受理された場合は$\boldsymbol{w}\rightarrow\boldsymbol{w}'$と遷移させ、されなかった場合は遷移させず$\boldsymbol{w}\rightarrow\boldsymbol{w}$とする
  5. 上記を繰り返す

実装

上のアルゴリズムPythonで書いてみた。

import numpy as np

def MH_method(w, X, q, sampling_q, f):
    '''
    MH法によるサンプリング関数
    
    Args:
        w: サンプリング対象のパラメータの現在値
        X: データ
        q: 提案分布関数
        sampling_q: 提案分布から候補点を取得する関数
        f: wとXとの同時確率を計算する関数
        
    Returns:
        w_new: 新たに遷移したパラメータ
        is_accept: 受諾されたかどうか
    '''
    
    # 1. 提案分布から候補点を取得
    w_new = sampling_q(w)
    # 2. 候補点の受理確率計算
    a = (f(w_new, X)*q(w, w_new))/(f(w, X)*q(w_new, w))
    # 3. 受理判定
    is_accept = np.random.uniform() < min(1, a)
    return w_new, is_accept
    
    
def sampling_MH(X, q, sampling_q, f, sample_num, warmup, w_init):
    '''
    MH法でサンプルリストを取得する関数
    
    Args:
        X, q, sampling_q, f: MH_method参照
        sample_num: サンプル数
        warmup: はじめの方で定常化していないサンプルを無視する数
        w_init: wの初期値

    Returns:
        samples: サンプルリスト
    '''
    
    w = w_init
    samples = []
    while(len(samples) < sample_num + warmup):
        w_new, is_accept = MH_method(w, X, q, sampling_q, f)
        # 4. 受理された場合はwを遷移
        if is_accept:
            w = w_new
            samples.append(w)
    samples = samples[warmup:]
    return samples

関数q, sampling_q, f は予め設定する必要がある。

qは提案分布$q(\boldsymbol{w}'|\boldsymbol{w})$の確率を計算する関数であり、PRMLによると平均$\boldsymbol{w}$と適当な標準偏差\sigma_qガウス分布N({\bf{w}}'|{\bf{w}},\sigma^{2}_{q})とすれば良い。 sampling_qもこのガウス分布から取得する関数にすれば良い。

標準偏差\sigma_qの目安としては事後分布$P(\boldsymbol{w}|\boldsymbol{X})$の標準偏差ベクトル要素の最小値と同じくらいのオーダーであれば良いらしい。

fは$\boldsymbol{w}$の事前分布と尤度関数から自動的に求められる。

実際にサンプリングしてみる

あるデータ$X={x_1, x_2, ..., x_N}$が平均$\mu_X$、 標準偏差\sigma_Xガウス分布に従うとし、$\mu$をMCMCにより推定してみる(\sigmaは既知とする)。

事前準備

まず、適当にデータ$X$を作成する。

from scipy.stats import norm

mu_X = 5.0
sigma_X = 3.0
N = 50

X = norm.rvs(mu_X, sigma_X, size=N)

次に$\mu$の事前分布を設定する。 事前分布は平均$\mu_0=1.0$、標準偏差\sigma_0=10ガウス分布とした。

mu_0 = 1.0
sigma_0 = 10

def prior(mu):
    return norm.pdf(mu, mu_0, sigma_0)

関数q, sampling_q, f を設定する。 標準偏差\sigma_qは0.5とした。

sigma_q = 0.5

# 提案分布関数
def q(mu, mu_cond):
    return norm.pdf(mu, mu_cond, sigma_q)

# 提案分布から候補点を取得する関数
def sampling_q(mu):
    return norm.rvs(mu, sigma_q)

# wとXとの同時確率を計算する関数
def f(mu, X):
    # 尤度x事前分布
    return np.prod(norm.pdf(X, mu, sigma_X))*prior(mu)

サンプリング実行

サンプル数2000、$\mu$の初期値0としてサンプリングを実行する。 さらにサンプルの初めの方は初期値の影響を受けるのではじめの500サンプルは無視するようにした。

sample_num = 2000
warmup = 500 # はじめに無視するサンプル数
mu_init = 0.0

mu_samples = sampling_MH(X, q, sampling_q, f, sample_num, warmup, mu_init)

サンプリング結果

MCMCにより得られた$\mu$の事後分布サンプル結果を出力してみる。

import matplotlib.pyplot as plt

# muの事後分布
plt.hist(mu_samples, bins=50)
print('mean:', np.average(mu_samples))
print('standard deviation:', np.std(mu_samples))

$\mu$の事後分布サンプルは以下のようなヒストグラムとなった。

f:id:ricrowl:20200822025636p:plain
$\mu$の事後分布サンプルのヒストグラム

また、事後分布の平均と分散は以下となった。

mean: 4.360513325828551
standard deviation: 0.41838366124240634

さらに、$\mu$の事後分布は以下の式で解析的に求められるので検証のため解析的にも求めてみた。


{\displaystyle
\begin{align}
P(\mu|X)=N(\mu|\mu_N,\sigma^2_N)\\
\mu_N=\frac{\sigma^2_X}{\sigma^2_X+N\sigma^2_0}\mu_0+\frac{N\sigma^2_0}{\sigma^2_X+N\sigma^2_0}\bar{x}\\
\sigma^2_N=\frac{\sigma^2_0\sigma^2_X}{\sigma^2_X+N\sigma^2_0}\\
\bar{x}=\frac{1}{N}\sum^N_i x_i
\end{align}}

コードは以下。

X_mean = np.average(X)
eq_mu_N = (sigma_X**2*mu_0 + N*sigma_0**2*X_mean)/(sigma_X**2 + N*sigma_0**2)
eq_sigma_N = np.sqrt((sigma_0**2*sigma_X**2)/(sigma_X**2 + N*sigma_0**2))
print('mean:', eq_mu_N)
print('standard deviation:', eq_sigma_N)

plt.hist(mu_samples, bins=50, density=True, alpha=0.5)
plt_arr = np.linspace(min(mu_samples), max(mu_samples), 1000)
plt.plot(plt_arr, norm.pdf(plt_arr, eq_mu_N, eq_sigma_N), c='r')

得られた事後分布を先程のヒストグラムに重ねてみた。

f:id:ricrowl:20200822034621p:plain
青色:事後分布サンプルのヒストグラム、赤色:解析的に求めた事後分布

また、解析的に求めた事後分布の平均と分散は以下となった。

mean: 4.369898737127923
standard deviation: 0.42388274575892587

事後分布サンプルが結構いい感じ推定できてる👍

stanと比較

stanでもサンプリングして比較してみた。

import pystan

# モデル設計
code = '''

data {
    int N;
    real X[N];
    real sigma;
}

parameters {
    real mu;
}

model {
    for (i in 1:N) {
        X[i] ~ normal(mu, sigma);
    }
}

'''

dat = {
    'N': len(X),
    'X': X,
    'sigma': sigma_X
}

sm = pystan.StanModel(model_code=code)

# サンプリング実行
fit = sm.sampling(data=dat, iter=2500, warmup=500, thin=1, chains=1)

# 事後分布サンプル結果出力
la = fit.extract(permuted=True)
mu_samples_stan = la['mu']
print('mean:', np.average(mu_samples_stan))
print('standard deviation:', np.std(mu_samples_stan))
plt.hist(mu_samples_stan, bins=50, density=True, alpha=0.5, color='orange')
plt_arr = np.linspace(min(mu_samples_stan), max(mu_samples_stan), 1000)
plt.plot(plt_arr, norm.pdf(plt_arr, eq_mu_N, eq_sigma_N), c='r')

結果は以下となった。

f:id:ricrowl:20200822034650p:plainf:id:ricrowl:20200822034621p:plain
左:stanのサンプリング結果、右:実装コードのサンプリング結果

stanの結果 実装コードの結果
mean 4.354558596529226 4.369898737127923
standard deviation 0.4375811297845626 0.42388274575892587

stanも実装コードと似た結果が得られた👍

サンプルの時系列を見てみる

実装コードのサンプル時系列を見てみた。

f:id:ricrowl:20200822041455p:plain
実装コードのサンプル時系列。緑:受理されたサンプル、赤:棄却されたサンプル。

事後分布の平均となる4.36に近づくほど多く受理されて、離れるほど棄却されているように見える。これで事後分布のサンプルが得られるのがなんとなくわかる。

また、はじめの無視するサンプル数を500としたけど実際はずっと早く初期値の影響がなくなってるのでもっと少なくてよかったみたい。

さらに今度は提案分布の標準偏差を0.5から0.05に小さくしてみた。

f:id:ricrowl:20200822042529p:plainf:id:ricrowl:20200822042606p:plain
提案分布の標準偏差を0.05にした結果。上:サンプルの時系列、下:サンプルのヒストグラム

時系列が長い相関を持ってしまっている。。。この結果、サンプルも真の事後分布から大きく外れたものが得られてしまった。

次に提案分布の標準偏差を0.5から5に大きくしてみた。

f:id:ricrowl:20200822043256p:plainf:id:ricrowl:20200822043312p:plain
提案分布の標準偏差を0.05にした結果。上:サンプルの時系列、下:サンプルのヒストグラム

今度はサンプルのばらつきが大きすぎて多くのサンプルが棄却されてしまっている。そしてサンプルも山の部分が真の事後分布よりも低めのヒストグラムが得られてしまった。

今回真の事後分布の標準偏差は0.42...なので提案分布の標準偏差は0.5などのオーダーが適切であることが実際に確認できた。しかし桁数が1つでも違うと推定ができなくなるとは。。。事後分布の標準偏差がどのくらいのオーダーか当たりをつけて提案分布の標準偏差を設定しなければならないから結構シビアに感じる。。。

自作サーバーの性能比較(歴代RaspberryPi+PC)

家にいくつかのバージョンのRaspberryPi(とPC)があったのでサーバー化して性能比較してみた。

比較条件

比較ハード

  • Raspberry Pi Model B, Revision 2.0
  • Raspberry Pi3 Model B
  • Raspberry Pi4 Model B, 4GB RAM
  • ノートPC

※Pi2はないです...

ハード CPU メモリ 有線LAN USB OS
Raspberry Pi BCM2835
(700 MHz 1コア)
512 MB 10/100Mbps
イーサネット
USB2.0 Debian Buster
4.19
Raspberry Pi3 BCM2837
(1.2 GHz 4コア)
1 GB 10/100 Mbps
イーサネット
USB2.0 Debian Buster
4.19
Raspberry Pi4 BCM2711
(1.5 GHz 4コア)
4 GB Gigabit
イーサネット
USB3.0 Debian Buster
4.19
ノートPC Corei7-3632QM
(4コア8スレッド)
16 GB Gigabit
イーサネット
USB3.0 Ubuntu
20.04 LTS

その他

  • サーバーのストレージ:全てのハードで同条件にするため、USB3.0 & セルフパワーの外付けHDDをストレージに
  • サーバー化方法:sambaを使用(方法)。
  • クライアントPC(サーバーにアクセスするPC):Windows10

性能比較

  • CrystalDiskMark
  • 実際のファイル書き込み&読み込み時間
  • 使用感

で比較。

CrystalDiskMark

CrystalDiskMark(ver.7)を使用した。

CrystalDiskMarkを管理者権限で実行するとネットワークドライブのベンチマークが取れないため、「このアプリがデバイスに変更を加えることを許可しますか?」の画面で「いいえ」を選択して実行する。

f:id:ricrowl:20200714222448p:plain

  • Raspberry Pi3

f:id:ricrowl:20200714222602p:plain

  • Raspberry Pi4

f:id:ricrowl:20200714223001p:plain

  • ノートPC

f:id:ricrowl:20200714223021p:plain

ファイルサーバーの場合は上から1, 2行目のSEQ(シーケンシャルアクセス)を参考にすればいいみたい。

また、Q8T1などはそれぞれキュー数とスレッド数に対応していて、この場合はキュー数8, スレッド数1を表している。


Read・WriteともにPi1<Pi3<<Pi4<PCの順で性能が良かった。 特にWriteの差が大きい。

やっぱりGigabitイーサネットの差がでかいように見える。

実際のファイル書き込み&読み込み時間

1GBのファイルを実際にコピーして完了までの時間を計測した。

  • 1GBファイル作成方法
    下記コマンドをPowerShellで実行。
fsutil file createnew dummy_1GB.dat (1024*1024*1024)
@echo off

rem 読み込み速度:サーバーからクライアントへコピー
echo read
echo %time%
copy Z:\dummy_1GB.dat dummy_1GB_read.dat
echo %time%

rem 書き込み速度:クライアントからサーバーへコピー
echo write
echo %time%
copy dummy_1GB_read.dat Z:\dummy_1GB_write.dat
echo %time%

3回計測した平均が以下。

ハード 読み込み速度[MB/s] 書き込み速度[MB/s]
Raspberry Pi 5.43 2.04
Raspberry Pi3 10.99 11.03
Raspberry Pi4 100.00 30.58
ノートPC 103.09 107.53

読み込み速度はPi1<Pi3<<Pi4=PCの順で性能が良かった。CrystalDiskMarkのSEQ1M Q8T1(一行目)と似た結果。

書き込み速度はPi1<Pi3<Pi4<PCの順で、これもCrystalDiskMarkのSEQ1M Q8T1(一行目)と似た結果。

Pi1<Pi3<Pi4までは性能差が出るけど、Pi4とPCの読み込み速度は差が出ないみたい。

使用感

各ハードをサーバー化して使った時の使用感を比較した。

  • Raspberry Pi
    普通にエクスプローラでフォルダを開くときに固まったりする。
    普段使うサーバーとして使いたくない。。。めったに使わないファイル保管所としてはありか。。。?

  • Raspberry Pi3
    エクスプローラの操作やファイルアクセスには特に不満がない。
    ただし、ファイルのコピー時に2,3秒アクセスが遅れることがある。
    サーバーとして使えるけどたまにイラっとする。

  • Raspberry Pi4
    ファイルをいくつかコピーしているときでもアクセスが快適。
    サーバーとして安定して使える。

  • ノートPC
    使用感はRaspberry Pi4と変わらない。
    コピー時に速いなと思うくらいかな。

追加検証

USB2.0と3.0の比較

外付けHDDをUSB2.0接続か3.0接続かでサーバーの性能差が出るのか検証した。

Raspberry Pi4のUSB2.0ポートに外付けHDDを接続して性能を調べた。

f:id:ricrowl:20200714223001p:plainf:id:ricrowl:20200718171049p:plain
左:Pi4(USB3.0),右:Pi4(USB2.0)

ハード 読み込み速度[MB/s] 書き込み速度[MB/s]
Pi4(USB3.0) 100.00 30.58
Pi4(USB2.0) 100.00 29.15

CrystalDiskMarkのSEQ1M Q8T1のWrite以外は性能は変わらない。

実際のファイル書き込み&読み込み時間はどちらも性能に差なし。

イーサネットボトルネックになっていてUSB2.0でも3.0でも変わらないっぽい。

比較したサーバーの性能ってどんなものなの?

比較したサーバーの性能ってどんなものなのか調べるためにクライアントPCの内蔵SSDとクライアントPCにUSB接続した外付けHDDの性能も調べた。

f:id:ricrowl:20200718173437p:plainf:id:ricrowl:20200718173435p:plainf:id:ricrowl:20200718173431p:plain
左:内蔵SSD,中央:USB3.0接続外付けHDD,右:USB2.0接続外付けHDD

内蔵SSD圧倒的ィ...

Pi4やノートPCのサーバーは直接クライアントPCにUSB2.0接続した外付けHDD以上の性能はある。

まとめ

Raspberry Pi1, Pi3, Pi4, ノートPCをサーバー化して性能比較した。

だいたいPi1<Pi3<<Pi4<=PCの性能差で、個人用途ではPi4で支障なさそう。

個人的にはPi4にするかなー。

RaspberryPiのサーバー化


Sambaのインストール

以下のコマンドでインストールする。

sudo apt update
sudo apt install -y --fix-missing samba samba-common-bin

以下のような画面が出てくるので<はい>を選択して続行。 f:id:ricrowl:20200613163251p:plain:w400

以下のコマンドでSambaのステータスについての情報が表示されればOK。

sudo systemctl status smbd

共有ディレクトリ作成

適当な場所に共有したいディレクトリを作成しておく。

今回はhome下に共有ディレクトリShareを作成した。

mkdir /home/user00/Share

Sambaの設定

Sambaの設定ファイルを編集する。

sudo vim.tiny /etc/samba/smb.conf

まず、ファイルの最後に以下を追記。

[Share00]
comment = Share
path = /home/user00/Share
public = yes
read only = no
browsable = yes
force user = user00
guest ok = no
  • [Share00]: このサーバーの名前
  • path = /home/user00/Share: 共有ディレクトリのパス
  • force user = user00: ラズパイのユーザー名で、サーバーアクセス時このユーザーとしてアクセス

次にファイル上の[global]の直後に以下を追記。

map to guest = Never
security = user
guest account = nobody
browsable = no

一応特定のユーザーしかアクセスできないようにしている。

誰でもアクセスできるようにする場合はファイルの最後の追記内容のguest ok = noの箇所と[global]の追記を消す。

追記後、以下のコマンドでサーバーアクセスのパスワードを設定。

sudo smbpasswd -a user00

そして以下のコマンドでSambaを再起動。

sudo systemctl restart smbd

Windows上で認識

Windows10ではエクスプローラのネットワークを右クリックするとネットワーク ドライブの割り当て(N)という欄があるのでこれをクリックする。

出てきた画面で

  • ドライブ: 好きなドライブ文字
  • フォルダー: \\[サーバーのipアドレス]\[サーバー名]

を設定して完了を押す。

f:id:ricrowl:20200613182833p:plain:w300
ネットワーク ドライブの割り当て画面

ネットワーク資格情報の入力の画面が出てくるのでSambaで設定したユーザー名(この記事ではuser00)とパスワード(sudo smbpasswd -a user00で設定したパスワード)を入力する。

これでWindowsからサーバーにアクセスできるようになった。

RaspberryPiセットアップ(OS書き込み・SSH・ユーザー変更)

OS書き込み

ラズパイのSDやmicroSDカードにOSを書き込む。今回書き込むOSはRaspbian。

まずRaspbianの公式ページよりOSをダウンロードする。

推奨ソフトウェア入り(with desktop and recommended software),普通(with desktop),最小構成(Lite)があるので好きなものを選ぶ。

SDへの書き込みにはbalenaEtcherというソフトウェアを使う。WindowsMacLinuxで使えるので便利。

f:id:ricrowl:20200613020125p:plain:w300
balenaEtcher起動画面
balenaEtcherの画面の、Select imageでダウンロードしたOSのzipファイル、Select targetで書き込み先を指定してFlash!を押せばOK。

ラズパイに書き込んだカードを挿入して電源を入れればRasbianが起動する。

SSH設定

はじめに有線か無線でネットワークに繋いでおいてからラズパイを起動する。

起動させたRasbian上で端末を開き、以下をコマンドする。

sudo systemctl start ssh
sudo systemctl enable ssh

SSHが使用可能になったかは以下のコマンドで確認できる。

slogin localhost

Are you sure you want to continue connecting (yes/no)?というメッセージが出たらSSH使用可能。

yesを打ち込んで続行する。

追記:ラズパイ上で設定せずにSSHを有効化する方法(ヘッドレスインストール)

ラズパイ用のディスプレイがないときに便利。

OS書き込みの直後、microSDカードがbootという名前になっているのでboot直下にsshという空のファイルを作成する。

これによってラズパイ上で上記の設定をしなくてもOS書き込みの際にSSHを有効化することができる。

SSH接続のipアドレスルーターの設定画面かなんかで接続デバイスMACアドレスの一覧を見て当たりをつける。

ipアドレス確認/固定

PCからSSH接続するためにラズパイのipアドレスを確認もしくは固定する。

ipアドレスifconfigをコマンドして調べる。

ipアドレスを固定する場合は以下のようにしてdhcpcd.confに固定ipを書き込む。

sudo vim.tiny /etc/dhcpcd.conf

dhcpcd.confへは以下を追記する。

interface eth0 #無線の場合はinterface wlan0
static ip_address=[固定ip]/24
static routers=[ルーターのipアドレス]
static domain_name_servers=[ルーターのipアドレス]

ルーターipアドレスを知らない場合はifconfigで見たラズパイのipアドレス(192.168.XX.YY)の末尾を1にして192.168.XX.1にすればOK。

これ以降はSSHでラズパイを操作する。

ユーザーとパスワード変更

Raspbian のデフォルトでは

  • ユーザー名: pi
  • パスワード: raspberry

となっていてこのままではややこしいしセキュリティ的によろしくないので変更する。

一時ユーザー作成

SSHでpiにログインし、 piユーザーの名前を変更するために一時ユーザー(tmp)を作成する。

sudo useradd -M tmp
sudo gpasswd -a tmp sudo
sudo passwd tmp

一時ユーザーなのでとりあえずパスワードはpassとした。

デフォルトではpiユーザーに自動でログインしてしまうので起動時にユーザーを選択してログインするように変更する。

まず以下コマンドでraspi-configを起動する。

sudo raspi-config

すると以下の画面が出てくる。

f:id:ricrowl:20200613013742p:plain:w400
raspi-config起動画面

この画面上で

3 Boot Options Configure options for start-up

-> B1 Desktop / CLI Choose whether to boot into a desktop environm

と進んでいき、B1 Console Text console, requiring user to login を選択する。

選択後、<Finish> -> Would you like to reboot now? -> はい で終了し、再起動する。

piを新ユーザーに変更&パスワード変更

SSHでtmpにログインし、 以下コマンドでpiユーザーから新ユーザーへ名前を変更する。

sudo usermod -l [新ユーザー名] pi
sudo usermod -d /home/[新ユーザー名] -m [新ユーザー名]
sudo groupmod -n [新ユーザー名] pi

次に以下のコマンドで新ユーザーのパスワードを変更する。

sudo passwd [新ユーザー名]

これで変更完了。

後は一度exitでログアウトして新ユーザーでログイン後、以下のコマンドで一時ユーザーを削除する。

sudo userdel tmp