詳解　確率ロボティクス　Ｐｙｔｈｏｎによる基礎アルゴリズムの実装 (ＫＳ理工学専門書)

作者:上田隆一
発売日: 2019/12/20
メディア: Kindle版

はじめに

前回の記事では、シンプルな正規分布に
従うセンサ値のモデル化とそれを計算する
Juliaサンプルコードについて紹介しました。

www.eureka-moments-blog.com

しかしながら実際は、こういった綺麗な
分布でモデル化できない複雑な現象も
沢山あります。

今回の記事では、複雑な分布を持つ
センサ値に対する前処理とモデル化、
ベイズの定理による原因推定を
行うJuliaサンプルコードについて紹介
します。

GitHubリポジトリ

本記事で紹介するJuliaコードは全てこちらの
リポジトリにて公開しています。

github.com

複雑な分布の例

まずは例として、JuliaAutonomy/dataディレクトリ
にあるサンプルデータの中から、sensor_data_600.txtを
読み込んで可視化してみましょう。

こちらのJuliaコードを実行すると、データを
読み込んでヒストグラムを描画してくれます。
src/prob_sats/complex_dist/histogram_600mm/draw_hist_600.jl

module DrawHist600
    using Plots, DataFrames, CSV
    pyplot()

    function main(is_test=false)
        data_path = joinpath(split(@__FILE__, "src")[1], "data/sensor_data_600.txt")
        df_600_mm = CSV.read(data_path, DataFrame, 
                             header=["date", "time", "ir", "lidar"],
                             delim=' ')
        
        bin_min_max = maximum(df_600_mm.lidar) - minimum(df_600_mm.lidar)
        histogram(df_600_mm.lidar, bins=bin_min_max, label="histogram")

        if is_test == false
            save_path = joinpath(split(@__FILE__, "src")[1], "src/prob_stats/complex_dist/histogram_600mm/histogram_600_mm.png")
            savefig(save_path)
        end
    end
end

そしてこちらが描画されるヒストグラムです。
f:id:sy4310:20210511221022p:plain

前回の記事で見たような正規分布と違い、
出現頻度のピークが二つあることが
分かります。
こういった分布は「マルチモーダルな(多峰性の)
分布」と呼ばれます。

ここから先は、このデータに対するモデル化に
ついた説明していきます。

データの時系列表示

複雑な分布を持つデータをモデル化する
場合、データを解析して何故こういった
分布になるのかを考察する必要があります。

そこで、こちらのJuliaコードを実行して、
全センサ値を時系列順にプロットした
グラフにしてみます。
src/prob_sats/complex_dist/time_series_600mm/draw_time_series_600.jl

module DrawTimeSeries600
    using Plots, DataFrames, CSV
    pyplot()

    function main(is_test=false)
        data_path = joinpath(split(@__FILE__, "src")[1], "data/sensor_data_600.txt")
        df_600_mm = CSV.read(data_path, DataFrame, 
                             header=["date", "time", "ir", "lidar"],
                             delim=' ')
        
        time_idx = Array(0:length(df_600_mm.lidar)-1)

        plot(time_idx, df_600_mm.lidar, lebel="time series")

        if is_test == false
            save_path = joinpath(split(@__FILE__, "src")[1], "src/prob_stats/complex_dist/time_series_600mm/time_series_600_mm.png")
            savefig(save_path)
        end
    end
end

データを読み込んだ際に作成される
データフレームでは縦に時系列順で
並んでいるので、それと同じ長さの
0始まりの連番配列を11行目で作成
しています。

そして、それを横軸としてセンサ値の
プロットはこちらのようになります。
f:id:sy4310:20210511235851p:plain

グラフを見ると、時間が経過する中で
センサ値が上下に変動していることが
分かります。
このことから、もう少し時間に注目して
解析してみると何かが分かりそうです。

時間帯ごとの平均値

先程の解析から時間帯が関係してそうな
ことが分かったので、今度は時間ごとに
センサ値をグループ分けして、各グループ
の平均値を計算してグラフにしてみます。

この解析をするためには、こちらの
Juliaコードを実行します。
src/prob_sats/complex_dist/group_by_hour/group_by_hour.jl

