詳解 確率ロボティクス 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
実行すると,以下のように返ってくるはずです(長いので途中を省略しています).

date time ir lidar
0 20180122 95819 305 214
1 20180122 95822 299 211
2 20180122 95826 292 199
3 20180122 95829 321 208
4 20180122 95832 298 212
5 20180122 95835 327 212
         
... ... ... ... ...
         
         
58983 20180124 120023 313 208
58984 20180124 120026 297 200
58985 20180124 120030 323 204
58986 20180124 120033 326 207
58987 20180124 120036 321 208
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\}

ここで,[tex:{\bf{Z}}_{LiDAR}]はリスト(ベクトル)を表し,[tex: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):母集団から集団の一部を抽出すること