詳解 確率ロボティクス 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)$の期待値を考えると
[tex:<f(z)>_{p(z)} = \int_{- \infty}^{\infty} f(z)p(z) dz]
と定義できます.この式は,ある確率も出るから別の確率モデルのパラメータを求める時に頻出します.