Probabilistic Robotics（確率・統計の基礎）

永井 忠一 2020.10.24

histogram の描画

実験データは、こちらにある『sensor_data/sensor_data_200.txt』を使用

ソースコード
PythonR言語
import pandas as pd

data = pd.read_table('sensor_data_200.txt', delimiter=' ',
                     header=None, names=('date', 'time', 'ir', 'lidar'))

import matplotlib.pyplot as plt

plt.clf()
plt.hist(data.lidar, bins=(max(data.lidar) - min(data.lidar) + 1))
plt.title('LiDAR (200mm)')
plt.xlabel('Z [mm]')
plt.ylabel('Frequency')
plt.savefig('fig1.png')
 
data <- read.table('sensor_data_200.txt', sep=' ')
colnames(data) <- c('date', 'time', 'ir', 'lidar')

png('fig2.png')
hist(data$lidar, main='LiDAR (200mm)', xlab='Z [mm]',
     breaks=seq(min(data$lidar), max(data$lidar)))
dev.off()

平均値、分散、標準偏差

DataFrame/Series の describe() と、R言語の summary() 関数を利用

IPythonEmacs Speaks Statistics (ESS)
In [1]: import pandas as pd                                                     

In [2]: data = pd.read_table('sensor_data_200.txt', delimiter=' ', header=None, 
   ...: names=('date', 'time', 'ir', 'lidar'))                                  

In [3]: data.lidar.describe()
Out[3]: 
count    58988.000000
mean       209.737133
std          4.838192
min        193.000000
25%        206.000000
50%        210.000000
75%        213.000000
max        229.000000
Name: lidar, dtype: float64


 > data <- read.table('sensor_data_200.txt', sep=' ',
 + col.names=c('date', 'time', 'ir', 'lidar'))
!> summary(data$lidar)
    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
   193.0   206.0   210.0   209.7   213.0   229.0

扱っているセンサーの平均値は 209.7[mm]、標準偏差は 4.8[mm]

プログラムと実行結果

PythonR言語
import pandas as pd

data = pd.read_table('sensor_data_200.txt', delimiter=' ',
                     header=None,
                     names=('date', 'time', 'ir', 'lidar'))

print('平均値:', round(data.lidar.mean(), 1))
print('中央値:', int(data.lidar.median()))
print('最頻値:', int(data.lidar.mode()))
print('分散:', round(data.lidar.var(), 1))
print('標準偏差:', round(data.lidar.std(), 1))
 
data <- read.table('sensor_data_200.txt', sep=' ')
colnames(data) <- c('date', 'time', 'ir', 'lidar')

print(paste('平均値:', round(mean(data$lidar), 1)))
print(paste('中央値:', median(data$lidar)))
print('最頻値:') # ?!?
print(paste('分散:', round(var(data$lidar), 1)))
print(paste('標準偏差:', round(sd(data$lidar), 1)))
 
平均値: 209.7
中央値: 210
最頻値: 211
分散: 23.4
標準偏差: 4.8

[1] "平均値: 209.7"
[1] "中央値: 210"
[1] "最頻値:"
[1] "分散: 23.4"
[1] "標準偏差: 4.8"

※ライブラリの分散は「不偏分散」を求める。（サンプル分散ではない）

（R言語には、最頻値を求める関数が見当たらない。）

確率分布

頻度ではなく割合を考える \( P_{\mathbf{z}\text{LiDAR}}(z) = N_z / N \)（\( N_z \): センサー値が \( z \) だった頻度、\( N \): 実験の回数）

ソースコード
PythonR言語
import pandas as pd

data = pd.read_table('sensor_data_200.txt', delimiter=' ',
                     header=None,
                     names=('date', 'time', 'ir', 'lidar'))

import matplotlib.pyplot as plt

plt.clf()
plt.hist(data.lidar, bins=(max(data.lidar) - min(data.lidar) + 1),
         density=True)
plt.title('P_zLiDAR(z)')
plt.xlabel('Z [mm]'); plt.ylabel('Density')
plt.savefig('fig3.png')
 
data <- read.table('sensor_data_200.txt', sep=' ')
colnames(data) <- c('date', 'time', 'ir', 'lidar')

png('fig4.png')
hist(data$lidar,
     main=expression(P[zLiDAR](z)), xlab='Z [mm]',
     breaks=seq(min(data$lidar), max(data$lidar)),
     freq=FALSE)
dev.off()

この関数を確率と呼ぶ（確率の定義）（縦軸の値が違うだけ）（histogram を正規化したもの）

期待値の性質

  • 線形性
    • \( \left\langle f(z) + \alpha g(z) \right\rangle_{p(z)} = \left\langle f(z) \right\rangle_{p(z)} + \alpha\left\langle g(z) \right\rangle_{p(z)} \)
    • \( \left\langle f(z) + \alpha \right\rangle_{p(z)} = \left\langle f(z) \right\rangle_{p(z)} + \alpha\left\langle 1 \right\rangle_{p(z)} = \left\langle f(z) \right\rangle_{p(z)} + \alpha \)
  • 平均値
    • \( \begin{align} \left\langle z \right\rangle_{p(z)} &= \mu \\ \left\langle z - \mu \right\rangle_{p(z)} &= 0 \end{align} \)
  • 分散
    • \( \left\langle (z - \mu)^2 \right\rangle_{p(z)} = \sigma^2 \)
  • \( x \): 確率変数、\( c \): 定数
    • \( \mathbb E[c] = c \)
    • \( \mathbb E[x + c] = \mathbb E[x] + c \)
    • \( \mathbb E[cx] = c\mathbb E[x] \)
  • \( x, y \): 確率変数
    • \( \mathbb E[x + y] = \mathbb E[x] + \mathbb E[y] \)
\[ \mathbb V[z] = \mathbb E\left[(z - \mu)^2\right] = \sigma^2 \]

マルチモーダル（多峰性）

こちらにあるデータ『sensor_data/sensor_data_600.txt』を使用

モード(mode)が二つあるバイモーダル(bimodal)な分布になっている（山が二つある）（→ 分布の unimodal/multimodal について）

ソースコード

Python
import pandas as pd

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

import matplotlib.pyplot as plt
plt.style.use('ggplot')

plt.clf()
plt.hist(data.lidar, bins=(max(data.lidar) - min(data.lidar) + 1))
plt.title('LiDAR (600mm)'); plt.xlabel('Z [mm]')
plt.savefig('fig5.png')

614, 635 の頻度がゼロになっている 

import pandas as pd

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

values = list(data.lidar.value_counts().index)

for i in range(min(values), max(values) + 1):
    if not i in values:
        print(i)
センサー内部の処理の都合だと考えられるが、ここでは扱わない

時系列でプロット

ソースコード 

import pandas as pd

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

import datetime