module GroupByHour
    using Plots, DataFrames, CSV, Statistics
    pyplot()

    function main(is_test=false)
        data_path = joinpath(split(@__FILE__, "src")[1], "data/sensor_data_600.txt")
        df_org = CSV.read(data_path, DataFrame, 
                      header=["date", "time", "ir", "lidar"],
                      delim=' ')
        
        hours = [Int64(floor(e/10000)) for e in df_org.time]

        # add hour array to df
        df_h = DataFrame(hour=hours)
        df_new = hcat(df_org, df_h)

        # sort
        df_sort = sort(df_new, "hour")

        # grouping
        df_grp = groupby(df_sort, "hour")
        
        # mean by group
        grps = Array(0:23)
        ovsbs = [mean(grp.lidar) for grp in df_grp]
        plot(grps, ovsbs, label="observation")
        
        if is_test == false
            save_path = joinpath(split(@__FILE__, "src")[1], "src/prob_stats/complex_dist/group_by_hour/group_by_hour.png")
            savefig(save_path)
        end
    end
end

データ読み込み時に作成するデータフレームには
センサ値の取得時刻を時分秒の6桁形式でtimeと
いう列に格納しています。

これを11行目のコードのように、10000で割って
整数に丸めた時の部分だけを残した値に変換し、
14, 15行目のコードによりhourという列として
新たに追加します。

そして作成された新しいデータフレームdf_new全体を
hour列の値に基づいて昇順にソート(18行目)、
グループ化(21行目)し、グループごとの平均値を
格納した配列を作ります(25行目)。

最後にその平均値をソートしたグループごとに
グラフにするとこちらのようになります。
f:id:sy4310:20210513234511p:plain

横軸はグループ化された時間帯ですが、
早朝の5時、6時あたりでセンサ値が最も
大きくなり、昼過ぎの14時、15時で
逆に最小になっていることが分かります。

時間帯での条件付け

先程の時間帯ごとの平均値のグラフにて
特徴的だった朝の6時台と昼過ぎの14時台
のセンサ値を抽出して、更に解析を進めて
みます。

今度はこちらのJuliaコードを使い、
グループ化したデータフレームから
6時台と14時台のセンサ値を抽出し、
ヒストグラムにしてみます。
src/prob_sats/complex_dist/hist_grp_by_hour/hist_grp_by_hour.jl

module HistGrpByHour
    using Plots, DataFrames, CSV, Statistics
    pyplot()

    function main(is_test=false)
        data_path = joinpath(split(@__FILE__, "src")[1], "data/sensor_data_600.txt")
        df_org = CSV.read(data_path, DataFrame, 
                      header=["date", "time", "ir", "lidar"],
                      delim=' ')
        
        hours = [Int64(floor(e/10000)) for e in df_org.time]

        # add hour array to df
        df_h = DataFrame(hour=hours)
        df_new = hcat(df_org, df_h)

        # sort
        df_sort = sort(df_new, "hour")

        # grouping
        df_grp = groupby(df_sort, "hour")

        # histogram
        df_grp_6 = df_grp[6]
        df_grp_14 = df_grp[14]
        histogram(df_grp_6.lidar, color=:blue, label="6")
        histogram!(df_grp_14.lidar, color=:orange, label="14")

        if is_test == false
            save_path = joinpath(split(@__FILE__, "src")[1], "src/prob_stats/complex_dist/hist_grp_by_hour/hist_grp_by_hour.png")
            savefig(save_path)
        end
    end
end

データを読み込んで時間帯ごとに
グループ化したデータフレームを
作るとこまでは、既に紹介した
コードと同じですが、24～25行目で
6時台と14時台にあたるデータを取り
出し、26～27行目でそのうちのlidarの
センサ値をヒストグラムにします。

f:id:sy4310:20210514221000p:plain
ヒストグラムにした結果、最初と
比べてそれぞれの時間帯のセンサ値が
ガウス分布に近い形状になりました。

このように、ある変数 $x$ で条件付けた
別の変数 $y$ の確率分布は、
f:id:sy4310:20210514222335p:plain
と表せます。この確率を条件付き確率と
呼ばれます。

同時確率分布

ここまでの解析結果から、時間帯ごとの
確率分布はガウス分布になると考えられる
ので、モデル化のために時間帯ごとの
ガウス分布を作ります。

その準備として、まずはセンサ値がa,
時刻がbということが同時に起こる
確率 $P(z=a, t=b)$ を求めます。
こちらのコードを実行してみましょう。
src/prob_stats/complex_dist/group_hour_heatmap/group_hour_heatmap.jl

