2019年11月5日火曜日

詳解 確率ロボティクス 2.4 複雑な分布

前回の投稿に引き続き,詳解 確率ロボティクスの2.4 複雑な分布を通して学んだことやコードを実行してみた結果のメモです.この投稿はGithubでもご覧になれます(こちら).

sensor_data_600.txt*(ロボットから600[mm]離してlidar, irのセンサ値を取得)を用います.
%matplotlib inline
import pandas as pd
import matplotlib.pyplot as plt

data  = pd.read_csv("sensor_data_600.txt", delimiter=" "
                    header=None, names = ("date","time","ir","lidar"))

data["lidar"].hist(bins = max(data["lidar"]) - min(data["lidar"]),align='left')

plt.show()

二つの値が欠けるのは,センサ内部でディジタル処理をする様々なセンサで起こりうることです.
  • マルチモーダル(多峰性):ピークが2個以上ある分布
c.f. モードが二つの分布は特にバイモーダルと呼ばれます.モード(最頻値:統計用語)

全センサ値を時系列順に描画します(横軸:取得順にセンサ値につけた番号,縦軸:センサ値).

data.lidar.plot()
plt.show()


センサ値を時間毎にグループ分けして,各グループの平均値をグラフ化します.
data["hour"] = [e//10000 for e in data.time]  #hourly_mean,Note that "//" is not "/"#
d = data.groupby("hour")
d.lidar.mean().plot()

plt.show()

1行目:データフレームに hour という列を追加.値は時刻(時分秒の6桁)を10000で割ったもの(時の部分だけを残している)./ ではなく //であることに注意が必要です(// は小数点以下を切り捨てる割り算).
2行目:hour 列の値ごとに各レコードをグループ分けした新たなデータフレーム d を生成している.

6時台(オレンジ)と14時台(青)のヒストグラムを描画してみます.
d.lidar.get_group(6).hist()
d.lidar.get_group(14).hist()

plt.show()

ある変数$x$で条件づけられる別の変数$y$の確率分布は
$P(y|x)$
と表現されます(条件付き確率).

$P(z|t)$において,時刻$t$が分布を変えている直接の原因を表しているわけではありません(条件付き確率は変数間の直接の因果を表すものではない).

ここまでで,時刻を指定すると$P(z|t)$はガウス分布になるので,時間帯毎のガウス分布を作成します(センサ値$z$は整数,時刻$t$は$t$時台を表す整数).

センサ値が$a$, 時刻が$b$ということが同時に起こるという事象に対する確率$P (z = a, t = b)$を考えます($z = a$と$t = b$の同時確率).

each_hour = {i : d.lidar.get_group(i).value_counts().sort_index() for i in range(24)} # data frames every hour
freqs = pd.concat(each_hour, axis = 1) # Connect with concat
freqs = freqs.fillna(0) # Replace missing values(NaN) with 0
probs = freqs/len(data) # Convert frequency to probability

probs



グループ毎にセンサ値の頻度を集計して,全体の頻度で割ると確率になります.
import seaborn as sns
sns.heatmap(probs)

plt.show()

seabornのheatmapを使用した表示(横軸:時間,縦軸:センサ値).

sns.jointplot(data["hour"], data["lidar"], data, kind = "kde")
plt.show()

dataを使って描画したヒストグラム(seabornのjointplotを使用).

同時確率分布(統合確率分布)
図の右側は$P(z, t)$をセンサ値毎に表したグラフ.
$P(t) = \sum_{z} P(z, t)$

図の上側は$P(z, t)$を時間帯毎に表したグラフ.
$P(z) = \sum_{t} P(z, t)$

同時確率の表 probs から$P(z)$の値を計算してみると以下のようになる.
p_t = pd.DataFrame(probs.sum())
p_t.plot()
p_t.transpose()

p_t.sum()



p_t.sum()
とすると
0    1.0
dtype: float64
と返ってきます.

同時確率の表 probs から$P(z)$の値を計算してみると以下のようになります.
p_z = pd.DataFrame(probs.transpose().sum())
p_z.plot()
p_z.transpose()


p_z.sum()
とすると
0    1.0
dtype: float64
と返ってきます.

セル[9][10]で出力されたp_tとp_zの値は,セル[8]の図の上側,右側の分布と同じです.
上記のように,二つの変数に対しては,一般的に
$P(x) = \sum_{y \in \mathcal{Y}} P(x, y)$ ($\mathcal{Y}$: $y$の定義域)
が成り立ちます.

連続的な変数$x, y \in \mathcal{R}$に対しては
$p(x) = \int_{- \infty}^{+ \infty} p(x, y) dy$
が成り立ちます.これら(↑)の関係は,確率の加法定理と呼ばれます.

$P(x, y)$から変数(次元)を落として作った$P(x)$や$P(y)$の値は周辺確率と呼ばれます.
この操作を周辺化と呼び,周辺化してできた分布$P(x)$, $P(y)$は周辺分布と呼ばれます.
周辺確率からは消し去った変数の情報が消えます.

同時確率
  • 同時確率:$P(z, t)$と条件付き確率$P(z|t)$の関係
e.g. $t = 0$と固定したとします.
  •  $P(z|t = 0)$の分布は$P(z, t)$の分布の$t = 0$である部分と形状が同じになります.
  • $P(z, t)の$t = 0$の部分は$P(z|t = 0)$に対して$P(t = 0)$の分だけ値が割り引かれています.
任意の$t$に対して
$P(z|t) = P(z, t) / P(t)$ 条件付き確率
$P(z|t) = P(z, t) / P(t)$ 同時確率

これらの関係は確率の乗法定理と呼ばれます.

一般に
$P(x, y) = P(x | y) P(y) = P(y | x) P(x)$

条件付き確率の式
$P(z|t) = P(z, t) / P(t)$
に基づき,$P(z, t)$(下のコードではprobs)から$P(z|t)$を作っています.
cond_z_t = probs/p_t[0]
cond_z_t

(cond_z_t[6]).plot.bar(color = "blue", alpha = 0.5)
(cond_z_t[14]).plot.bar(color = "orange", alpha = 0.5)
plt.show()

上図は,変数 cond_z_t から確率分布$P(z|t = 6)$,$P(z|t =14)$を出力したものです.

同時確率$P(z, t)$と条件付き確率$P(z|t)$の関係

e.g. $t = 0$と固定する.
  • $P(z|t = 0)$の分布は$P(z, t)$の分布の$t = 0$である部分と形状が同じになる.
  • $P(z, t)の$t = 0$の部分は$P(z|t = 0)$に対して$P(t = 0)$の分だけ値が割り引かれている.
任意の$t$に対して
$P(z|t) = P(z, t) / P(t)$ 条件付き確率
$P(z|t) = P(z, t) / P(t)$ 同時確率

これらの関係は確率の乗法定理と呼ばれます.

一般に
$P(x, y) = P(x | y) P(y) = P(y | x) P(x)$

が成り立ち,この関係を使うと,下式の周辺化の計算
$P(x) = [P(x, y)]_{y}$

を以下のように行うことができます.
$P(x) = [P(x, y)]_{y} = [P(x|y) P(y)]_{y} = <P(x|y)>_{P(y)}$

同様に,
$p(x) = [p(x, y)]_{y}$

に対しても
$p(x) = [p(x, y)]_{y} = [p(x|y)]_{p(y)}$

と計算できます.

独立,従属,条件付き独立
乗法定理の式
$P(x, y) = P(x | y) P(y) = P(y | x) P(x)$
において,もし,事象$y$が事象$x$に全く影響を与えない場合,
$P(x|y) = P(x)$
が成り立ちます.

互いに無関係な事象は,
$x \Perp y$
上記のような記号で表現され,事象$x$と事象$y$が独立と呼ばれます(独立でない事象は従属と呼ばれる).

$\Prep$はTeXの"txfonts"において独立性を表す記号だが,うまく表示されていません....

ある条件(今回の場合は時間帯)に限ると独立が成り立つ場合は条件付き確率と同様に縦棒 | を使って
$x \Perp y | t$
と表記します.

このような条件付きの独立関係は条件付き独立と呼ばれ,上式($x \Perp y | t$)が成り立つ場合は
$P(z, y |t) = P(z|t)P(y|t)$
が成り立ちます.

確率分布の性質を利用した計算
加法定理,乗法定理,独立を利用すると確率分布の絡んだ積分を簡単にすることができます.

例えば,$x \Prep y$,$x \in R$,$y \in R$のとき
$\bf{z} = \left(
    \begin{array}{c}
      x \\
      y \\
    \end{array}
  \right)$
というベクトル${\bf{z}} \in R^{2}$を考えて,次の積分を考えます.
$\int_{{\bf{z}} \in R^{2}} p({\bf{z}}) \{ f(x) + \alpha g(y) \} d{\bf{z}} = [ p({\bf{z}}) \{ f(x) + \alpha g(y) \} ]_{{\bf{z}}}$
この式は$x$と$y$の独立性を考えると,
$[p({\bf{z}}) \{ f(x) + \alpha g(y) \}]_{{\bf{z}}} = < f({\bf{z}}) + \alpha g(y)>_{p({\bf{z}})} = <f(x) +\alpha g(y)>_{p(x)} = <f(x)>_{p(x)p(y)} + \alpha g(y)>_{p(x)p(y)}$
と変形できます.
最後の式を見ると,左側の項の$p(x)$と,右側の項の$p(x)$は消去でき,
$[p({\bf{z}}) \{ f(x) + \alpha g(y) \}]_{{\bf{z}}} = < f(x) >_{p(x)} + \alpha <g(y)>>_{p(y)}$
となります.

分布を消去できる理由
期待値$<f(x)>_{p(x)p(y)}$は,分布$p(x)p(y)$から$x, y$を無限にドローして,$f(x)$の平均値を取った値です.
この時,ドローされた$y$の値は一切使われません.
このこと(↑)を考えると,期待値$<f(x)>_{p(x)p(y)}$は$x$だけ$p(x)$からドローして求めた期待値$<f(x)>_{p(x)}$と一致します.

i.e.
$𝑥 \Prep 𝑦 \Rightarrow <f(x)>_{p(x)p(y)} = <f(x)>_{p(x)}$

続いて,以下のような積分を考えます.
$\int_{R^{2}} p({\bf{z}}) f(x) g(y) d {\bf{z}} = < f(x) g(y) >_{p({\bf{z}})}$
ここで,${\bf{z}} = ( x \ y )^{T}$

${\bf{z}}$の積分を$x, y$の積分に分解して良いならば,左辺の積分は,
$\int_{R^{2}} p({\bf{z}}) f(x) g(y) d {\bf{z}} = \int_{R} \int_{R} p(x) p(y) f(x) g(y) dy dx$

$= \int_{R} p(x) \int_{R} p(y) f(x) g(y) dy dx$
$= \int_{R} p(x) <f(x) g(y)>_{p(y)} dx$
$= <<f(x) g(y)>_{p(y)}>_{p(x)}$
$= <f(x) <g(y)>_{p(y)} >_{p(x)}$
$= <g(y)>_{p(y)} <f(x)>_{p(x)}$

となります.

この結果を$x$, $y$の順を変えて遡ると
$<g(y)>_{p(y)} <f(x)>_{p(x)} = <<f(x)g(y)>_{p(x)}>_{p(y)}$
となります.
これらの計算から,$x \Prep y$のとき
$<g(y)>_{p(y)} <f(x)>_{p(x)} = <<f(x) g(y)>_{p(x)}>_{p(y)}$
$= <<f(x) g(y)>_{p(y)}>_{p(x)}$
$= <f(x) g(y)>_{p(x)p(y)}$
が成り立ちます.
この式は,$x$, $y$のどちらか一方,あるいは両方をベクトルに変えても成立します.

ベイズの定理
乗法定理の式
$P(x, y) = P(x|y)P(y) = P(y|x)P(x)$
の中辺,右辺からは下式が導出できる(ベイズの定理).
$P(x | y) = \cfrac{P(y|x)P(x)}{P(y)}$

$x$の取り得る値の集合を$\mathcal{X}$とすると,加法定理から
$P(y) = \sum_{x \in \mathcal{X}} P(x, y)$
が成り立つので,これを
$P(x | y) = \cfrac{P(y|x)P(x)}{P(y)}$
に代入すると,
$P(x | y) = \cfrac{P(y|x)P(x)}{\sum_{x' \in \mathcal{X}} P(x', y)} = \cfrac{P(y|x)P(x)}{\sum_{x' \in \mathcal{X}} P(y | x') P(x')} = \cfrac{P(y|x)P(x)}{<P(y|x')>_{P(x')}}$
が成り立ちます(ベイズの定理).

ベイズの定理は$x$,$y$が連続でも成立ち,さらにはベクトルでも成立ちます.

連続なベクトル${\bf{x}}$,${\bf{y}}$を使ってベイズの定理を表記すると$x \in \mathcal{X}$のとき
$p({\bf{x}} | {\bf{y}}) = \cfrac{p({\bf{y}} | {\bf{x}}) p({\bf{x}})}{{\bf{y}}} = \cfrac{p({\bf{y}} | {\bf{x}}) p({\bf{x}})}{\int_{\mathcal{X}} p({\bf{y}} | {\bf{x'}}) p({\bf{x'}}) d{\bf{x'}}} = \cfrac{p(\bf{y}|\bf{x})p(\bf{x})}{<p(\bf{y}|\bf{x'})>_{p(\bf{x'})}}$
となる.

さらに,$y$が決まっていると
$P(x | y) = \cfrac{P(y|x)P(x)}{\sum_{x' \in \mathcal{X}} P(x', y)} = \cfrac{P(y|x)P(x)}{\sum_{x' \in \mathcal{X}} P(y | x') P(x')} = \cfrac{P(y|x)P(x)}{<P(y|x')>_{P(x')}}$
の分母は定数なので,
$P(x | y) = \eta P(y|x)P(x)$
と表記することもあります.
上記の式
$P(x | y) = \cfrac{P(y|x)P(x)}{\sum_{x' \in \mathcal{X}} P(x', y)} = \cfrac{P(y|x)P(x)}{\sum_{x' \in \mathcal{X}} P(y | x') P(x')} = \cfrac{P(y|x)P(x)}{<P(y|x')>_{P(x')}}$
の場合は,
$P({\bf{x}} | {\bf{y}}) = \eta P({\bf{y}}|{\bf{x}})P({\bf{x}})$
となる.

$\eta$に関しては,最終的に分布を積分したときに1になるように$\eta$を選べば良い((確率分布の)正規化)という考え方です.
$\eta$は正規化定数と呼ばれます.

分布の形状に影響を与えない変数,定数を全て$\eta$に押し込むという使われ方なので,途中で値が変わることもあります.

$P(z = 630 | t = 13)$について,セル[12]の cond_z_t で計算したものと,ベイズの定理(下式)から計算したものを比較します.
$P(x | y) = \cfrac{P(y|x)P(x)}{P(y)}$

cond_t_z = probs.transpose() / probs.transpose().sum()

print("P(z = 630) = ", p_z[0][630])
print("P(t = 13) = ", p_t[0][13])
print("P(t = 13 | z = 630) =", cond_t_z[630][13])
print("Bayes P(z = 630 | t = 13) = ", cond_t_z[630][13] * p_z[0][630] / p_t[0][13])
print("answer P(z = 630 | t = 13) = ", cond_z_t[13][630])

とすると,

P(z = 630) =  0.06694936878045224
P(t = 13) =  0.043024993620976656
P(t = 13 | z = 630) = 0.023230490018148822
Bayes P(z = 630 | t = 13) =  0.036147980796385204
answer P(z = 630 | t = 13) =  0.036147980796385204

という結果が返ってきて,計算結果が同じであることがわかります.

ベイズの定理を使うと,得られているデータから原因を推定することが可能になります.
sensor_dataディレクトリの中には様々なセンサ値(距離)が含まれています.

続いて,「センサ値から計測した時間帯を推定する問題」を考えます.
以下の関数は,時間帯tごとに$P(t)$の値を受け取り,4行目で$P(t|z) = \eta P(z|t)P(t)$の$P(z|t)P(t)$を計算して,6行目で正規化して$P(t|z)$を返す関数です.

def bayes_estimation(sensor_value, current_estimation): 
    new_estimation = []
    for i in range(24):
        new_estimation.append(cond_z_t[i][sensor_value]*current_estimation[i])
        
    return new_estimation / sum(new_estimation)

センサ値 630 が得られた時の$P(t|z)$は次のように計算できます.

estimation = bayes_estimation(630, p_t[0])
plt.plot(estimation)


630の次に,632, 636というセンサ値が得られたとして,これらの値を$z_{1}$,$z_{2}$,$z_{3}$と表すと,ベイズの定理より
$P(t | z_{1}, z_{2}, z_{3}) = \eta P(z_{1}, z_{2}, z_{3} | t) P(t)$
となります($z_{1}, z_{2}, z_{3}$をひとまとめの事象とみなしている).

もし,$z_{1}, z_{2}, z_{3}$が時間帯$t$を固定したときに互いに独立だとすると,下式
$P(z, y |t) = P(z|t)P(y|t)$
から,
$P(t | z_{1}, z_{2}, z_{3}) = \eta P(z_{1}, z_{2}, z_{3} | t)P(t) = \eta P(z_{1} | t) P(z_{2}, z_{3} | t) P(t) = \eta P(z_{1} | t) P(z_{2} | t) P(z_{3} | t) P(t)$
となり,$P(z|t)$を掛け算していくだけで,$z_{1}$,$z_{2}$,$z_{3}$を得た値$t$を推定していくことができます.

value_5 = [630, 632, 636]

estimation = p_t[0]
for v in value_5:
    estimation = bayes_estimation(v, estimation)

plt.plot(estimation)

sensor_data_600.txt からある日の11時台のデータを三つ連続で選んでみます.

values_11 = [617, 624, 619]

estimation = p_t[0]
for v in values_11:
    estimation = bayes_estimation(v, estimation)
    
plt.plot(estimation)


2019年11月4日月曜日

詳解 確率ロボティクス 2.3 確率モデル

前回の投稿に引き続き,詳解 確率ロボティクスの2.3 確率モデルを通して学んだことやコードを実行してみた結果のメモです.この投稿はGithubでもご覧になれます(こちら).

ガウス分布*は(例えば)センサ値zがa以上b未満($[a, b)$)に入る確率を
$P (a \leq z < b) = \int_a^b p(z) dz$
ここで,
$p(z) = \cfrac{1}{\sqrt{2\pi\sigma^{2}}}\exp \left\{ - \cfrac{(z - \mu)^{2}}{2 \sigma^{2}} \right\}$
を表します. また,$\sigma^{2}$は分散,$\mu$は平均値を表します.

 *ガウス分布と正規分布は同義です.

センサ値の平均値$\mu = 209.7[mm]$,分散$\sigma^{2} = 23.4$を代入して$P, p$を描画するには以下のようにします.
import math

def p(z, mu = 209.7, dev = 23.4):
    return math.exp(-(z -mu)**2 / (2*dev))/math.sqrt(2*math.pi*dev)


zs = range(190, 230)
ys = [p(z) for z in zs]

import matplotlib.pyplot as plt
%matplotlib inline
plt.plot(zs, ys)

plt.show()
zsは横軸の数値のリスト,ysは縦軸の関数pの値のリストです.

さらに,$p$を積分して,センサ値が整数に限定される場合の確率分布を作ってみます.
センサ値$x$に対して,区間$[ x - 0.5, x + 0.5)$の範囲で積分することにします(ここでは台形公式で近似).
def prob(z, width = 0.5):
    return width*(p(z - width) + p(z + width))

zs = range(190, 230)
ys = [prob(z) for z in zs]

import pandas as pd
data = pd.read_csv("sensor_data_200.txt", delimiter = " ",
                   header = None, names = ("data", "time", "ir", "lidar"))
freqs = pd.DataFrame(data["lidar"].value_counts())
freqs["probs"] = freqs["lidar"] / len(data["lidar"])

plt.bar(zs, ys, color = "red", alpha = 0.3) # Make graphs transparent with alpha
f = freqs["probs"].sort_index()
plt.bar(f.index, f.values, color = "blue", alpha = 0.3)

plt.show()
ガウス分布は$\mu$, $\sigma^{2}$を決めると形状が決まるので,式を記述する必要がなければ
$\mathcal{N}(z|\mu, \sigma^{2})$ \ あるいは $\mathcal{N}(\mu, \sigma^{2})$
などと略記します.
  • モデル化:ある現象を説明するために適切な確率分布の数式を用いてパラメータを決めること
  • 確率モデル:モデル化で分布に当てはめられる数式
ガウス分布の確率密度関数(probability density function, pdf)
$p(z) = \cfrac{1}{\sqrt{2\pi\sigma^{2}}}\exp \left\{ - \cfrac{(z - \mu)^{2}}{2 \sigma^{2}} \right\}$
  • 確率密度関数は,積分すると確率になる関数
  • 確率密度関数の値は密度と呼ばれる
ガウス分布などの確率密度関数を扱うにはSciPyが便利です.
このモジュール下にあるサブモジュールstatsにはガウス分布の確率密度関数norm.pdfが実装されています.
mean1 = sum(data["lidar"].values)/len(data["lidar"].values)

zs = data["lidar"].values
mean = sum(zs) / len (zs)
diff_square = [(z -mean)**2 for z in zs]
sampling_var = sum(diff_square)/(len(zs))
stddev1 = math.sqrt(sampling_var)

from scipy.stats import norm

zs = range(190, 230)
ys = [norm.pdf(z, mean1, stddev1) for z in zs]

plt.plot(zs, ys)
plt.show()
変数$z$が実数のとき,確率密度関数$p$を積分したものは累積分布関数(cumulated distribution function, cdf)と呼ばれます.
$p(z < a) = \int_{-\infty}^{a} p(z) dz$
zs = range(190, 230)
ys = [norm.cdf(z, mean1, stddev1) for z in zs]

plt.plot(zs, ys, color = "red")
plt.show()
先に台形公式で近似計算を行なった下式の計算は,確率の計算に置き換えられます.
$P(a \leq z < b) = \int_{a}^{b} p(z) dz = \int_{-\infty}^{b} p(z) dz - \int_{-\infty}^{a} p(z) dz = P(z < b) -P(z < a)$
上式を使って確率分布を描くには以下のようにします.
zs = range(190, 230)
ys = [norm.cdf(z + 0.5, mean1, stddev1) - norm.cdf(z - 0.5, mean1, stddev1) for z in zs]

plt.bar(zs, ys)
plt.show()
確率密度関数は,確率質量関数が大文字$P$で表記されるのに対して,多くの場合,区別のために小文字$p$と表記されます.


  • 期待値:確率分布$P$について$z \sim P(z)$を無限に繰り返した場合に,$z$の平均値がどれくらいになるかを表す値(確率分布$P$から$z$を発生させる).

$z$が離散的な場合,期待値は
$\sum_{z = \infty}^{\infty} z P(z)$
で計算できます.

$z$が連続的な場合は
$\int_{- \infty}^{\infty} zP(z) dz$
で計算できます.
期待値は具体的な値をドローしなくても分布が決まっていれば,下式のいずれか(離散的な場合は上の式,連続的な場合は下の式)
$\sum_{z = \infty}^{\infty} z P(z)$

$\int_{- \infty}^{\infty} zP(z) dz$
で計算できます(定義通りに何回もドローして値をサンプリングし,平均を取ることでも期待値を近似的に求めることが可能).
$z \sim p(z)$や,$z \sim P(z)$のとき,$z$の期待値は
$E_{p(z)} [z], \ E_{P(z)} [z]$ あるいは $_{p(z)}, \ _{P(z)}$
と表記されます.

期待値をさらに一般化して,$z \sim p(z)$から計算される関数の値$f(z)$の期待値を考えると
$<f(z)>_{p(z)} = \int_{- \infty}^{\infty} f(z)p(z) dz$
と定義できます.この式は,ある確率も出るから別の確率モデルのパラメータを求める時に頻出します.

2019年11月1日金曜日

詳解 確率ロボティクス 2.2 度数分布と確率分布

UASの自律飛行にあたって,ロボティクスに興味を持ちました.上田隆一先生の書籍 詳解 確率ロボティクスが出版され,非常に分かり易そうだったので,詳解 確率ロボティクスの2.1センサデータの収集とJupyter Notebook場での準備,2.2 度数分布と確率分布を通して学んだことやコードを実行してみた結果のメモです.この投稿はGithubでもご覧になれます(こちら).
  • 頻度:ヒストグラムの縦軸の値のこと.
  • 雑音(ノイズ):計測値(センサ値)の変動(e.g. LiDARの計測値).例えば,LiDARの場合は,外乱光や電気回路中の電圧や電流を乱す"何か(原因は複雑にあって,互いに影響し合っていることが多い)"が影響して雑音が発生する.
  • 誤差:何らかの測りたいものの"真の値"とセンサ値の差のこと.
  • 偶然誤差(accidental error, random error):雑音によって発生する誤差のこと.
  • 偏り(バイアス):計測機器の取り付け位置がずれていたり,取り付けた本体が傾いていたりするときに生じるずれのこと.
  • 系統誤差(Systematic error):バイアスによって生じる定常的な誤差のこと.系統誤差の量はセンサ値から推定することができないので,別のセンサや計測方法で突き止める必要がある.しかし,別の計測方法やセンサにも雑音やバイアスが存在する.
  • 系統誤差は,アルゴリズムの出力にも悪影響を及ぼす.対策は取りにくいが,バイアスや系統誤差の存在を頭の隅に置いて,ロボットのアルゴリズムを考えていく必要がある.
書籍 詳解 確率ロボティクスの中では,Python 3系とJupyter Notebookを用いてコードの実行結果を確認する方法が取られています.コードは著者の上田隆一先生のGitHubで公開されています.書籍中のコードの閲覧はこちらから可能です.完成したコードのみを見たい場合はこちらから閲覧が可能です.

以下は,詳解 確率ロボティクス 2.2.2 頻度,雑音,バイアス にあるコードを実行してみた際のメモです.

まずは,Pythonのバージョン等のチェック(Python 2系だとコードが動かないため,3系であることを確認)
import sys
sys.version    # Check Python version

実行すると,以下のように返ってくるはずです(バージョンなどは環境によって異なります).

'3.7.5 (default, Nov 1 2019, 02:16:32) \n[Clang 11.0.0 (clang-1100.0.33.8)]'

続いて,Pandasモジュールの"read_csv"関数を使って data という変数に既存の計測データ(ファイル名:sensor_data_200.txt, 計測した日付,時間,光センサの計測値,LiDARの計測値)*を読み込んでいます.
import pandas as pd
data = pd.read_csv("sensor_data_200.txt", delimiter = " "
                  header = None, names = ("date", "time", "ir", "lidar"))
data

実行すると,以下のように返ってくるはずです(長いので途中を省略しています).

datetimeirlidar
02018012295819305214
12018012295822299211
22018012295826292199
32018012295829321208
42018012295832298212
52018012295835327212

...............

5898320180124120023313208
5898420180124120026297200
5898520180124120030323204
5898620180124120033326207
5898720180124120036321208
58988 rows × 4 columns

続いて,sensor_data_200.txt における,LiDARのセンサ値値の規則性の有無を調べるために,センサ値のヒストグラムの描画を行います.
%matplotlib inline    #necessary if the graph is not displayed.
import matplotlib.pyplot as plt
data["lidar"].hist(bins = max(data["lidar"]) - min(data["lidar"]), align = 'left')
plt.show()

実行すると,以下のようなヒストグラムが描画されます.


* 計測データの中身は以下のようになっています.
20180122 095819 305 214
20180122 095822 299 211
20180122 095826 292 199
20180122 095829 321 208
20180122 095832 298 212
...途中省略...
20180124 120023 313 208
20180124 120026 297 200
20180124 120030 323 204
20180124 120033 326 207
20180124 120036 321 208

雑音の原因を突き止めたり除去するのは困難な場合が多いので,原因はわからないままにして,雑音の傾向を把握する(分からないことは放っておく).
傾向を把握するためには以下の値を求めます.
  1. センサ値の平均値を求める.
  2. センサ値の分散を求める.
  3. センサ値の標準偏差を求める.
なお,分散と標準偏差は各センサ値のばらつきを表します.
${\bf{Z}}_{LiDAR} = \{ z_{i} | i = 0, 1, 2, \cdots, N-1\}$
ここで,${\bf{Z}}_{LiDAR}$はリスト(ベクトル)を表し,$z_{i}$の$i$は(0からの)要素番号を表し,要素数は$N$となります.
  • 分散:集合の要素の値を全て足して,要素数で割ったもの
$\mu = \cfrac{1}{N} \sum_{i = 0}^{N-1} z_{i}$
  • 分散には標本分散と不偏分散の二つがある.
  • 標本分散:各値と平均値($\mu$)の差を二乗したものの平均値(各値と平均値($\mu$)が離れているほど大きい)
$\sigma^{2} = \cfrac{1}{N} \sum_{i = 0}^{N-1} (z_{i} - \mu)^{2} \ \ \ (N > 0)$
  • 不偏分散標本分散の割り算の分母が N ではなく N-1 としたもの.
$(s^{2} = )\sigma^{2} = \cfrac{1}{N-1} \sum_{i = 0}^{N-1} (z_{i} - \mu)^{2} \ \ \ (N > 1)$
  • 標本分散と不偏分散は,Nの値が小さい時に違いが現れる.
上記の分散を実際に計算するコードは以下のようになります.
# Calculate from definition
zs = data["lidar"].values
mean = sum(zs) / len(zs)
diff_square = [(z - mean)**2 for z in zs]

sampling_var = sum(diff_square) / (len(zs))     # sample variance
unbiased_var = sum(diff_square) / (len(zs) - 1) # unbiased sample variance

print(sampling_var)
print(unbiased_var)

# Use Pandas
pandas_sampling_var = data["lidar"].var(ddof = False) # sample variance
pandas_default_var = data["lidar"].var()              # default(unbiased sample variance)

print(pandas_sampling_var)
print(pandas_default_var)

# Use Numpy
import numpy as np

numpy_default_var = np.var(data["lidar"])            # default(sample variance)
numpy_unbiased_var = np.var(data["lidar"], ddof = 1) # unbiased sample variance

print(numpy_default_var)
print(numpy_unbiased_var)

最初の# Calculate from definition以下は,定義から計算するもので,次の# Use Pandas以下はPandasを使って計算するもの,最後の# Use Numpy以下はNumpyを使って計算するものです.各計算結果を表示するように命令してあり,実行結果は以下のようになります.
この結果からPandasとNumpyを使った結果は一致していることがわかります.Pythonでは,上記のように定義から計算するよりもPandasやNumpyを使って計算することが推奨されており,計算速度もそちらの方が速いようです.
計算結果は以下のようになります.

23.407709770274106
23.40810659855441
23.4077097702742
23.408106598554504
23.4077097702742
23.408106598554504

続いて標準偏差を求めます(標準偏差は分散の正の平方根).
import math

# Calculate from definition
stddev1 = math.sqrt(sampling_var)
stddev2 = math.sqrt(unbiased_var)

# Use Pandas
pandas_stddev = data["lidar"].std() # Pandas calculates standard deviation from unbiased variance

print(stddev1)
print(stddev2)

print(pandas_stddev)

計算結果は以下のようになります.
4.838151482774605
4.83819249292072
4.838192492920729


計算結果から,Pandasでは不偏分散を用いて標準偏差を求めていることがわかります.
  • 確率:値の出やすさを数値化したもの
まず,各センサ値の頻度を集計してみます.以下では,value_countsでlidar列のセンサ値の頻度を数え上げて,pd.DataFrameでデータフレームにしています.
freqs = pd.DataFrame(data["lidar"].value_counts())
freqs.transpose() # Output horizontally

集計結果は以下のようになります.

1 rows × 35 columns

続いて,変数freqsに確率の列を追加してみます.1行目では,lidar列に入っているそれぞれの頻度を,dataの要素数で割っています.
freqs["probs"] = freqs["lidar"]/len(data["lidar"])
freqs.transpose()

2 rows × 35 columns
(表は途中-右側-を省略しています)

確率の合計が1になっていることを確認します.
sum(freqs["probs"])

合計の計算結果は以下のようになります.
1.0

出力結果を並べ替えて,横軸にセンサ値を,縦軸に確率を描いてみます.
ヒストグラムと似ていますが,縦軸は頻度ではなく確率であることに注意が必要です.
freqs["probs"].sort_index().plot.bar()
plt.show() # The vertical axis changes from frequency to establishment

描画結果は以下のようになります.

  • 確率質量関数:個別の確率$P(z)$を与える関数$P$
  • 確率分布:各変数に対して確率がどのように分布するのかを表す実体
最後に,確率分布を用いたシミュレーションを行ってみた結果です.
シミュレーションでは,先に求めた確率分布($P_{{\bf{Z}}LiDAR}$と記する)に従って$z$を選びます.Pandasでは,sampleメソッドを使用することで確率分布から値を選ぶことができます.
上記の処理を数式で表すと
$z_{N} \sim P_{{\bf{Z}}LiDAR}$
となります.
ここで,左辺の$z_{N}$は実際に選ばれた値を表します.

def drawing(): # Define as a function
    return freqs.sample(n = 1, weights = "probs").index[0]

drawing() # meaning of execution, not drawing graph

実行結果は以下のようになります.
211

sampleの引数nが選ぶ個数,weights = "probs"がデータフレームの"probs"の列に選ぶときの確率が入っていることを意味しています.sampleの後ろの.index[0]は,データフレームのレコードの名前(この場合はセンサ値)を取り出すためのもので,これによってセンサ値を得ることができます.
上記の処理を数式で表すと
$z_{N} \sim P_{{\bf{Z}}LiDAR}$

となります.
ここで,左辺の$z_{N}$は実際に選ばれた値を表します.

Pandasでsamplesメソッドを使って確率分布から値を選ぶことが可能です.N-1回目までのセンサ値で作った分布から$z_{N}$を発生させるには以下のようにします.

samples = [drawing() for i in range(len(data))]
# samples = [drawing() for i in range(100)]
simulated = pd.DataFrame(samples, columns = ["lidar"])
p = simulated["lidar"]
p.hist(bins = max(p) - min(p), color = "orange", align = 'left')
plt.show()

描画結果は以下のようになります.
  • ドローイング(drawing):母集団から個々のものを抽出すること
  • サンプリング(sampling):母集団から集団の一部を抽出すること