skymatix Developers Blog

株式会社スカイマティクスの開発チームによるDevelopers Blogです。

ポアソン分布とありえなさの原理

Lead Researcherの藤田です。

クラウド型ドローン測量・現地管理DXツール「くみき」には日々ユーザー様が撮影したドローンなどの空撮画像がアップロードされます。我々はそれら空撮画像を品質良く3次元復元するために日々取り組んでいるのですが、今回は「くみき」にアップロードされる画像が1セット何枚くらいになっているか見てみましょう。

下のヒストグラムはある期間に「くみき」にアップロードされた空撮画像の1セットごとの画像枚数を50枚区切りにまとめたヒストグラムです。

アップロード枚数のヒストグラム

サンプル数は926、平均は222.0、中央値は:139.5となっています。下記がPythonコード。

import numpy as np
import matplotlib.pyplot as plt
# ファイルの読み込み
path = 'data.csv'
data = np.loadtxt(path)
print(f"n:{len(data)}, mean:{data.mean()}, median:{np.median(data)}")
# レンジの設定
range=(0, 2000)
n_bins=40
# グラフ描画
result = plt.hist(data, range=range, bins=int(n_bins))
hist, bins = result[0], result[1]

撮影画像1セットごとの「くみき」へのアップロードと言うイベントをポアソン分布と仮定して分布を求めてみましょう。 ポアソン分布とは、あるイベントが一定期間の間に発生する確率を表した分布で、平均と分散が同じ値であり、不良品や交通事故の発生確率など、稀に発生するイベントに対してよく当てはまる分布です (詳しくはWikipediaを参照)。 先にPythonコードを記載します(ここでy軸は件数から確率になっている)。

from scipy.stats import poisson
# ヒストグラムのbinごとの発生件数にするため、binを新しく作り直す
new_bins = np.arange(n_bins)
# 平均もbinごとの発生件数の平均として求め直す
new_mean = np.sum(hist * (new_bins + 0.5) / len(data))
# ポアソン分布の作成
y = poisson.pmf(new_bins, new_mean)
# グラフ描画
plt.bar(new_bins, hist / len(data), label='actual probability')
plt.plot(y, c='orange', label='poission probability distortion')
# 各binの両端の平均をbinの代表値とすしてx軸ラベルに割り当てる
xlabels = list(map(str, ((bins[:-1] + bins[1:]) * 0.5).astype(np.int)))
plt.xticks(new_bins, xlabels, rotation=90)
plt.xticks(new_bins, xlabels, rotation=90)
plt.legend()

早速描画されたヒストグラムを見てみると。。。

撮影データヒストグラムとポアソン分布

あまり合っていない。。。だいぶ少ない枚数に偏っていますね。また、ポアソン分布ではほとんど発生しないはずの600枚以上の撮影セットも予想に反して多数アップロードされていることがわかります。世の中起こり得ないはずの出来事は意外と頻繁に起こるものですね。記事作成中に検索して何度か出てきた書籍、「偶然」の統計学を読んでみようと思います。

終わりに

株式会社スカイマティクスでは「くみき」の三次元復元エンジンを世界一にしたく、研究開発に取り組んでいます。今後も三次元復元に関する技術的記事を書いて行こうと思います。

www.slideshare.net