module GrpHrHeatmap
    using DataFrames, CSV, FreqTables, Plots
    pyplot()

    function main(is_test=false)
        data_path = joinpath(split(@__FILE__, "src")[1], "data/sensor_data_600.txt")
        df_org = CSV.read(data_path, DataFrame, 
                      header=["date", "time", "ir", "lidar"],
                      delim=' ')
        
        hours = [Int64(floor(e/10000)) for e in df_org.time]

        # add hour array to df
        df_h = DataFrame(hour=hours)
        df_new = hcat(df_org, df_h)

        # pivot
        pivot = freqtable(df_new, :lidar, :hour)
        println(pivot)

        #heatmap
        x = names(pivot, 2)
        y = names(pivot, 1)
        freqs = [value for(name, value) in enumerate(pivot)]
        probs = freqs ./ length(df_org.lidar)
        heatmap(x, y, probs)
        
        if is_test == false
            save_path = joinpath(split(@__FILE__, "src")[1], "src/prob_stats/complex_dist/group_hour_heatmap/group_hour_heatmap.png")
            savefig(save_path)
        end
    end
end

まずは2行目でFreqTablesパッケージを追加で
取り込みます。これにより、データフレームに
ある2つのデータを指定して、2次元の頻度値
テーブルを作ってくれます。詳しい使い方は
こちらの参照ください。
github.com

18, 19行目のコードにより時間帯ごとのセンサ値の
頻度テーブルを作成したものはこちらのようになります。
f:id:sy4310:20210514233242p:plain

そしてこの頻度値から時間帯ごとの確率を計算し、
それをヒートマップで表現してみます。
22, 23行目でテーブルの縦軸データと横軸データを
取り出していますが、これらはNamedArrayという
文字列配列なので、names()を通すことで整数の
配列に変換する必要があります。

次に、テーブルの頻度値を24, 25行目で1次元の
配列として取り出し、その各要素を全データ数で
割ることで確率を計算します。
最後に、26行目でカラーマップに確率を設定
したヒートマップを作成します。
f:id:sy4310:20210515114553p:plain

横軸が時間、縦軸がセンサ値となっており、
色が明るい部分が確率が高いことを示します。

また、頻度値テーブルではなく、こちらの
コードを使うことで確率分布のカーネル密度
グラフを作成することができます。

marginalkde(df_new.hour, df_new.lidar, c=:ice, 
            xlabel="Hour", ylabel="LiDAR")

f:id:sy4310:20210515134843p:plain
club.informatix.co.jp

このグラフの右側と上側に表示されて
いるのは、センサ値ごと、時間帯ごとに
$P(z,t)$ を合計したグラフであり、
それぞれ $P(z)$ 、 $P(t)$ となります。
f:id:sy4310:20210515135339p:plain

周辺確率

先程の頻度値テーブルから、 $P(z)$ 、 $P(t)$ の
値を計算してみます。

まず $P(t)$ はこちらのコードで計算できます。
src/prob_stats/complex_dist/prob_sum_hour/prob_sum_hour.jl

module ProbSumHour
    using DataFrames, CSV, Plots, FreqTables
    pyplot()

    function main(is_test=false)
        data_path = joinpath(split(@__FILE__, "src")[1], "data/sensor_data_600.txt")
        df_org = CSV.read(data_path, DataFrame, 
                      header=["date", "time", "ir", "lidar"],
                      delim=' ')
        
        hours = [Int64(floor(e/10000)) for e in df_org.time]

        # add hour array to df
        df_h = DataFrame(hour=hours)
        df_new = hcat(df_org, df_h)

        # pivot
        freqs = freqtable(df_new, :lidar, :hour)
        probs = freqs / length(df_org.lidar)

        # sum
        probs_sum = [sum(probs[begin:end, i]) for i in Array(1:24)]
        plot(Array(0:23), probs_sum, label="Probability")
        println("Sum of probability = $(sum(probs))")
        
        if is_test == false
            save_path = joinpath(split(@__FILE__, "src")[1], "src/prob_stats/complex_dist/prob_sum_hour/prob_sum.png")
            savefig(save_path)
        end

        return sum(probs)
    end
end