data['datetime'] = [ datetime.datetime(int(str(row.date)[0:4]),
                                       int(str(row.date)[4:6]),
                                       int(str(row.date)[6:8]),
                                       hour=row.time//10000,
                                       minute=(row.time//100)%100,
                                       second=row.time%100) for idx, row in data.iterrows() ]

assert data.datetime[0] == data.datetime.min(), 'timestamp error!'
data['timedelta'] = data.datetime - data.datetime[0]

import matplotlib.pyplot as plt
plt.style.use('ggplot')

plt.clf()
plt.plot([ int(s.total_seconds()) for s in data.timedelta.to_list() ], data.lidar.to_list())
plt.title('time series'); plt.xlabel('[s]'); plt.ylabel('LiDAR')
plt.savefig('fig6.png', dpi=125)

散布図

ソースコード

PythonR言語
import pandas as pd

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

data['my_time'] = ((data.time//10000)*60*60
                   + ((data.time//100)%100)*60
                   + data.time%100)

import matplotlib.pyplot as plt
plt.style.use('ggplot')

plt.clf()
for my_day in range(data.date.min(), data.date.max() + 1):
    view = data[data.date == my_day]
    assert len(view) > 0, 'data error!'
    plt.scatter(view.my_time.to_list(), view.lidar.to_list(), alpha=0.5)
plt.title('time vs LiDAR'); plt.xlabel('time'); plt.ylabel('LiDAR')
plt.savefig('fig7.png', dpi=125)
 
data <- read.table('sensor_data_600.txt', sep=' ')
colnames(data) <- c('date', 'time', 'ir', 'lidar')

data$my.time <- (as.integer(data$time/10000)*60*60
    + (as.integer(data$time/100)%%100)*60
    + data$time%%100)

png('fig8.png', width=800, height=600)
plot(data$my.time, data$lidar, type='n',
     main='time vs LiDAR', xlab='time', ylab='LiDAR')
i = 1
pallet <- rainbow(length(min(data$date):max(data$date)))
for (my.day in min(data$date):max(data$date)) {
    view <- data[data$date == my.day, ]
    points(view$my.time, view$lidar, col=pallet[i])
    i <- i + 1
}
dev.off()

時間で条件付けした histogram

６時台と14時台の histogram

（上図：ビンの数はライブラリのデフォルト、下図：ビンの幅を１とした（条件付きの）histogram）

ソースコード

Python
import pandas as pd

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

data['hour'] = data.time//10000

import matplotlib.pyplot as plt
plt.style.use('ggplot')

plt.clf()
for hour in [ 6, 14 ]:
    # view = data[(hour*10000 <= data.time) & (data.time < (hour + 1)*10000)]
    view = data.groupby('hour').get_group(hour)
    plt.hist(view.lidar, alpha=0.6, label=str(hour) + " o'clock")
plt.title('LiDAR (600mm)'); plt.xlabel('Z [mm]')
plt.legend(loc='best')
plt.savefig('fig9.png')

plt.clf()
for hour in [ 6, 14 ]:
    view = data.groupby('hour').get_group(hour)
    plt.hist(view.lidar, alpha=0.6, label=str(hour) + " o'clock",
             bins=(max(view.lidar) - min(view.lidar) + 1))
plt.title('LiDAR (600mm)'); plt.xlabel('Z [mm]')
plt.legend(loc='best')
plt.savefig('fig9_2.png')

条件付き確率

６時台のセンサー値の確率分布の表記の例\[ P(z|t \in \text{6 o'clock}) \]

一般的な条件付き確率の表記

同時確率

「時刻 \( t \) でセンサー値 \( z \) となる確率」。同時確率の表記\[ P(z, t) \]（確率なので \( \sum_z\sum_t P(z, t) = 1 \) となる。）

使用するライブラリのインストール（私の環境では、Linux apt）

Linux apt
$ apt-cache search seaborn
python3-seaborn - statistical visualization library for Python3
$ sudo apt install python3-seaborn

同じデータを heatmap で描画

（※欠損値 614, 635 については無視）

ソースコード

Python
import pandas as pd

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

data['hour'] = data.time//10000

frequency = pd.concat({
    i: data.groupby('hour').lidar.get_group(i).value_counts() for i in range(24)
}, axis=1)
frequency.fillna(0, inplace=True)
density = frequency/len(data)
density.sort_index(ascending=False, inplace=True) # descending

import matplotlib.pyplot as plt
import seaborn as sns

plt.clf()
sns.heatmap(density)
plt.title('P(z, t)'); plt.xlabel('hour'); plt.ylabel('LiDAR')
plt.savefig('fig10.png', dpi=125)

周辺化

（周辺化、周辺分布、周辺確率について）

seaborn の jointplot() を利用

ソースコード
Python
import pandas as pd

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

data['hour'] = data.time//10000

import matplotlib.pyplot as plt
import seaborn as sns

plt.clf()
sns.jointplot(data.hour, data.lidar, data, kind='kde')
plt.savefig('fig11.png', dpi=125)
（jointplot(..., kind='kde') でカーネル密度推定を行っている）

周辺確率の計算

同時確率（の表） \( P(z, t) \) から、周辺確率 \( P(z) \)、\( P(t) \) を（プログラムで）求める（周辺化する）

（※欠損値 614, 635 については無視）

ソースコード

Python
import pandas as pd

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

data['hour'] = data.time//10000

frequency = pd.concat({
    h: data.groupby('hour').lidar.get_group(h).value_counts() for h in range(24)
}, axis=1)
frequency.fillna(0, inplace=True)
density = frequency/len(data)

# joint probability
P_z_t = density # P(z, t)

# marginal probability
P_t = P_z_t.sum() # P(t)
P_z = P_z_t.sum(axis=1) # P(z)
# P_z = P_z_t.transpose().sum()
assert round(P_t.sum(), 7) == 1.0, 'normalize error!'
assert round(P_z.sum(), 7) == 1.0, 'normalize error!'

import matplotlib.pyplot as plt

plt.clf()
# [ plt.figure().get_figwidth(), plt.figure().get_figheight() ]
plt.figure(figsize=(12.8, 4.8))

plt.subplot(121)
P_t.plot()
plt.title('P(t)') # P(t|600mm)
plt.xlabel('time [h]'); plt.ylabel('Density')

plt.subplot(122)
P_z.plot()
plt.title('P(z)') # P(z|600mm)
plt.xlabel('z [mm]'); plt.ylabel('Density')

plt.savefig('fig12.png', dpi=100)

確率の加法定理・乗法定理

確率の加法定理

  • \( P(z) = \sum_{t = -\infty}^{\infty}P(z, t) \)、\( P(t) = \sum_{z = -\infty}^{\infty}P(z, t) \)
  • \( p(z) = \int_{-\infty}^\infty p(z, t) dt \)、\( p(t) = \int_{-\infty}^\infty p(z, t) dz \)

確率の乗法定理

  • \( P(z, t) = P(z|t)P(t) = P(t|z)P(z) \)
  • \( p(z, t) = p(z|t)p(t) = p(t|z)p(z) \)

確率変数が３変数以上の乗法定理

  • \( p(x, y, z) = p(x, y|z)p(z) \)
  • \( p(x, y, z) = p(x|y, z)p(y, z) \)
  • \( p(x, y|z) = p(x|y, z)p(y|z) \)

（３変数以上の場合、確率変数がベクトルになっていると考える）
（式は変わらない）（計算できる）

独立、条件付き独立

（確率変数どうしの関係性）

独立。条件付き確率において、条件 \( y \) が \( x \) の確率に何も影響を与えない場合、次式が成り立つ\[ p(x|y) = p(x) \]つまり、\( y \) が \( x \) に対して何も情報を与えない

乗法定理に上式を代入\[ p(x, y) = p(x|y)p(y) \Longrightarrow p(x, y) = p(x)p(y) \]この関係を \( x, y \) とが互いに独立と表現

（表記。独立の記号が TEX に無い）

条件付き独立。\( z \) が分かっている（\( z \) が観測された）とき、\( x \) に対して \( y \) が何も情報を与えない

（条件付き独立について）《復習》

確率の性質を利用した計算

（ルベグ積分（ルベーグ積分）、無限と集合（無限の集合）、測度論、フビニの定理など）《積み残し》

（確率の性質だけで計算できる）期待値の計算

\( \mathbf z = (x\ y)^\mathrm T \)、\( x \) と \( y \) は独立とする

（和の期待値、期待値の和）\[ \begin{align} \langle f(x) + \alpha g(y) \rangle_{p(\mathbf z)} &= \langle f(x) \rangle_{p(\mathbf z)} + \langle\alpha g(y) \rangle_{p(\mathbf z)} \\ &= \langle f(x) \rangle_{p(\mathbf z)} + \alpha\langle g(y) \rangle_{p(\mathbf z)} \\ &= \langle f(x) \rangle_{p(x)p(y)} + \alpha\langle g(y) \rangle_{p(x)p(y)} \\ &= \langle f(x) \rangle_{p(x)} + \alpha\langle g(y) \rangle_{p(y)} \end{align} \]

（積の期待値、期待値の積）\[ \begin{align} \langle f(x)g(y) \rangle_{p(\mathbf z)} &= \langle f(x)g(y) \rangle_{p(x)p(y)} \\ &= \left\langle\langle f(x)g(y) \rangle_{p(y)}\right\rangle_{p(x)} \\ &= \left\langle f(x)\langle g(y) \rangle_{p(y)}\right\rangle_{p(x)} \\ &= \langle g(y) \rangle_{p(y)}\langle f(x)\rangle_{p(x)} \end{align} \] 得られる関係式\[ \langle g(y) \rangle_{p(y)}\langle f(x) \rangle_{p(x)} = \left\langle \langle f(x)g(y) \rangle_{p(x)} \right\rangle_{p(y)} = \left\langle \langle f(x)g(y) \rangle_{p(y)} \right\rangle_{p(x)} = \langle f(x)g(y) \rangle_{p(x)p(y)} \]

ベイズの定理

乗法定理 \( p(z, t) = p(z|t)p(t) = p(t|z)p(z) \) から導出 \[ \begin{align} p(z|t)p(t) &= p(t|z)p(z) \\ p(z|t) &= \frac{p(t|z)p(z)}{p(t)} = \eta p(t|z)p(z) \end{align} \] \( \eta \) は正規化定数（積分して１にするため）

ベイズの定理を使った計算の例

条件付き確率 \( P(z|t) \) を描画。 \( P(z|t=6) \) と \( P(z|t=14) \)

（※欠損値 614, 635 については無視）

ソースコード

Python
import pandas as pd

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

data['hour'] = data.time//10000

frequency = pd.concat({
    h: data.groupby('hour').lidar.get_group(h).value_counts() for h in range(24)
}, axis=1)
frequency.fillna(0, inplace=True)
density = frequency/len(data)

# joint probability
P_z_t = density # P(z, t)

# marginal probability
P_t = P_z_t.sum() # P(t)

def conditional_probability(cond_t):
    global P_z_t, P_t
    return P_z_t[cond_t]/P_t[cond_t]

import matplotlib.pyplot as plt

plt.clf()
conditional_probability(6).plot.bar(color="blue", alpha=0.5, label="6 o'clock")
conditional_probability(14).plot.bar(color="orange", alpha=0.5, label="14 o'clock")
plt.legend(loc='best')
plt.title('P(z|t)');  plt.xlabel('Z mm'); plt.ylabel('Density')
plt.savefig('fig13.png', dpi=125)

ベイズの定理 \( P(z|t) = {P(t|z)P(z)}/{P(t)} \) が成立していることをプログラムで確認

プログラムと実行結果

Python
import pandas as pd

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

data['hour'] = data.time//10000

frequency = pd.concat({
    h: data.groupby('hour').lidar.get_group(h).value_counts() for h in range(24)
}, axis=1)
frequency.fillna(0, inplace=True)
density = frequency/len(data)

# joint probability
P_z_t = density # P(z, t)

# marginal probability
P_t = P_z_t.sum() # P(t)
P_z = P_z_t.sum(axis=1) # P(z)

def conditional_probability_z(cond_t): # P(z|t)
    global P_z_t, P_t
    return P_z_t[cond_t]/P_t[cond_t]

def conditional_probability_t(cond_z): # P(t|z)
    global P_z_t, P_z
    return P_z_t.loc[cond_z]/P_z[cond_z]

if True: # P.43
    z = 630; t = 13
    print('P(z=%d) ='%z, P_z[z])
    print('P(t=%d) ='%t, P_t[t])
    print('P(t=' + '%d'%t + '|z=' + '%d'%z + ') =', conditional_probability_t(z)[t])
    print('P(z=' + '%d'%z + '|t=' + '%d'%t + ') =', conditional_probability_t(z)[t]*P_z[z]/P_t[t], '// Bayes')
    print('P(z=' + '%d'%z + '|t=' + '%d'%t + ') =', conditional_probability_z(t)[z], '// answer')
    print('')
    
if True:
    error = []
    for z in range(P_z_t.index.min(), P_z_t.index.max() + 1):
        if z in (614, 635): # missing value
            continue
        for t in range(0, 24):
            # Bayes' theorem
            lvalue = conditional_probability_z(t)[z] # P(z|t)
            rvalue = conditional_probability_t(z)[t]*P_z[z]/P_t[t] # P(t|z)P(z)/P(z)
            error.append(abs(lvalue - rvalue))
    print("Bayes' theorem:", round(sum(error), 7) == 0.0)
P(z=630) = 0.06694936878045224
P(t=13) = 0.043024993620976656
P(t=13|z=630) = 0.023230490018148822
P(z=630|t=13) = 0.036147980796385204 // Bayes
P(z=630|t=13) = 0.036147980796385204 // answer

Bayes' theorem: True

ベイズ推定

センサー値が \( z_1 = 630 \)、\( z_2 = 632 \)、\( z_3 = 636 \) と続けて得られたときの \( P(t|\mathbf z) \) を計算

ベイズの定理\[ \begin{align} P(t|z)P(z) &= P(z|t)P(t) \\ P(t|z) &= \frac{P(z|t)P(t)}{P(z)} = \eta P(z|t)P(t) \end{align} \]より \[ P(t|z_1, z_2, z_3)P(z) = \eta P(z_1, z_2, z_3|t)P(t) \]

さらに、条件付き独立より \[ \begin{align} P(t|z_1, z_2, z_3)P(z) &= \eta P(z_1, z_2, z_3|t)P(t) \\ &= \eta P(z_1|t)P(z_2|t)P(z_3|t)P(t) \end{align} \] （隣り合うサンプル値のように \( t \) がほぼ等しい場合は、 i.i.d. を仮定できる）

\( P(z_1|t) \)、\( P(z_2|t) \)、\( P(z_3|t) \) を逐次掛けていくことで、\( t \) を推定することができる（→ ガウス分布のベイズ推定について）

逐次推定結果

\( bel \): prior
observationlikelihood\( bel \): posterior/prior
\[ z_1 = 630 \]
\[ z_2 = 632 \]
observationlikelihood\( bel \): posterior
\[ z_3 = 636 \]

ソースコード

Python
import pandas as pd

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

data['hour'] = data.time//10000

frequency = pd.concat({
    h: data.groupby('hour').lidar.get_group(h).value_counts() for h in range(24)
}, axis=1)
frequency.fillna(0, inplace=True)
density = frequency/len(data)

# joint probability
P_z_t = density # P(z, t)

# marginal probability
P_t = P_z_t.sum() # P(t)

def conditional_probability(cond_t):
    global P_z_t, P_t
    return P_z_t[cond_t]/P_t[cond_t]

import matplotlib.pyplot as plt

plt.clf()
bel = P_t # bel(belief)
bel.plot()
plt.title('prior: ' + r'$ P(t) $'); plt.xlabel('time [h]'); plt.ylabel('Density')
plt.savefig('fig14-0.png')

z = [630, 632, 636] # observation

def bayesian_inference(prior, likelihood):
    posterior = prior*likelihood
    posterior /= posterior.sum() # normalize
    return posterior

# sequential estimation
plt.clf()
bel.plot()
plt.title('prior: ' + r'$ P(t) $' + ', likelihood: ' + r'$ P(z_1|t) $'); plt.xlabel('time [h]'); plt.ylabel('Density')
yyaxis = plt.gca().twinx()
yyaxis.plot([ None ], [ None] ) # adhoc
P_z_1_cond_t = pd.Series([ conditional_probability(t)[z[0]] for t in range(0, 24) ])
P_z_1_cond_t.plot(ax=yyaxis)
plt.savefig('fig14-1.png')

plt.clf()
bel = bayesian_inference(bel, P_z_1_cond_t) # Bayesian inference
bel.plot()
plt.title('posterior: ' + r'$ \eta P(z_1|t)P(t) $'); plt.xlabel('time [h]'); plt.ylabel('Density')
plt.savefig('fig14-2.png')

plt.clf()
bel.plot()
plt.title('prior: ' + r'$ \eta P(z_1|t)P(t) $' + ', likelihood: ' + r'$ P(z_2|t) $'); plt.xlabel('time [h]'); plt.ylabel('Density')
yyaxis = plt.gca().twinx()
yyaxis.plot([ None ], [ None] ) # adhoc
P_z_2_cond_t = pd.Series([ conditional_probability(t)[z[1]] for t in range(0, 24) ])
P_z_2_cond_t.plot(ax=yyaxis)
plt.savefig('fig14-3.png')

plt.clf()
bel = bayesian_inference(bel, P_z_2_cond_t) # Bayesian inference
bel.plot()
plt.title('posterior: ' + r'$ \eta P(z_2|t)P(z_1|t)P(t) $'); plt.xlabel('time [h]'); plt.ylabel('Density')
plt.savefig('fig14-4.png')

plt.clf()
bel.plot()
plt.title('prior: ' + r'$ \eta P(z_2|t)P(z_1|t)P(t) $' + ', likelihood: ' + r'$ P(z_3|t) $'); plt.xlabel('time [h]'); plt.ylabel('Density')
yyaxis = plt.gca().twinx()
yyaxis.plot([ None ], [ None] ) # adhoc
P_z_3_cond_t = pd.Series([ conditional_probability(t)[z[2]] for t in range(0, 24) ])
P_z_3_cond_t.plot(ax=yyaxis)
plt.savefig('fig14-5.png')

# posterior
plt.clf()
bel = bayesian_inference(bel, P_z_3_cond_t) # Bayesian inference
bel.plot()
plt.title('posterior: ' + r'$ \eta P(z_3|t)P(z_2|t)P(z_1|t)P(t) $'); plt.xlabel('time [h]'); plt.ylabel('Density')
plt.savefig('fig14-6.png')

if False: # animation
    import os
    os.system('convert -delay 100 -loop 0 fig14-*.png fig14.gif')

（推定の結果も分布として求まる）

（尤度(likelihood)について。\( P(z_1 | t) \)、\( P(z_2 | t) \)、\( P(z_3 | t) \) が確率の場合には、\( z_1 \)、\( z_2 \)、\( z_3 \) は確率変数。既に \( z_1 \)、\( z_2 \)、\( z_3 \) が観測された状況で、\( t \) の方を変数だと思って \( P(\boldsymbol{z} | t) \) を見た場合には、\( P(z_1 | t) \)、\( P(z_2 | t) \)、\( P(z_3 | t) \) は尤度（尤度関数）となる）

多次元のガウス分布

別の実験データ『sensor_data/sensor_data_700.txt』を可視化して、分布を確認

700[mm]、12～15時台

ソースコード

Python
import pandas as pd

data = pd.read_table('sensor_data_700.txt', delimiter=' ', header=None, names=('date', 'time', 'IR', 'LiDAR'))
data = data[(120000 <= data.time) & (data.time < 160000)]

import matplotlib.pyplot as plt

if False: # P.46
    import seaborn as sns
    plt.clf()
    sns.jointplot(data.IR, data.LiDAR, data, kind='kde')
    plt.savefig('fig15.png')

plt.clf()
plt.axes().set_aspect('equal')
plt.scatter(data.IR.tolist(), data.LiDAR.tolist(), alpha=0.1)
plt.title('IR vs LiDAR'); plt.xlabel('IR'); plt.ylabel('LiDAR')
plt.savefig('fig15_2.png')

（同時分布 \( P(\mathbf z) = P(x, y) \) は（ほぼ丸い）、どちらの確率変数で周辺化してもガウス分布の形をしている）

\( n \) 次元ガウス分布

（多次元のガウス分布は、どの変数を残して周辺化してもガウス分布になる）

\[ \begin{align} \mathcal N(\mathbf z|\boldsymbol\mu, \Sigma) &= \frac{1}{(2\pi)^\frac{n}{2}\sqrt{|\Sigma|}}\exp\left\{-\frac{1}{2}(\mathbf z - \boldsymbol\mu)^\mathrm T\Sigma^{-1}(\mathbf z - \boldsymbol\mu) \right\} \\ &= \eta\exp\left\{ -\frac{1}{2}(\mathbf z - \boldsymbol\mu)^\mathrm T\Sigma^{-1}(\mathbf z - \boldsymbol\mu) \right\} \end{align} \]\( \boldsymbol\mu \)：（平均値、分布の）中心（\( n \) 次元）、\( \Sigma \)：共分散行列（\( n \times n \) 対称行列）（\( n \) 次対称行列）

（\( \eta \) は、積分して１にするための正規化定数）

分布の形は \( e \) の指数乗。\( \exp \)

確率の値を一定にして図を描くと（等高線を引くと）楕円になる（楕円の式とも言える）

（別の形で式を書く）\[ \mathcal N(\mathbf z|\boldsymbol\mu, \Sigma) = \eta e^{ -\frac{1}{2}(\mathbf z - \boldsymbol\mu)^\mathrm T\Sigma^{-1}(\mathbf z - \boldsymbol\mu) } \] １次元のガウス分布との比較

cf.
\[ \mathcal N(z|\mu, \sigma) = \eta e^{ -\frac{(z - \mu)^2}{2\sigma^2} } = \eta e^{ -\frac{1}{2}\frac{(z - \mu)^2}{\sigma^2} } \]

形は、\( e^{ -\frac{1}{2}\frac{(z - \mu)^2}{\sigma^2} } \)。対応を見ると、\( \Sigma^{-1} \) は \( \frac{1}{\sigma^2} \)。スカラーの二乗 \( (z - \mu)^2 \) は、ベクトルになると２次形式(quadratic form) \( (\mathbf z - \boldsymbol\mu)^\mathrm T\Sigma^{-1}(\mathbf z - \boldsymbol\mu) \) の形になる

（１次元を２次元（多次元）に拡張している。式は難しく見えるが、実は変わっていない）

楕円の式（２次形式）

\( a \) を 5、\( b \) を 3 とした楕円

ソースコード

Maxima
load(implicit_plot);

a: 5;
b: 3;

implicit_plot(x^2/a^2 + y^2/b^2 = 1, [x, -6, 6], [y, -5, 5]);

楕円の式を行列で書く\[ \begin{align} \frac{x^2}{a^2} + \frac{y^2}{b^2} &= 1 \\ \begin{pmatrix} x & y \end{pmatrix}\begin{pmatrix} \frac{1}{a^2} & 0 \\ 0 & \frac{1}{b^2} \end{pmatrix}\begin{pmatrix} x \\ y \end{pmatrix} &= 1 \\ \begin{pmatrix} x & y \end{pmatrix}\begin{pmatrix} a^2 & 0 \\ 0 & b^2 \end{pmatrix}^{-1}\begin{pmatrix} x \\ y \end{pmatrix} &= 1 \\ \begin{pmatrix} x \\ y \end{pmatrix}^\mathrm T\begin{pmatrix} a^2 & 0 \\ 0 & b^2 \end{pmatrix}^{-1}\begin{pmatrix} x \\ y \end{pmatrix} &= 1 \end{align} \] 式の形が \( (\mathbf z - \boldsymbol\mu)^\mathrm T\Sigma^{-1}(\mathbf z - \boldsymbol\mu) \) と同じになっている。楕円の形の式だと分かる

式変形の検算

Maxima
(%i1) kill(a, b);
(%o1)                                done
(%i2) (x^2/a^2 + y^2/b^2) - matrix([x, y]).matrix([1/a^2, 0], [0, 1/b^2]).matrix([x], [y]);
(%o2)                                  0

(%i3) (x^2/a^2 + y^2/b^2) - matrix([x, y]).invert(matrix([a^2, 0], [0, b^2])).matrix([x], [y]);
(%o3)                                  0

(%i4) (x^2/a^2 + y^2/b^2) - transpose(matrix([x], [y])).invert(matrix([a^2, 0], [0, b^2])).matrix([x], [y]);
(%o4)                                  0
（Maxima で「.」演算子は行列の積）

共分散行列（分散共分散行列）

２次元の場合

\[ \Sigma = \begin{pmatrix} \sigma_x^2 & \sigma_{xy} \\ \sigma_{yx} & \sigma_y^2 \end{pmatrix} \]\( \sigma_x^2,\ \sigma_y^2 \)：\( x,\ y \) の分散、\( \sigma_{xy},\ \sigma_{yx} \)：共分散

共分散\[ \sigma_{xy} = \sigma_{yx} = \frac{1}{N}\sum_{i = 0}^{N - 1}(x_i - \mu_x)(y_i - \mu_y) \]

（共分散行列は、必ず対称行列になる）

（不偏分散なら \( \frac{1}{N - 1} \) とする。\( \frac{1}{N} \) は、標本分散）（母集団の分散（母分散）を推定しているのか、サンプルの分散（ばらつき）そのものを計算しているのかで使い分ける。\( N \) が大きい場合は誤差で無視できるが、\( N \) が小さい場合には値が変わってくる)

２次元ガウス分布の当てはめ

（パラメータ推定）

センサー値『sensor_data/sensor_data_700.txt』から共分散行列の値を計算

プログラムと実行結果

PythonR言語
import pandas as pd

data = pd.read_table('sensor_data_700.txt', delimiter=' ',
                     header=None, names=['date', 'time', 'IR', 'LiDAR'])
data = data[(120000 <= data.time) & (data.time < 160000)][['IR', 'LiDAR']]

print('variance of IR:', data.IR.var())
print('variance of LiDAR:', data.LiDAR.var())
print('covariance:',
      ((data.IR - data.IR.mean())*(data.LiDAR - data.LiDAR.mean())).sum()/(len(data) - 1))
 
data <- read.table('sensor_data_700.txt', sep=' ',
                   col.names=c('date', 'time', 'IR', 'LiDAR'))
data <- data[(120000 <= data$time) & (data$time < 160000), c('IR', 'LiDAR')]

print(paste('variance of IR:', var(data$IR)))
print(paste('variance of LiDAR:', var(data$LiDAR)))
print(paste('covariance:',
            sum((data$IR - mean(data$IR))*(data$LiDAR - mean(data$LiDAR)))/(dim(data)[1] - 1)))
 
variance of IR: 42.117126367701594
variance of LiDAR: 17.702026469211457
covariance: -0.3167780338543711

[1] "variance of IR: 42.1171263677016"
[1] "variance of LiDAR: 17.7020264692115"
[1] "covariance: -0.316778033854371"

ライブラリ関数による共分散行列の計算

IPythonEmacs Speaks Statistics (ESS)
In [1]: %run covariance_matrix.py
variance of IR: 42.117126367701594
variance of LiDAR: 17.702026469211457
covariance: -0.3167780338543711

In [2]: data.mean()
Out[2]: 
IR        19.860247
LiDAR    729.311958
dtype: float64

In [3]: data.cov()
Out[3]: 
              IR      LiDAR
IR     42.117126  -0.316778
LiDAR  -0.316778  17.702026


 > source('covariance_matrix.R')
 [1] "variance of IR: 42.1171263677016"
 [1] "variance of LiDAR: 17.7020264692115"
 [1] "covariance: -0.316778033854371"
 > c(mean(data$IR), mean(data$LiDAR))
 [1]  19.86025 729.31196
!> var(data)
              IR     LiDAR
 IR    42.117126 -0.316778
 LiDAR -0.316778 17.702026

推定されたパラメタ（平均値と共分散行列）\[ \boldsymbol\mu = \begin{pmatrix} 19.86 \\ 729.312 \end{pmatrix},\ \Sigma = \begin{pmatrix} 42.117 & -0.317 \\ -0.317 & 17.702 \end{pmatrix} \]

\( \boldsymbol\mu \) と \( \Sigma \) の値を使って、ガウス分布を描画

ソースコード

PythonR言語
import pandas as pd

data = pd.read_table('sensor_data_700.txt', delimiter=' ',
                     header=None, names=['date', 'time', 'IR', 'LiDAR'])
data = data[(120000 <= data.time) & (data.time < 160000)][['IR', 'LiDAR']]

mu = data.mean().values
Sigma = data.cov().values

from scipy.stats import multivariate_normal
import numpy as np

x,y = np.mgrid[0:40,710:(710 + 40)]
pos = np.dstack((x, y))

import matplotlib.pyplot as plt

plt.clf()
plt.axes().set_aspect('equal')
plt.contour(x, y,
            multivariate_normal(mean=mu, cov=Sigma).pdf(pos)).clabel(
                # fmt='%1.1e'
            )
plt.title('IR vs LiDAR'); plt.xlabel('IR'); plt.ylabel('LiDAR')
plt.savefig('fig17.png', dpi=125)
 
data <- read.table('sensor_data_700.txt', sep=' ',
                   col.names=c('date', 'time', 'IR', 'LiDAR'))
data <- data[(120000 <= data$time) & (data$time < 160000), c('IR', 'LiDAR')]

mu <- c(mean(data$IR), mean(data$LiDAR))
Sigma <- matrix(var(data), nrow=2)

library(mvtnorm)

x <- seq(0, 40 - 1)
y <- seq(710, (710 + 40) - 1)

png('fig18.png', width=600, height=600)
contour(x, y,
        outer(x, y, function(x, y) dmvnorm(cbind(x, y), mu, Sigma)),
        asp=1,
        main='IR vs LiDAR', xlab='IR', ylab='LiDAR')
dev.off()

多次元のガウス分布には、SciPy の multivariate_normal() 関数と、R言語の mvtnorm パッケージの dmvnorm() 関数を利用

（確率密度関数(probability density function, PDF)）

R言語の mvtnorm パッケージのインストール（私の環境では Linux apt）

Linux apt
$ sudo apt install r-cran-mvtnorm

２次元のカーネル密度推定結果との比較

カーネル密度推定２次元のガウス分布

ソースコード

R言語
data <- read.table('sensor_data_700.txt', sep=' ', col.names=c('date', 'time', 'IR', 'LiDAR'))
data <- data[(120000 <= data$time) & (data$time < 160000), c('IR', 'LiDAR')]

library(MASS)

window <- c(bandwidth.nrd(data$IR), bandwidth.nrd(data$LiDAR))
dense <- kde2d(data$IR, data$LiDAR, window, n=512,
               lims=c(0, 40 - 1, 710, (710 + 40) - 1))

png('fig19.png', width=512, height=512)
image(dense, asp=1, xlab='IR', ylab='LiDAR')
contour(dense, # drawlabels=FALSE,
        add=TRUE)
dev.off()

mu <- c(mean(data$IR), mean(data$LiDAR))
Sigma <- matrix(var(data), nrow=2)

library(mvtnorm)

x <- seq(0, 40 - 1, length=512)
y <- seq(710, (710 + 40) - 1, length=512)

png('fig20.png', width=512, height=512)
image(dense, asp=1, xlab='IR', ylab='LiDAR')
contour(x, y,
        outer(x, y, function(x, y) dmvnorm(cbind(x, y), mu, Sigma)),
        add=TRUE)
dev.off()

２次元のカーネル密度推定には、R言語の MASS パッケージの kde2d() 関数を利用

R言語の MASS パッケージのインストール（私の環境では Linux apt）

Linux apt
$ sudo apt install r-cran-mass

使用したデータの散布図を再掲載

ソースコード
R言語
data <- read.table('sensor_data_700.txt', sep=' ', col.names=c('date', 'time', 'IR', 'LiDAR'))
data <- data[(120000 <= data$time) & (data$time < 160000), c('IR', 'LiDAR')]

png('fig21.png')
plot(data, asp=1, type='n', main='IR vs LiDAR')
points(data, col=rgb(0, 0, 1, alpha=0.1), pch=19, cex=1.6)
dev.off()

共分散行列の意味

それぞれ\[ \mathcal N\left[\boldsymbol\mu = \begin{pmatrix} 50\\ 50 \end{pmatrix}, \Sigma = \begin{pmatrix} 50 & 0 \\ 0 & 100 \end{pmatrix} \right],\ \mathcal N\left[\boldsymbol\mu = \begin{pmatrix} 100\\ 50 \end{pmatrix}, \Sigma = \begin{pmatrix} 125 & 0 \\ 0 & 25 \end{pmatrix} \right],\ \mathcal N\left[\boldsymbol\mu = \begin{pmatrix} 150\\ 50 \end{pmatrix}, \Sigma = \begin{pmatrix} 100 & -25\sqrt{3} \\ -25\sqrt{3} & 50 \end{pmatrix} \right] \]のガウス分布を描画

（共分散行列が、対角行列の場合と一般の場合 → 共分散行列）

ソースコード

Python
import math
import numpy as np
from scipy.stats import multivariate_normal
import matplotlib.pyplot as plt

x,y = np.mgrid[0:200, 0:100]
pos = np.dstack((x, y))

a = {
    'mu': [50, 50],
    'Sigma': [[50, 0],
              [0, 100]]
}
b = {
    'mu': [100, 50],
    'Sigma': [[125, 0],
              [0, 25]]
}
c = {
    'mu': [150, 50],
    'Sigma': [[100, -25*math.sqrt(3)],
              [-25*math.sqrt(3), 50]]
}

plt.clf()
plt.gca().set_aspect('equal')
for parameter in [a, b, c]:
    plt.contour(x, y,
                multivariate_normal(mean=parameter['mu'], cov=parameter['Sigma']).pdf(pos))
plt.xlabel('x'); plt.ylabel('y')
plt.savefig('fig22.png')

分散 \( \sigma_x^2,\ \sigma_y^2 \) が各軸の分布の広がりに対応し、共分散 \( \sigma_{xy} \) に値があると分布が回転することが分かる

ガウス分布の傾き。共分散行列 \( \Sigma \) の固有値と固有ベクトルを求める（→ 固有値と固有ベクトル）

プログラムと実行結果

Python
import math
import numpy as np
from scipy.stats import multivariate_normal

mu = [150, 50]
Sigma = [[100, -25*math.sqrt(3)],
         [-25*math.sqrt(3), 50]]

c = multivariate_normal(mean=mu, cov=Sigma)

val,vec = np.linalg.eig(c.cov)

print('eigenvalue:', val)
print('eigenvector:', vec)
print('eigenvector:', vec[:,0], vec[:,1])
eigenvalue: [125.  25.]
eigenvector: [[ 0.8660254  0.5      ]
 [-0.5        0.8660254]]
eigenvector: [ 0.8660254 -0.5      ] [0.5       0.8660254]

計算結果を図示

ソースコード
Python
import math
import numpy as np
from scipy.stats import multivariate_normal
import matplotlib.pyplot as plt

x,y = np.mgrid[100:200, 0:100]
pos = np.dstack((x, y))

mu = [150, 50]
Sigma = [[100, -25*math.sqrt(3)],
         [-25*math.sqrt(3), 50]]

c = multivariate_normal(mean=mu, cov=Sigma)

val,vec = np.linalg.eig(c.cov)

v = [math.sqrt(val[0])*vec[:,0], math.sqrt(val[1])*vec[:,1]]

plt.clf()
plt.gca().set_aspect('equal')
plt.contour(x, y, c.pdf(pos))
plt.quiver(c.mean[0], c.mean[1], v[0][0], v[0][1], color='red', angles='xy', scale_units='xy', scale=0.5)
plt.quiver(c.mean[0], c.mean[1], v[1][0], v[1][1], color='blue', angles='xy', scale_units='xy', scale=0.5)
plt.xlabel('x'); plt.ylabel('y')
plt.savefig('fig23.png')

共分散行列 \( \Sigma \) と、固有値、固有ベクトルの関係\[ \Sigma = V L V^{-1} \]ここでは、\( V = \begin{pmatrix} \boldsymbol{v}_1 & \boldsymbol{v}_2 \end{pmatrix}\)、\( L = \begin{pmatrix} l_1 & 0 \\ 0 & l_2 \end{pmatrix} \)。計算結果\[ \begin{pmatrix} 100 & -25\sqrt{3} \\ -25\sqrt{3} & 50 \end{pmatrix} = \begin{pmatrix} 0.8660254 & 0.5 \\ -0.5 & 0.8660254 \end{pmatrix} \begin{pmatrix} 125.0 & 0 \\ 0 & 25.0 \end{pmatrix} \begin{pmatrix} 0.8660254 & 0.5 \\ -0.5 & 0.8660254 \end{pmatrix}^{-1} \]

数値計算による確認。プログラムと実行結果

Python
import math
import numpy as np
from scipy.stats import multivariate_normal

mu = [150, 50]
Sigma = [[100, -25*math.sqrt(3)],
         [-25*math.sqrt(3), 50]]

c = multivariate_normal(mean=mu, cov=Sigma)

val,vec = np.linalg.eig(c.cov)

if False:
    print('v_1:', vec[:,0])
    print('v_2:', vec[:,1])
    print('l_1:', val[0])
    print('l_2:', val[1])

print((c.cov - vec@np.diag(val)@np.linalg.inv(vec)).round(7))
[[ 0.  0.]
 [-0. -0.]]
（Python では、Python 3.5 から行列の積のために「@」演算子が提供されている）

線形代数の復習

対称行列の対角化
diagonalization
  • \( n\times n \) の行列 \( A \) が実対称行列であれば、重解も含めて \( n \) 個の固有値が存在する。
  • 行列 \( A \) の固有値を \( \lambda_1, \lambda_2, \cdots, \lambda_n \) とし、対応する固有ベクトルを \( \boldsymbol{u}_1, \boldsymbol{u}_2, \cdots, \boldsymbol{u}_n \) とする。
  • それぞれの固有ベクトルは、大きさが \( 1 \) となるように正規化されている。\( \boldsymbol{u}_i = 1 \)
  • 実対称行列の固有ベクトルは互いに直交する。（一次独立、線形独立）\[ \boldsymbol{u}_i^\mathrm T\boldsymbol{u}_j = \begin{cases} 0 & (i \neq j)\\ 1 & (i = j) \end{cases} \]
  • \( n \) 次元の固有ベクトルを \( n \) 個並べた行列 \( U \) を作ると、\( U \) は直交行列となる。\[ U = \begin{bmatrix} \boldsymbol{u}_1, \boldsymbol{u}_2, \cdots, \boldsymbol{u}_n \end{bmatrix},\ U^\mathrm T = U^{-1} \]
  • 行列 \( U \) を用いて、行列 \( A \) を対角化できる。\[ U^\mathrm T A U = \Lambda,\ \Lambda = \begin{bmatrix} \lambda_1 & & & 0 \\ & \lambda_2 & & \\ & & \ddots & \\ 0 & & & \lambda_n \end{bmatrix} \]
  • \( U \) は直交行列。\( U^\mathrm T U = U^{-1} U = I \)\[ \begin{align} U^\mathrm T A U &= \Lambda \\ U(U^\mathrm T A U) &= U\Lambda \\ U(U^\mathrm T A U)U^\mathrm T &= U\Lambda U^\mathrm T \\ (U U^\mathrm T)A(U U^\mathrm T) &= U\Lambda U^\mathrm T \\ A &= U\Lambda U^\mathrm T \end{align} \]（行列の演算は「結合法則」が成り立つ）（行列 \( A \) を２つの直交行列と１つの対角行列の積に分解(decomposition)できる）
相似変換
similarity transformation
  • 正則な行列 \( P \) を使って、行列 \( A \) に対して\[ P^{-1}A P \]のような変換を行う。このような変換を「相似変換」という。
  • 相似変換をしても、「行列式」の値は変化しない。\[ |P^{-1}A P| = |P^{-1}|\cdot|A|\cdot|P| = \frac{1}{|P|}\cdot|A|\cdot|P| = |A| \]
  • 相似変換に対して「固有値」も変化しない。「固有多項式」が変化しないことから説明。\[ |P^{-1}A P - \lambda I| = |P^{-1}A P - P^{-1}\lambda P| = |P^{-1}(A - \lambda I)P| = |P^{-1}|\cdot |A - \lambda I|\cdot |P| = \frac{1}{|P|}\cdot |A - \lambda I|\cdot |P| = |A - \lambda I| \]
  • 相似変換後の行列 \( P^{-1}A P \) は、もとの行列 \( A \) と「よく似た性質を持つ」。（行列 \( A \) の性質（行列式、固有値）を引き継いでいる）

共分散行列は対角行列を回転行列（と回転行列の逆行列）で挟んだ形に（分解）できる《証明は教科書》

２次元の回転行列\[ R_\theta = \begin{pmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{pmatrix} \]

計算結果を回転行列で表現\[ \Sigma = \begin{pmatrix} 0.8660254 & 0.5 \\ -0.5 & 0.8660254 \end{pmatrix} \begin{pmatrix} 125.0 & 0 \\ 0 & 25.0 \end{pmatrix} \begin{pmatrix} 0.8660254 & 0.5 \\ -0.5 & 0.8660254 \end{pmatrix}^{-1} = R_{-\pi/6} \begin{pmatrix} 125.0 & 0 \\ 0 & 25.0 \end{pmatrix} R_{-\pi/6}^{-1} \]

数値計算による確認。プログラムと実行結果

Python
import math
import numpy as np
from scipy.stats import multivariate_normal

mu = [150, 50]
Sigma = [[100, -25*math.sqrt(3)],
         [-25*math.sqrt(3), 50]]

c = multivariate_normal(mean=mu, cov=Sigma)
val,vec = np.linalg.eig(c.cov)

theta = -np.pi/6
R = np.array([[np.cos(theta), -np.sin(theta)],
              [np.sin(theta), np.cos(theta)]])

if False:
    print(vec); print(R)

print((vec - R).round(7))
[[ 0.  0.]
 [-0.  0.]]

回転行列によるガウス分布の回転

ソースコード
Python
import math
import numpy as np
from scipy.stats import multivariate_normal
import matplotlib.pyplot as plt

x,y = np.mgrid[100:200, 0:100]
pos = np.dstack((x, y))

mu = [150, 50]
Sigma = np.array([[125, 0],
                  [0, 25]])

for angle in range(0, 360, 15):
    theta = -(angle*(np.pi/180))
    R = np.array([[np.cos(theta), -np.sin(theta)],
                  [np.sin(theta), np.cos(theta)]])
    c = multivariate_normal(mean=mu, cov=R@Sigma@R.T)
    val,vec = np.linalg.eig(c.cov)
    v = [math.sqrt(val[0])*vec[:,0], math.sqrt(val[1])*vec[:,1]]
    plt.clf()
    plt.gca().set_aspect('equal')
    plt.contour(x, y, c.pdf(pos))
    palet = ['red' if np.linalg.norm(v[0]) > np.linalg.norm(v[1]) else 'blue']
    palet.append('blue' if palet[0] == 'red' else 'red')
    plt.quiver(c.mean[0], c.mean[1], v[0][0], v[0][1], color=palet[0], angles='xy', scale_units='xy', scale=0.5)
    plt.quiver(c.mean[0], c.mean[1], v[1][0], v[1][1], color=palet[1], angles='xy', scale_units='xy', scale=0.5)
    plt.title(('-' if angle != 0 else '') + str(angle) + ' [degree]'); plt.xlabel('x'); plt.ylabel('y')
    plt.savefig('fig24-' + '%03d'%angle + '.png', dpi=125)

if False: # animation
    import os
    os.system('convert -delay 100 -loop 0 fig24-*.png fig24.gif')
（回転行列の性質（直交行列）\( R_\theta^{-1} = R_\theta^\mathrm T \)）

共分散行列と誤差楕円

ガウス分布を楕円で描く。たとえば、長軸、短軸を標準偏差を３倍して楕円を描いたもの。標準偏差で正規化されているので、楕円の大きさで分布を比較することができる

３倍の意味。長軸、短軸方向とも \( 3\sigma \)。積分すると

になる（確率が決まっている。図示できる）

ガウス分布の合成

確率変数がガウス分布に従う場合

確率変数の和

ガウス分布同士の和

（式の導出。１次元の場合）《導出は教科書》

シミュレーション

ソースコード
Python
import numpy as np
from scipy.stats import multivariate_normal
import matplotlib.pyplot as plt

x,y = np.mgrid[0:200, 0:100]
pos = np.dstack((x, y))

count = 0
if True:
    bel = {
        'mu': np.array([50, 50]),
        'Sigma': 25*np.eye(2)
    }
    u = {
        'mu': np.array([10, 0]),
        'Sigma': 25*np.eye(2) # isotropic
    }
    for count in range(10):
        plt.clf(); plt.gca().set_aspect('equal')
        plt.contour(x, y, multivariate_normal(mean=bel['mu'], cov=bel['Sigma']).pdf(pos))
        plt.title(r'$ \Sigma $' + ': isotropic'); plt.xlabel('x'); plt.ylabel('y')
        plt.savefig('fig25-' + '%02d'%count + '.png', dpi=125)
        bel['mu'] += u['mu']; bel['Sigma'] += u['Sigma'] # update bel

if True:
    bel = {
        'mu': np.array([50, 50]),
        'Sigma': 25*np.eye(2)
    }
    u = {
        'mu': np.array([10, 0]),
        'Sigma': np.diag([5, 25]) # diagonal
    }
    for count in range(count, count + 10):
        plt.clf(); plt.gca().set_aspect('equal')
        plt.contour(x, y, multivariate_normal(mean=bel['mu'], cov=bel['Sigma']).pdf(pos))
        plt.title(r'$ \Sigma $' + ': diagonal'); plt.xlabel('x'); plt.ylabel('y')
        plt.savefig('fig25-' + '%02d'%count + '.png', dpi=125)
        bel['mu'] += u['mu']; bel['Sigma'] += u['Sigma'] # update bel

if True:
    bel = {
        'mu': np.array([50, 50]),
        'Sigma': 25*np.eye(2)
    }
    theta = -60*(np.pi/180)
    R = np.array([[np.cos(theta), -np.sin(theta)],
                  [np.sin(theta), np.cos(theta)]]) 
    u = {
        'mu': np.array([10, 0]),
        'Sigma': R@np.diag([5, 25])@R.T # general form
    }
    for count in range(count, count + 10):
        plt.clf(); plt.gca().set_aspect('equal')
        plt.contour(x, y, multivariate_normal(mean=bel['mu'], cov=bel['Sigma']).pdf(pos))
        plt.title(r'$ \Sigma $' + ': general form'); plt.xlabel('x'); plt.ylabel('y')
        plt.savefig('fig25-' + '%02d'%count + '.png', dpi=125)
        bel['mu'] += u['mu']; bel['Sigma'] += u['Sigma'] # update bel
        
if False: # animation
    import os
    os.system('convert -delay 100 -loop 0 fig25-*.png fig25.gif')

（復習）

状態空間、制御出力

\[ \boldsymbol{x} \in \mathcal X,\ \boldsymbol{u} \in \mathcal U \]（状態変数で表現されたものしか扱えない）（フレーム問題）

離散時間を考える

\[ t \in \mathcal T \]（状態遷移には時間が掛かる）

（雑音を考慮した）状態方程式

線形な状態方程式\[ \boldsymbol{x}_t = A\boldsymbol{x}_{t-1} + B\boldsymbol{u}_t + \boldsymbol{\epsilon}_t \]非線形な状態方程式\[ \boldsymbol{x}_t = f(\boldsymbol{x}_{t-1}, \boldsymbol{u}_t) + \boldsymbol{\epsilon}_t \]（\( f() \) は状態遷移関数）

状態遷移の確率密度関数(PDF)

状態遷移を、確率の式として書く\[ \boldsymbol{x}_t \sim p(\boldsymbol{x} | \boldsymbol{x}_{t-1}, \boldsymbol{u}_t) \]

状態に関する確率密度関数(PDF)

\[ bel_t(\boldsymbol{x}) = \int_\mathcal{X} p(\boldsymbol{x} | \boldsymbol{x}’, \boldsymbol{u}_t) bel_{t-1}(\boldsymbol{x}’) d\boldsymbol{x}’ \]（\( \boldsymbol{x}’ \) は \( \boldsymbol{x}_{t-1} \) のある１サンプル点）（\( bel \) は belief の略。信念）

積分計算の意味を、シミュレーションで理解

ソースコード

Python
import numpy as np
from scipy.stats import multivariate_normal
import matplotlib.pyplot as plt

x,y = np.mgrid[0:200, 0:100]
pos = np.dstack((x, y))

bel = {
    'mu': np.array([50, 50]),
    'Sigma': 25*np.eye(2) # isotropic
}
u = {
    'mu': np.array([20, 0]),
    'Sigma': np.diag([5, 25]) # diagonal
}

np.random.seed(0)

frame = 0
for t in range(0, 5 + 1):
    plt.clf(); plt.gca().set_aspect('equal')
    state = multivariate_normal(mean=bel['mu'], cov=bel['Sigma'])
    plt.contour(x, y, state.pdf(pos))
    plt.title('t = ' + str(t)); plt.xlabel('x'); plt.ylabel('y')
    plt.savefig('fig25_2-' + '%02d'%frame + '.png', dpi=125); frame += 1
    samples = state.rvs(size=32) # draw
    plt.scatter(samples[:,0].tolist(), samples[:,1].tolist())
    plt.savefig('fig25_2-' + '%02d'%frame + '.png', dpi=125); frame += 1
    for sample in samples:
        plt.plot([sample[0], (sample + u['mu'])[0]], [sample[1], (sample + u['mu'])[1]], ls='--')
        plt.contour(x, y, multivariate_normal(mean=(sample + u['mu']), cov=u['Sigma']).pdf(pos), levels=3)
    plt.savefig('fig25_2-' + '%02d'%frame + '.png', dpi=125); frame += 1
    bel['mu'] += u['mu']; bel['Sigma'] += u['Sigma'] # update bel
    plt.contour(x, y, multivariate_normal(mean=bel['mu'], cov=bel['Sigma']).pdf(pos))
    plt.savefig('fig25_2-' + '%02d'%frame + '.png', dpi=125); frame += 1

if False: # animation
    import os
    os.system('convert -delay 100 -loop 0 fig25_2-*.png fig25_2.gif')
R言語
library(mvtnorm)

bel.mu <- c(50, 50)
bel.Sigma <- 25*diag(rep(1, 2)) # isotropic

u.mu <- c(20, 0)
u.Sigma <- diag(c(5, 25)) # diagonal

x <- seq(0, 200 - 1)
y <- seq(0, 100 - 1)

set.seed(0)

initial.belief <- function(t) {
    contour(x, y,
            outer(x, y, function(x, y) dmvnorm(cbind(x, y), bel.mu, bel.Sigma)),
            drawlabels=FALSE, asp=1,
            main=paste('t =', t), xlab='x', ylab='y')
}

x.dash <- function(samples) {
    points(samples)
}

motion <- function(samples) {
    for (i in seq(dim(samples)[1])) {
        my.sample <- samples[i,]
        lines(c(my.sample[1], (my.sample + u.mu)[1]), c(my.sample[2], (my.sample + u.mu)[2]), lty='dashed')
        contour(x, y,
                outer(x, y, function(x, y) dmvnorm(cbind(x, y), my.sample + u.mu, u.Sigma)), nlevels=3,
                drawlabels=FALSE, asp=1,
                add=TRUE)
    }
}

belief.after.motion <- function() {
    contour(x, y,
            outer(x, y, function(x, y) dmvnorm(cbind(x, y), bel.mu, bel.Sigma)),
            drawlabels=FALSE, asp=1,
            add=TRUE)
}

my.frame <- 0
get.fname <- function() {
    fname <- paste(paste('fig25_3', formatC(my.frame, width=2, flag='0'), sep='-'), 'png', sep='.')
    my.frame <<- my.frame + 1
    return(fname)
}

for (t in 0:5) {
    png(get.fname(), width=800, height=600)
    initial.belief(t)
    dev.off()

    samples <- rmvnorm(32, bel.mu, bel.Sigma) # draw

    png(get.fname(), width=800, height=600)
    initial.belief(t)
    x.dash(samples)
    dev.off()

    png(get.fname(), width=800, height=600)
    initial.belief(t)
    x.dash(samples)
    motion(samples)
    dev.off()

    png(get.fname(), width=800, height=600)
    initial.belief(t)
    x.dash(samples)
    motion(samples)
    bel.mu <- bel.mu + u.mu; bel.Sigma <- bel.Sigma + u.Sigma # update bel
    belief.after.motion()
    dev.off()
}

if (FALSE) # animation
    system("convert -delay 100 -loop 0 fig25_3-*.png fig25_3.gif")

１次元の問題の例

状態は１次元 \( x \in \Re \)。制御出力 \( u \) は状態の変位 \( \Delta_t \)（制御出力も状態と同じ次元にある一般的ではない場合）\[ x_t = x_{t - 1} + \Delta_t + \epsilon \]雑音 \( \epsilon \) は標準偏差 \( \delta_t \) でガウス分布に従う

この \( bel(x_t) \) の一般解を求める\[ bel(x_t) = \int_{\mathcal X} p(x_t | x_{t-1}, \Delta_t) bel(x_{t-1}) dx_{t-1} \]

以下から導出\[ p(x | x_{t-1}, \Delta_t) = \frac{1}{\sqrt{2\pi\delta^2_t}}\exp\left\{ -\frac{(x - x_{t-1} - \Delta_t)^2}{2\delta^2_t} \right\},\ bel(x_{t-1}) = \frac{1}{\sqrt{2\pi\sigma^2_{t-1}}}\exp\left\{ -\frac{(x_{t-1} - \hat{x}_{t-1})^2}{2\sigma^2_{t-1}} \right\} \]\( \hat{x}_{t-1} \) は \( t - 1 \) における分布の中心（状態の推定値）、\( \sigma_{t-1} \) は \( t - 1 \) における分布の標準偏差

\( bel \) の式にガウス分布を代入\[ bel(x_t) = \int_{-\infty}^{\infty} \frac{1}{\sqrt{2\pi\delta^2_t}}\exp\left\{ -\frac{(x - x_{t-1} - \Delta_t)^2}{2\delta^2_t} \right\} \frac{1}{\sqrt{2\pi\sigma^2_{t-1}}}\exp\left\{ -\frac{(x_{t-1} - \hat{x}_{t-1})^2}{2\sigma^2_{t-1}} \right\} dx_{t-1} \]

Maxima で記号計算

Maxima
(%i1) assume(sigma[t-1] > 0);
(%o1)                          [sigma      > 0]
                                     t - 1
(%i2) assume(delta[t] > 0);
(%o2)                            [delta  > 0]
                                       t
(%i3) tex(integrate(1/sqrt(2*pi*delta[t]^2)*exp(-1/(2*delta[t]^2)*(x - x[t-1] - Delta[t])^2)*1/sqrt(2*pi*sigma[t-1]^2)*exp(-1/(2*sigma[t-1]^2)*(x[t-1] - \\hat\{x\}[t-1])^2), x[t-1], minf, inf));
\[{{\sqrt{\pi}\,e^ {- {{x^2+\left(-2\,\Delta_{t}-2\,{\it \hat{x}}_{t- 1}\right)\,x+\Delta_{t}^2+2\,{\it \hat{x}}_{t-1}\,\Delta_{t}+ {\it \hat{x}}_{t-1}^2}\over{2\,\delta_{t}^2+2\,\sigma_{t-1}^2}} } }\over{\sqrt{2}\,\pi\,\sqrt{\delta_{t}^2+\sigma_{t-1}^2}}}\]

正規化定数 \( \eta \) を導入して式を整理\[\eta\,e^ {- {{x^2+\left(-2\,\Delta_{t}-2\,{\it \hat{x}}_{t- 1}\right)\,x+\Delta_{t}^2+2\,{\it \hat{x}}_{t-1}\,\Delta_{t}+ {\it \hat{x}}_{t-1}^2}\over{2\,\delta_{t}^2+2\,\sigma_{t-1}^2}} }\]（※Maxima の計算結果の TEX を手作業で修正）

Maxima の factor() 関数を使って \( e \) の肩の式の形を（ガウス分布として）見やすいように整理

Maxima
(%i1) tex(eta*exp(-factor(x^2 + (-2*Delta[t] - 2*\\hat\{x\}[t-1])*x + Delta[t]^2 + 2*\\hat\{x\}[t-1]*Delta[t] + \\hat\{x\}[t-1]^2)/factor(2*delta[t]^2 + 2*sigma[t-1]^2)));
\[\eta\,e^ {- {{\left(x-\Delta_{t}-{\it \hat{x}}_{t-1}\right)^2 }\over{2\,\left(\delta_{t}^2+\sigma_{t-1}^2\right)}} }\]

結局 \( bel(x_t) \) は\[ bel(x_t) = \eta\exp\left\{ -\frac{(x - \hat{x}_{x-1} - \Delta_t)^2}{2(\sigma^2_{t-1} + \delta^2_t)} \right\} = \mathcal{N}(\hat{x}_{t-1} + \Delta_t, \sigma^2_{t-1} + \delta^2_t) \]再びガウス分布となる

導出された結論の式の意味の理解

  • 分布の中心：もとの移動量に制御出力の値を加えた値
  • 分散：もとの（\( bel(x_{t-1}) \) の）分散に移動による雑音の分散を足した値

シミュレーション

ソースコード

Python
import math
import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt

bel = {
    'mu': 50,
    'sigma': math.sqrt(20)
}
Delta_t = 10; delta_t = math.sqrt(25)
u = {
    'mu': Delta_t,
    'sigma': delta_t
}

x = np.linspace(0, 200 - 1, 800)
plt.clf()
for t in range(10):
    y = norm.pdf(x, bel['mu'], bel['sigma'])
    plt.plot(x, y, label=r'$ bel(x_' + str(t) + ') $', alpha=0.8)
    bel['mu'] += u['mu']; bel['sigma'] = math.sqrt(bel['sigma']**2 + u['sigma']**2) # motion update
plt.legend(loc='best')
plt.title('t = ' + str(list(range(10)))); plt.xlabel('x'); plt.ylabel('Density')
plt.savefig('fig25_4.png', dpi=125)

（カルマンフィルターの移動更新 (motion update) の式）（1次元で、いろいろ省略されている）

パーティクルフィルター

（ガウス分布に従わない場合、数値計算で解く）（モンテカルロ法）

ソースコード

Python
import math
import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt

bel_0 = {
    'mu': 50,
    'sigma': math.sqrt(20)
}
Delta_t = 10; delta_t = math.sqrt(25)
u = {
    'mu': Delta_t,
    'sigma': delta_t
}

t = 0 # initial belief
plt.clf()
x = np.linspace(0, 200 - 1, 800)
y = norm.pdf(x, bel_0['mu'], bel_0['sigma'])
plt.plot(x, y)
plt.title(r'$ bel(x_' + str(t) + ') $'); plt.xlabel('x'); plt.ylabel('Density')
frame = 0
plt.savefig('fig25_5-' + '%02d'%frame + '.png', dpi=125); frame += 1
particle = np.random.normal(bel_0['mu'], bel_0['sigma'], 4096) # draw(sampling)
plt.plot(particle, len(particle)*[0], marker='.', linestyle='None', alpha=0.1)
plt.savefig('fig25_5-' + '%02d'%frame + '.png', dpi=125); frame += 1
my_xlim = plt.gca().get_xlim(); my_ylim = plt.gca().get_ylim()
t += 1

for t in range(t, 10):
    plt.clf(); plt.xlim(my_xlim); plt.ylim(my_ylim)
    plt.plot(x, y)
    particle += np.random.normal(u['mu'], u['sigma'], len(particle)) # motion update
    plt.plot(particle, len(particle)*[0], marker='.', linestyle='None', alpha=0.1)
    plt.title(r'$ bel(x_' + str(t) + ') $'); plt.xlabel('x'); plt.ylabel('Density')
    plt.savefig('fig25_5-' + '%02d'%frame + '.png', dpi=125); frame += 1
    plt.clf(); plt.xlim(my_xlim); plt.ylim(my_ylim)
    y = norm.pdf(x, particle.mean(), particle.std())
    plt.plot(x, y)
    plt.plot(particle, len(particle)*[0], marker='.', linestyle='None', alpha=0.1)
    plt.title(r'$ bel(x_' + str(t) + ') $'); plt.xlabel('x'); plt.ylabel('Density')
    plt.savefig('fig25_5-' + '%02d'%frame + '.png', dpi=125); frame += 1

if False:
    import os
    os.system('convert -delay 100 -loop 0 fig25_5-*.png fig25_5.gif')

（自己位置推定に使われている、モンテカルロローカライゼーションの移動の式）（乱数を使っている。パーティクル（点）の数が少ないと劣化する）

（他にも（格子（格子地図）を使った \( bel \) の計算）「ヒストグラムフィルター」など）

参照資料

千葉工業大学の上田隆一先生のスライド『確率ロボティクス2016第1回、第2回』を拝読させていただきました。m(_o_)m

確率変数の積

ガウス分布同士の積

精度行列 \( \Lambda \) を導入（共分散行列 \( \Sigma \) の逆行列）

精度行列を使って式を書き直す

（式の導出。１次元の場合で、新たな観測の数 \( N \) が一般の場合、と \( N = 1 \) の場合に限定して式変形したもの）《復習》

結論の式の意味の理解。精度は足し算になる。平均値は、精度の重み付き平均になっている（結果は、精度は必ず上がり（分散が小さくなる）、平均値（中心）は精度の良い分布の方に寄る）（ガウス分布同士の掛け算は、結局、ガウス分布になる）

シミュレーション（→ １次元のベイズ推定）

ソースコード
Python
import numpy as np
from scipy.stats import multivariate_normal
import matplotlib.pyplot as plt

x,y = np.mgrid[0:200, 0:100]
pos = np.dstack((x, y))

bel = {
    'mu': np.array([50, 50]),
    'Sigma': 25*np.eye(2)
}
plt.clf(); plt.gca().set_aspect('equal')
plt.contour(x, y, multivariate_normal(mean=bel['mu'], cov=bel['Sigma']).pdf(pos))
plt.title('prior'); plt.xlabel('x'); plt.ylabel('y')
plt.savefig('fig26-0.png', dpi=125)

z = { # observation, likelihood
    'mu': np.array([150, 50]),
    'Sigma': np.diag([50, 125]) # diagonal
}
plt.clf(); plt.gca().set_aspect('equal')
plt.contour(x, y, multivariate_normal(mean=bel['mu'], cov=bel['Sigma']).pdf(pos))
plt.contour(x, y, multivariate_normal(mean=z['mu'], cov=z['Sigma']).pdf(pos))
plt.title('prior ' + r'$ \times $' + ' observation, likelihood'); plt.xlabel('x'); plt.ylabel('y')
plt.savefig('fig26-1.png', dpi=125)

def update(bel, z):
    Sigma = np.linalg.inv(np.linalg.inv(bel['Sigma']) + np.linalg.inv(z['Sigma']))
    mu = Sigma@(np.linalg.inv(bel['Sigma'])@bel['mu'] + np.linalg.inv(z['Sigma'])@z['mu'])
    return { 'mu': mu, 'Sigma': Sigma }

bel = update(bel, z) # Bayesian inference
    
plt.clf(); plt.gca().set_aspect('equal')
plt.contour(x, y, multivariate_normal(mean=bel['mu'], cov=bel['Sigma']).pdf(pos))
plt.title('posterior'); plt.xlabel('x'); plt.ylabel('y')
plt.savefig('fig26-2.png', dpi=125)

if False: # animation
    import os
    os.system('convert -delay 100 -loop 0 fig26-*.png fig26.gif')

実験データの可視化

（こちらの『sensor_data/』ディレクトリにあるファイルを使わせていただきました）

（教科書「問 2.8」）

最小二乗法による最尤推定

（確率ロボティクスは、観測変数が \( Z \)、潜在変数が \( X \)）

破線が最尤推定解(maximum likelihood)

モデル選択（多項式の次数）は手作業による試行錯誤。IR センサーは３次式、LiDAR センサーは１次式に決定。また、IR センサーは（近い距離の実験データが足りなく）上手くモデルが作れなかったため、定義域 \( ( 0 \leq z \leq 75 ) \) を適当に設定

ソースコード

Python
import glob
import re
import pandas as pd

data_files = glob.glob('sensor_data_*.txt')

count = 0
for data_file in data_files:
    m = re.match(r'sensor_data_(\d{3,4})\.txt', data_file)
    assert m, 'data file name error!'
    distance = m.groups()
    assert len(distance) == 1 and (len(distance[0]) == 3 or len(distance[0]) == 4), 'data file name error!'
    distance = int(distance[0])
    temp = pd.read_table(data_file, delimiter=' ', header=None, usecols=[2, 3], names=['ir', 'lidar'])
    if count == 0:
        data = temp.assign(distance=distance)
        distances = [distance]
    else:
        data = data.append(temp.assign(distance=distance), ignore_index=True)
        distances.append(distance)
    count += 1

import matplotlib.pyplot as plt
plt.style.use('ggplot')

fig1 = plt.figure(1); fig1.clf()
fig2 = plt.figure(2); fig2.clf()

distances.sort()
for distance in distances:
    view = data[data.distance == distance]
    plt.figure(1)
    plt.scatter(view.ir, len(view)*[distance], alpha=0.1)
    plt.figure(2)
    plt.scatter(view.lidar, len(view)*[distance], alpha=0.1)
plt.figure(1)
plt.title('IR'); plt.xlabel('z'); plt.ylabel('x [mm]')
plt.figure(2)
plt.title('LiDAR'); plt.xlabel('z'); plt.ylabel('x [mm]')

import numpy as np

model_comparison = 3 # heuristics
view = data[data.distance >= 300] # pre-process
w = np.polyfit(view.ir.tolist(), view.distance.tolist(), model_comparison)
assert len(w) == model_comparison + 1, 'model error!'
poly_ir = np.poly1d(w)

x = np.linspace(0,
                75, # adhoc
                800)

plt.figure(1)
plt.plot(x, poly_ir(x), label='[' + ', '.join([str(round(coef, 4)) for coef in w]) + ']', ls='--', c='k')
plt.legend(loc='best')

model_comparison = 1 # heuristics

w = np.polyfit(data.lidar.tolist(), data.distance.tolist(), model_comparison)
assert len(w) == 2, 'model error!'
poly_lidar = np.poly1d(w)

x = np.linspace(200, 1000, 800)

plt.figure(2)
plt.plot(x, poly_lidar(x), label=r'$' + str(round(w[0], 4)) + 'z + ' + str(round(w[1], 4)) + r'$', ls='--', c='k')
plt.legend(loc='best')

fig1.savefig('fig27.png', dpi=125)
fig2.savefig('fig28.png', dpi=125)

（NumPy の polyfit() 関数と poly1d() 関数を利用）

推定結果（\( x_{\mathrm{ML}} \)）

IPython
In [2]: print(poly_ir)                                                                                      
            3          2
-0.0007894 x + 0.2244 x - 21.93 x + 1052

In [3]: print(poly_lidar)                                                                                   
 
0.9517 x + 3.488

（過学習、ベイズ線形回帰、予測分布について）《復習》

教科書

© 2020 Tadakazu Nagai