18, 19行目で頻度値テーブルの作成、
確率テーブルへの変換を行い、22~24
行目で各時間帯での確率の合計を
計算しています。それをプロットした
グラフがこちらです。
f:id:sy4310:20210515142144p:plain

次に、上記のコードの22~24行目を
こちらのコードに書き換えれば、
$P(z)$ を計算し、グラフをプロット
できます。

keys = names(probs, 2)
probs_sum = [sum(probs[begin:end, Name(key)]) for key in keys]
plot(keys, probs_sum, label="Probability")
println("Sum of probability = $(round(sum(probs)))")

f:id:sy4310:20210515155114p:plain

これらのように、 $P(z, t)$ から
次元を落とした $P(z)$ や $P(t)$ の
値を周辺確率といい、それを計算
するための操作を周辺化、周辺化に
よってできた分布 $P(z)$ や $P(t)$ を
周辺分布と呼びます。

確率の乗法定理

ここまでに説明してきた同時確率 $P(z,t)$ と
条件付き確率 $P(z|t)$ の間には、次のような
関係が成り立ちます。
f:id:sy4310:20210515164058p:plain

これを踏まえて、先程のコードで
計算した $P(z,t)$ から $P(z|t)$ を
求めてみます。
こちらのコードを実行すると、各時間帯に
おける $P(z|t)$ が計算できます。 src/prob_stats/complex_dist/cond_z_t_bar/cond_z_t_bar.jl

# pivot
freqs = freqtable(df_new, :lidar, :hour)
probs = freqs / length(df_org.lidar)

# sum
keys = names(probs, 2)
probs_t = [sum(probs[begin:end, Name(key)]) for key in keys]

# P(z|t)
cond_z_t = probs / probs_t[1] # P(t=0)

コードにおけるprobsが $P(z,t)$ 、probs_tが $P(t)$
であり、先程の関係式から計算した $P(z|t)$ が
cond_z_tです。

そして次のコードにより、先程注目した
時間帯である6時と14時の条件付き確率
分布である $P(z|t=6)$ , $P(z|t=14)$ を
プロットしてみます。

z = names(probs, 1)

# bar plot
bar(cond_z_t[begin:end, Name(6)], xticks=z, alpha=0.3, color=:blue, label="Hour:6")
bar!(cond_z_t[begin:end, Name(14)], xticks=z, alpha=0.3, color=:red, label="Hour:14")

f:id:sy4310:20210515173646p:plain

以上で、時間帯ごとの $P(z|t)$ が導出
できました。これに加えて時間帯 $t$ が
分かると、センサ値 $z$ に対するより
詳細な予測が可能になります。

これは確率の乗法定理と呼ばれ、一般に
f:id:sy4310:20210515175634p:plain
が成り立ちます。

ベイズの定理

上記の乗法定理の式より、
f:id:sy4310:20210515222738p:plain
という式が得られます。

また、 $x$ のとりうる値の集合を
$\chi$ とすると、加法定理から
f:id:sy4310:20210515223154p:plain
が成り立ちます。

2つ目の式を1つ目の式に代入すると、
f:id:sy4310:20210515230829p:plain
が成り立ち、ベイズの定理と呼ばれます。

次のコードで、これまでのように
頻度テーブルから $P(z|t)$ を計算する
のと、ベイズの定理から $P(z|t)$ を
計算するのを比べてみましょう。
src/prob_stats/complex_dist/probs_bayes/probs_bayes.jl

今回は例として、 $P(z=630|t=13)$ の値を
求めてみます。ちなみに、データフレームを
作成するところはこれまでに紹介したコードと
同じなので割愛します。

まずは頻度テーブルから同時確率 $P(z, t)$ を求め、
そこから周辺確率 $P(t)$ , $P(z)$ をこちらの
コードで計算します。

# P(t)
freqs = freqtable(df_new, :lidar, :hour)
probs = freqs / length(df_org.lidar)
keys = names(probs, 2)
p_t = [sum(probs[begin:end, Name(key)]) for key in keys]

# P(z)
probs_trans = transpose(probs)
keys_trans = names(probs_trans, 2)
p_z = [sum(probs_trans[begin:end, Name(key)]) for key in keys_trans]

次に、ベイズの定理に従って $P(z|t)$ を計算
するのに必要な $P(t|z)$ をこちらのコードで
計算します。

# P(t|z)
for (i, key) in enumerate(keys_trans)
    probs_trans[begin:end, Name(key)] /= p_z[i]
end
cond_t_z = probs_trans

最初に求めた同時確率 $P(z,t)$ を乗法定理に
従って $P(z)$ で割ることで $P(t|z)$ を求めて
います。

ここで例として求めたいのは $P(z=630|t=13)$ なので、
ここまでに計算した各種確率分布から、 $P(z=630)$ ,
$P(t=13)$ , $P(t=13|z=630)$ を取り出します。

p_z_630 = p_z[findfirst(keys_trans .== 630)]
p_t_13 = p_t[findfirst(keys .== 13)]
p_t_13_z_630 = cond_t_z[Name(13), Name(630)]

そして最後に、これらをベイズの定理に
当てはめたこちらのコードにより、
$P(z=630|t=13)$ を求めます。

bayes = p_t_13_z_630 * p_z_630 /p_t_13

ちなみにこの値は、先に説明したこちらの
式からも求めることができます。
f:id:sy4310:20210516113537p:plain
そのコードはこちらのようになります。

cond_z_t = probs / p_t[findfirst(keys .== 13)]
answer = cond_z_t[Name(630), Name(13)]

それぞれのやり方で計算した値を
出力して比較してみましょう。
f:id:sy4310:20210516114017p:plain

下から2番目の出力がベイズの定理による
もの、一番下が同時確率から計算したもの
ですが、どちらも同じものを計算している
ので、値が一致していることが分かります。

ベイズの定理による原因推定

ベイズの定理を利用して、得られている
データから原因を推定することができます。
ここでは例として、「センサ値から計測
した時間帯を推定する」問題を扱います。

センサ値を $z$ , 時間帯を $t$ とすると、
確率 $P(t|z)$ を求めれば、特定のセンサ値が
どの時間帯に得られた確率が高いかを推定
できます。
そのためにこちらのコードを用意します。

function bayes_estimation(sensor_value, cond_z_t, current_estimation)
    new_estimation = zeros(24)
    for i in Array(1:24)
        new_estimation[i] = cond_z_t[Name(sensor_value), Name(i - 1)] * current_estimation[i]
    end
    return new_estimation / sum(new_estimation)
end

# P(t)
freqs = freqtable(df_new, :lidar, :hour)
probs = freqs / length(df_org.lidar)
keys = names(probs, 2)
p_t = [sum(probs[begin:end, Name(key)]) for key in keys]

# P(z|t)
cond_z_t = probs / p_t[findfirst(==(0), keys)]

# estimate
estimation = bayes_estimation(630, cond_z_t, p_t)
plot(estimation, label="P(t|z=630)")

f:id:sy4310:20210516123318p:plain
グラフを見ると、朝の3時頃に一番値が大きく
なっているので、センサ値 630は深夜に得られた
確率が高いと推定できます。

また、これは複数のセンサ値が続けて
得られた場合でも同様に考えることが
できます。
例えば、センサ値が630, 632, 636だと
すると、ベイズの定理では
f:id:sy4310:20210516135837p:plain
となり、各センサ値における $P(z|t)$ を
掛け合わせればいいことになります。

# estimate
observations_5 = [630, 632, 636]
estimation = p_t
for observ in observations_5
    estimation = bayes_estimation(observ, cond_z_t, estimation)
end
plot(estimation, label="P(t=5|z=630,632,636)")

このコードを実行することで、例における
3つのセンサ値が得られた場合の確率がこの
ように推定できます。
f:id:sy4310:20210516140444p:plain

例に挙げた3つのセンサ値は、実際には5時台に
得られたものなのですが、グラフを見るとその
辺りで最も確率が高くなっているので、正しい
推定ができていると考えられます。

まとめ

今回は、複雑な形状の分布をモデル化する
ための解析や前処理、確率分布モデルの
計算、ベイズの定理による原因推定の
考え方とそのためのJuliaサンプルコードを
紹介しました。

ここまでは主に1次元のデータのみを扱って
きましたが、次の記事では2次元のデータの
扱い方について紹介していきます。

EurekaMoments

ロボットや自動車の自律移動に関する知識や技術、プログラミング、ソフトウェア開発について勉強したことをメモするブログ

Juliaで学ぶ確率ロボティクス～複雑な分布のモデル化とベイズの定理～

目次