コンプガチャのシミュレーション

前回の記事で、コンプガチャの期待値と分散を求めましたが、いまいち自信がなかったのでシミュレーションしてみました。
参考: 全種類の景品を集めるのに必要な回数の期待値

おさらいしておくと、$n$種類の景品があるクジを景品が全種類揃うまで引く回数は、
期待値が$n\sum_{k=1}^{n}\frac{1}{k}$, 分散が$n\sum_{k=1}^{n-1}\frac{k}{(n-k)^2}$です。

実際にそのような結果になるのか、仮に$n=20$として、プログラムで繰り返し実行してみましょう。
ちなみに、$n=20$の場合の期待値と分散は次ようになります。


import numpy as np


n = 20
# 期待値
print(sum(n/np.arange(1, n+1)))
# 71.95479314287363

# 分散
print(sum([n*k/(n-k)**2 for k in range(1, n)]))
# 566.5105044223357

シミュレーションに使うために、景品が全種類揃うまでクジを引く関数を実装します。


def complete_gacha(n):
    # 揃ったアイテムの配列
    item_list = []

    # n種類揃うまでクジを引く
    while len(set(item_list)) < n:
        item_list.append(np.random.randint(n))

    return item_list

ためしに$n=5$で実行すると次のような結果が得られます。


print(complete_gacha(5))
# [4, 2, 0, 2, 4, 1, 3]

それでは、この関数を100000回実行し、回数(=返された配列の長さ)のリストを作って、平均値と不変分散を出してみましよう。


result_list = np.array([len(complete_gacha(20)) for _ in range(100000)])

# 期待値
print(result_list.mean())
# 71.86858

# 不偏分散
print(result_list.var(ddof=1))
# 566.9716985005849

試行回数がかなり大きいのもあって、理論値にかなり近い結果が得られましたね。
どうやら前回の記事の結果は一応正しそうです。

全種類の景品を集めるのに必要な回数の期待値

いわゆるコンプガチャの問題です。

$n$種類の景品があるクジにおいて、全ての景品を揃えるためには何回程度引けば良いのか(=期待値)を考えていたところ、
うまく解けたので記事にすることにしました。

まず改めて問題の前提条件を整理しておきます。
– クジには$n$種類の景品がある。
– どの景品も当たる確率は等しく$1/n$である。
– クジは無限にあり、過去の景品は将来のクジの確率に影響しない。
– 全集類の景品を最低1回引くまでクジを引き続け、その回数の期待値を求める。

3番目の条件は重要です。要するに景品Aがなかなか出なかったからといって、そのあと景品Aを引ける確率が上がったりしないということです。

当初、場合分けを色々考えてアプローチしていたのですが、次のように考えるとすんなり解けました

この問題は、全ての景品が揃うまでのクジの回数を次のように分解して考えます。

全ての景品が揃うまでの回数 =
1種類目の景品を引くまでの回数
+ 1種類持っている状態から2種類目の景品を引くまでの回数
+ 2種類持っている状態から3種類目の景品を引くまでの回数
・・・
+ $n-1$種類持っている状態から$n$種類目の景品を引くまでの回数

すると、期待値の線形性から次のようになります。
$$
E(\text{全ての景品が揃うまでの回数}) = \sum_{k=1}^{n}E(\text{k-1種類持っている状態からk種類目の景品を引くまでの回数})
$$

あとは、$k-1$種類持っている状態から$k$種類目の景品を引くまでの回数の期待値がわかれば良いです。
ここで、$k-1$種類持っているということは、持っていない景品は$n-k+1$種類であり、
これは、$\frac{n-k+1}{n}$の確率で当たり(=まだ持ってない景品)を引けるクジと考えることができます。
そして、そのあたりを引くまでの回数は、$p=\frac{n-k+1}{n}$の幾何分布に従うので、
前回の記事で見た通り、その期待値は$\frac{1}{p}=\frac{n}{n-k+1}$となります。

よって、
$$
E(\text{k-1種類持っている状態からk種類目の景品を引くまでの回数}) = \frac{n}{n-k+1}
$$
ですから、
$$
E(\text{全ての景品が揃うまでの回数}) = \sum_{k=1}^{n}\frac{n}{n-k+1} = \frac{n}{n} + \frac{n}{n-1} + \cdots + \frac{n}{1}
$$
となり、和の順番を入れ替えて$n$で括ると、
$$
\begin{align}
E(\text{全ての景品が揃うまでの回数}) &= n\left(1+\frac{1}{2}+\frac{1}{3}+\cdots \frac{1}{n}\right)\\
&= n\sum_{k=1}^{n}\frac{1}{k}
\end{align}
$$
となります。
シンプルで美しい結果になりましたね。

ついでにですが、分散も求めておきましょう。
$X$と$Y$が独立の時、期待値同様に分散も$V(X+Y)=V(X)+V(Y)$と分解できることに注意すると以下のようになります。
$$
V(\text{全ての景品が揃うまでの回数}) = \sum_{k=1}^{n}V(\text{k-1種類持っている状態からk種類目の景品を引くまでの回数}).
$$
パラメーターが$p$の幾何分布の分散は、$\frac{1-p}{p^2}$ですから、$p=\frac{n-k+1}{n}$を代入すると、
$$
\begin{align}
V(\text{k-1種類持っている状態からk種類目の景品を引くまでの回数}) &= \frac{1-\frac{n-k+1}{n}}{\left(\frac{n-k+1}{n}\right)^2}\\
&=\frac{n(k-1)}{(n-k+1)^2}
\end{align}
$$
となります。
よって、求めたい分散は、
$$
\begin{align}
V(\text{全ての景品が揃うまでの回数}) &= \sum_{k=1}^{n}\frac{n(k-1)}{(n-k+1)^2}\\
&= n\sum_{k=1}^{n-1}\frac{k}{(n-k)^2}
\end{align}
$$
となります。

幾何分布の期待値と分散

この次の記事で、幾何分布の性質(期待値)を使いたいのでおさらいしておきます。

おさらい:
成功確率が$p$である独立なベルヌーイ試行を繰り返す時、初めて成功するまでの試行回数$X$が従う確率分布を幾何分布と言います。
(本やサイトによっては、初めて成功する回数ではなく、初めて成功するまでの失敗回数で定義することもあります。
負の二項分布との関係が明確になったり、台が0始まりになったりするので実は個人的にはそちらの方が好みです。)

確率質量関数は次のようになります。
$$
P(X=k) = (1-p)^{k-1}p \quad (k = 1, 2, 3, \cdots).
$$
確率$1-p$で発生する失敗を$k-1$回続けた後に、確率$p$で発生する成功を1回、と考えれば自明ですね。

期待値は$\frac{1}{p}$、分散は$\frac{(1-p)}{p^2}$です。

モーメント母関数は次の式になります。
$$
M_X(t) = \frac{pe^t}{1-(1-p)e^t} \qquad(t<-\log{(1-p)}). $$ 期待値と分散はこのモーメント母関数から算出することもできるのですが、見ての通り、この関数の微分、2回微分を計算していくのは結構手間です。 なので、幾何分布に関しては、期待値と分散は直接計算するのがおすすめです。 (といってもこれもそこそこトリッキーなことをするのですが。) では、期待値から導出していきましょう。まず期待値の定義です。 $$ E(X) = \sum_{k=1}^{\infty}k(1-p)^{k-1}p. $$ ここで、$\frac{1}{1-y}$のマクローリン展開を考えます。 $$ \frac{1}{1-y} = 1+y+y^2+\cdots = \sum_{k=0}^{\infty}y^k. $$ これを両辺$y$で微分すると次の式になります。(しれっと微分と極限の順序交換をしています。数学科の学生さんなどはこの辺りも厳密に議論することをお勧めします。) $$ \frac{1}{(1-y)^2} = \sum_{k=1}^{\infty}ky^{k-1}. $$ この式に、$y=1-p$を代入すると次のようになります。 $$ \sum_{k=1}^{\infty}k(1-p)^{k-1} = \frac{1}{(1-1+p)^2} = \frac{1}{p^2}. $$ 両辺に$p$をかけることで、 $$ E(X) = \frac{1}{p} $$ が導かれました。 成功率が$p=\frac{1}{n}$のベルヌーイ試行は、平均$\frac{1}{p}=n$回で成功する、と考えると直感とよくあいますね。 続いて、分散$V(X)$を導出するために、$E(X^2)$を計算していきましょう。 先ほどのマクローリン展開の微分の式から始めます。 $$ \frac{1}{(1-y)^2} = \sum_{k=1}^{\infty}ky^{k-1}. $$ この両辺に、$y$を掛けます。 $$ \sum_{k=1}^{\infty}ky^{k} = \frac{y}{(1-y)^2} $$ そして、この両辺をもう一回$y$で微分します。 $$ \sum_{k=1}^{\infty}k^2y^{k-1} = \frac{1+y}{(1-y)^3}. $$ $y=1-p$とすると、 $$ \sum_{k=1}^{\infty}k^2(1-p)^{k-1} = \frac{2-p}{p^3}. $$ 両辺に$p$をかけて、 $$ \sum_{k=1}^{\infty}k^2(1-p)^{k-1}p = \frac{2-p}{p^2}. $$ この左辺は$E(X^2)$ですね。 よって、分散$V(X)$は、次のように求まります。 $$ \begin{align} V(X) &=E(X^2) - E(X)^2\\ &=\frac{2-p}{p^2} - \left(\frac{1}{p}\right)^2\\ &=\frac{1-p}{p^2}. \end{align} $$

勤務先のテックブログの宣伝

現在、週2回のペースでこの(私用の)ブログを更新していますが、実はここ以外にも勤務先であるオープンワーク社のテックブログにも記事を投稿しています。
投稿頻度はこのブログに比べて非常に低く、僕はまだ4記事しか投稿していないのですが、
その分、1記事1記事は丁寧に時間をかけて書いてきました。

こちらのブログには書いてこなかった、実際に仕事でやっている内容と密接した記事も書いていますので、
もしご興味のあるかたがいらっしゃいましたらこれらの記事も読んでいただけると嬉しいです。

僕が投稿した記事は以下の4記事になります(新しい順。)

OpenWorkの年齢別年収機能の裏側
ABテストの目的と分析時にアナリストが考えていること
企業の”採用力”を指標化しようとして失敗した話
オープンワークのアナリストが分析していること

次の更新は未定ですが、また新しい記事を更新したらこちらのブログでも紹介させていただこうと思います。

二項分布のモーメント母関数とそれを用いた期待値と分散の導出

前回の記事で二項分布の期待値と分散を直接計算したわけですが、記事中でも述べている通り、
二項分布の期待値や分散を導出するのはモーメント母関数を使った方が楽です。
マイナーな方法だけ紹介しているというのも変なので、この記事で二項分布のモーメント母関数について紹介します。

早速ですが、確率関数は
$$P[X=k] = {}_{n}\mathrm{C}_{k}p^k(1-p)^{n-k}$$
なので、モーメント母関数は次のようになります。
$$
\begin{align}
M_X(t) &= E(e^{tX})\\
&= \sum_{k=0}^{n} e^{tk}{}_{n}\mathrm{C}_{k}p^k(1-p)^{n-k}\\
&= \sum_{k=0}^{n} {}_{n}\mathrm{C}_{k} (pe^{t})^k(1-p)^{n-k}\\
&= (pe^t + 1 -p)^n.
\end{align}
$$
最後の行の式変形は二項定理を使いました。

モーメント母関数を$t$で1回微分すると次式になります。
$$
\frac{d}{dt} M_X(t) = npe^t(pe^t + 1 -p)^{n-1}.
$$
$t=0$を代入することで期待値が得られます。
$$
\begin{align}
E(X) &= \left.\frac{d}{dt} M_X(t)\right|_{t=0}\\
&=np(p+1-p)\\
&=np.
\end{align}
$$
前回の記事の直接計算するのに比べて若干楽なのが感じられると思います。

続いて分散です。モーメント母関数の2回微分は次のようになります。
$$
\frac{d^2}{dt^2} M_X(t) = npe^t(pe^t + 1 -p)^{n-1} + n(n-1)p^2e^{2t}(pe^t + 1 -p)^{n-2}.
$$
正確には$n=1$の場合と$n\geq2$の場合でそれぞれ計算しないといけないのですが、結局どちらの場合も上の式で表されることが証明できます。

ここから二項分布の2次のモーメントが次のように計算できます。
$$
\begin{align}
E(X^2) &= \left.\frac{d^2}{dt^2} M_X(t)\right|_{t=0}\\
&= np(p+1-p)^{n-1} + n(n-1)p^2(p+1-p)^{n-2}\\
&= np +n^2p^2-np^2
\end{align}
$$

よって、二項分布の分散は次のように導出されます。
$$
\begin{align}
V(X) &= E(X^2)-E(X)^2\\
&= np +n^2p^2-np^2 – (np)^2\\
&= np(1-p).
\end{align}
$$

分散に関しては、直接計算するに比べてモーメント母関数を使った方がはるかに楽に導出できましたね。

二項分布の期待値と分散を定義から計算してみた

おさらい:
成功確率が$p$のベルヌーイ試行を独立に$n$回行った時の成功回数を確率変数とする分布を二項分布と呼び、$B(n, p)$と書きます。
確率関数は次の式になります。
$$
P[X=k] = {}_{n}\mathrm{C}_{k}p^k(1-p)^{n-k}.
$$

期待値$E(X)$と分散$V(X)$は次の式で表されることが知られています。
$$
\begin{align}
E(X) &= np.\\
V(X) &= np(1-p).
\end{align}
$$

色々なテキストを見ると期待値や分散の導出はモーメント母関数を使われているのをよく見かけます。
最近復習と計算の練習を兼ねて、これらをモーメント母関数を使わずに定義から直接算出してみたところ、思ったより手こずったので記事に残すことにしました。

式変形の途中で二項係数の次の性質を使いますので注意してみてください。
$x\geq1$の時、$x\cdot{}_{n}\mathrm{C}_{n} = n\cdot {}_{n-1}\mathrm{C}_{x-1}$です。

証明 $x\geq1$とすると、
$$
\begin{align}
x\cdot{}_{n}\mathrm{C}_{x} &= x\frac{n!}{x!(n-x)!}\\
& = n\frac{(n-1)!}{(x-1)!((n-1) – (x-1))!}\\
& = n \cdot {}_{n-1}\mathrm{C}_{x-1}.
\end{align}
$$

それでは、本題に戻って期待値$E(X)$から算出していきます。
$$
\begin{align}
E(X) &= \sum_{x=0}^{n}x\cdot{}_{n}\mathrm{C}_{x}p^{x}(1-p)^{n-x}\\
&= \sum_{x=1}^{n}x\cdot{}_{n}\mathrm{C}_{x}p^{x}(1-p)^{n-x}\\
&= \sum_{x=1}^{n}n\cdot{}_{n-1}\mathrm{C}_{x-1}p^{x}(1-p)^{n-x} \qquad& (\text{冒頭の二項係数の性質から})\\
&= np\sum_{x=1}^{n}{}_{n-1}\mathrm{C}_{x-1}p^{x-1}(1-p)^{n-x}\\
&= np\sum_{x=0}^{n-1}{}_{n-1}\mathrm{C}_{x}p^{x}(1-p)^{n-1-x} \qquad&(\text{x-1をxに置き換え})\\
&= np\{p+(1-p)\}^{n-1}\\
&= np.
\end{align}
$$

以上で、$B(n, p)$の期待値が$np$であることが証明できました。
つぎは分散$V(X)$ですが、$V(X)=E(X^2)-E(X)^2$を利用して算出するので、$E(X^2)$を計算していきます。
$$
\begin{align}
E(X^2) &= \sum_{x=0}^{n}x^2\cdot{}_{n}\mathrm{C}_{x}p^{x}(1-p)^{n-x}\\
&= \sum_{x=1}^{n}x^2\cdot{}_{n}\mathrm{C}_{x}p^{x}(1-p)^{n-x}\\
&= \sum_{x=1}^{n}xn\cdot{}_{n-1}\mathrm{C}_{x-1}p^{x}(1-p)^{n-x}\\
\end{align}
$$
ここで、$\sum$の中の最初の$x$を、$x=(x-1)+1$と変形して、2項にわけます。
$$
\begin{align}
E(X^2) &= \sum_{x=1}^{n}(x-1)n\cdot{}_{n-1}\mathrm{C}_{x-1}p^{x}(1-p)^{n-x}+\sum_{x=1}^{n}n\cdot{}_{n-1}\mathrm{C}_{x-1}p^{x}(1-p)^{n-x}\\
&= np\sum_{x=1}^{n}(x-1)\cdot{}_{n-1}\mathrm{C}_{x-1}p^{x-1}(1-p)^{n-x}+np\sum_{x=1}^{n}{}_{n-1}\mathrm{C}_{x-1}p^{x-1}(1-p)^{n-x}\\
&= np\sum_{x=0}^{n-1}x\cdot{}_{n-1}\mathrm{C}_{x}p^{x}(1-p)^{n-1-x}+np\sum_{x=0}^{n-1}{}_{n-1}\mathrm{C}_{x}p^{x}(1-p)^{n-1-x}
\end{align}
$$
ここで、$\sum_{x=0}^{n-1}x\cdot{}_{n-1}\mathrm{C}_{x}p^{x}(1-p)^{n-1-x}$は$B(n-1, p)$の期待値なので、$(n-1)p$です。
さらに、$\sum_{x=0}^{n-1}{}_{n-1}\mathrm{C}_{x}p^{x}(1-p)^{n-1-x}$は$B(n-1, p)$の確率関数の全体の和なので$1$になります。
よって、
$$
E(X^2) = n(n-1)p^2 + np
$$
となります。

あとはこれを使って、
$$
\begin{align}
V(X) &= E(X^2) – E(X)^2\\
&= n^2p^2-np^2+np-(np)^2\\
&= np(1-p)
\end{align}
$$
が導出されました。

2021年上半期(1月~6月)によく読まれた記事

2021年もあっという間に半分が終わってしまいました。
ここで恒例(?)のよく読まれた記事ランキングを掲載したいと思います。

参考ですが、昨年1年間のよく読まれた記事ランキングはこちらです。
参考: 2020年のまとめ

では早速発表していきます。
集計期間は2021年1月から6月まで。pvでソートしています。

  1. matplotlibのグラフを高解像度で保存する
  2. ネットワークグラフの中心性
  3. pyenvで作成した環境を消す
  4. TensorflowやKerasでJupyterカーネルが落ちるようになってしまった場合の対応
  5. Pythonで連続した日付のリストを作る
  6. numpyのpercentile関数の仕様を確認する
  7. matplotlibのデフォルトのフォントを変更する
  8. INSERT文でWITH句を使う
  9. kerasのto_categoricalを使ってみる
  10. scipyで階層的クラスタリング

相変わらず、プログラミングのちょっとした小ネタのような記事が人気を集めていますね。
このブログ自体そういう記事が多いので、やむを得ないことですが。

データサイエンティストのブログらしい記事としては、
今回ネットワークグラフの中心性の記事が2位にランクインしました。
これを書いたのは昨年なのですが、当時ネットワーク解析について色々勉強して書いた記事だったので、
ニーズがあって嬉しいです。

Amazon Rekognitionで物体の検出

前回に引き続いて、Amazon Rekognitionの話です。
Rekognitionでは顔だけでなく、画像に写っている物体やシーンについてのラベル情報を得ることができます。
参考: オブジェクトおよびシーンの検出 – Amazon Rekognition

boto3で使う方法は、顔検出の時とよく似ていて、
detect_faces() の代わりに、 detect_labels() を呼び出すだけです。
参考: Rekognition — Boto3 Docs 1.17.95 documentation

サンプルとして、ドキュメントに掲載されている車の並んだ道でスケーボーやっている人の画像(ファイル名: skateboard.jpg)でやってみます。
顔検出の場合と同様にローカルの画像ファイルを読み込む方法と、S3にアップロードされた画像を使う方法があります。

ローカルのファイルを使う場合は次のようにします。


import boto3

with open("./skateboard.jpg", "rb") as f:
    img = f.read()

client = boto3.client('rekognition')
response = client.detect_labels(Image={'Bytes': img})

S3にアップロードしたデータを使う場合は次のようにします。


client = boto3.client('rekognition')
response = client.detect_labels(
    Image={
        'S3Object': {
            'Bucket': '{バケット名}',
            'Name': 'skateboard.jpg',
        }
    },
)

結果は辞書型で戻ってきます。
キー: Labels の値がメインの結果です。


response.keys()
# dict_keys(['Labels', 'LabelModelVersion', 'ResponseMetadata'])

Labels の値は、検出できたものの配列になっています。
試しに二つほど表示すると次のようになります。


import json

print(json.dumps(response["Labels"][4], indent=4))
"""
{
    "Name": "Person",
    "Confidence": 98.37577819824219,
    "Instances": [
        {
            "BoundingBox": {
                "Width": 0.1903613954782486,
                "Height": 0.27238351106643677,
                "Left": 0.43754446506500244,
                "Top": 0.3520295023918152
            },
            "Confidence": 98.37577819824219
        },
        {
            "BoundingBox": {
                "Width": 0.037608712911605835,
                "Height": 0.06765095144510269,
                "Left": 0.9162867665290833,
                "Top": 0.50001460313797
            },
            "Confidence": 86.00637817382812
        }
    ],
    "Parents": []
}
"""
print(json.dumps(response["Labels"][6], indent=4))
"""
{
    "Name": "Pedestrian",
    "Confidence": 97.18687438964844,
    "Instances": [],
    "Parents": [
        {
            "Name": "Person"
        }
    ]
}
"""

Name にラベル名が格納され、 Instances にそれが画像のどこに含まれていたが示されています。
Instances は空の配列のこともあります。要するに画像のどこかに写っているけど、場所は不明ということです。
位置が出力されるものとそうでないものにどんな規則性があるのかはいまいちわかりませんでした。
このほか、Parentsという属性があり、親概念になるラベル名が取得されます。

さて、検出されたラベル名と、 Instances の数、 Parentsの一覧を出力してみましょう。
結構色々検出されていますね。


for label in response["Labels"]:
    print(f'ラベル名: {label["Name"]}', f'インスタンス数: {len(label["Instances"])}')
    if len(label["Parents"]) > 0:
        print(f'親ラベル: {label["Parents"]}')


"""
ラベル名: Car インスタンス数: 14
親ラベル: [{'Name': 'Vehicle'}, {'Name': 'Transportation'}]
ラベル名: Automobile インスタンス数: 0
親ラベル: [{'Name': 'Vehicle'}, {'Name': 'Transportation'}]
ラベル名: Vehicle インスタンス数: 0
親ラベル: [{'Name': 'Transportation'}]
ラベル名: Transportation インスタンス数: 0
ラベル名: Person インスタンス数: 2
ラベル名: Human インスタンス数: 0
ラベル名: Pedestrian インスタンス数: 0
親ラベル: [{'Name': 'Person'}]
ラベル名: Skateboard インスタンス数: 1
親ラベル: [{'Name': 'Sport'}, {'Name': 'Person'}]
ラベル名: Sport インスタンス数: 0
親ラベル: [{'Name': 'Person'}]
ラベル名: Sports インスタンス数: 0
親ラベル: [{'Name': 'Person'}]
ラベル名: Road インスタンス数: 0
ラベル名: Wheel インスタンス数: 10
親ラベル: [{'Name': 'Machine'}]
ラベル名: Machine インスタンス数: 0
ラベル名: Path インスタンス数: 0
ラベル名: Downtown インスタンス数: 0
親ラベル: [{'Name': 'City'}, {'Name': 'Urban'}, {'Name': 'Building'}]
ラベル名: City インスタンス数: 0
親ラベル: [{'Name': 'Urban'}, {'Name': 'Building'}]
ラベル名: Urban インスタンス数: 0
ラベル名: Building インスタンス数: 0
ラベル名: Town インスタンス数: 0
親ラベル: [{'Name': 'Urban'}, {'Name': 'Building'}]
ラベル名: Tarmac インスタンス数: 0
ラベル名: Asphalt インスタンス数: 0
ラベル名: Parking Lot インスタンス数: 0
親ラベル: [{'Name': 'Car'}, {'Name': 'Vehicle'}, {'Name': 'Transportation'}]
ラベル名: Parking インスタンス数: 0
親ラベル: [{'Name': 'Car'}, {'Name': 'Vehicle'}, {'Name': 'Transportation'}]
ラベル名: Intersection インスタンス数: 0
親ラベル: [{'Name': 'Road'}]
ラベル名: Architecture インスタンス数: 0
親ラベル: [{'Name': 'Building'}]
ラベル名: Office Building インスタンス数: 0
親ラベル: [{'Name': 'Building'}]
ラベル名: Sidewalk インスタンス数: 0
親ラベル: [{'Name': 'Path'}]
ラベル名: Pavement インスタンス数: 0
親ラベル: [{'Name': 'Path'}]
ラベル名: Neighborhood インスタンス数: 0
親ラベル: [{'Name': 'Urban'}, {'Name': 'Building'}]
ラベル名: Street インスタンス数: 0
親ラベル: [{'Name': 'City'}, {'Name': 'Road'}, {'Name': 'Urban'}, {'Name': 'Building'}]
ラベル名: Coupe インスタンス数: 0
親ラベル: [{'Name': 'Sports Car'}, {'Name': 'Car'}, {'Name': 'Vehicle'}, {'Name': 'Transportation'}]
ラベル名: Sports Car インスタンス数: 0
親ラベル: [{'Name': 'Car'}, {'Name': 'Vehicle'}, {'Name': 'Transportation'}]
ラベル名: Sedan インスタンス数: 0
親ラベル: [{'Name': 'Car'}, {'Name': 'Vehicle'}, {'Name': 'Transportation'}]
"""

結果が多すぎる場合は、 MaxLabels 引数で出力される結果の数を絞り込むこともできます。

さて、ラベル名の一覧だけみてもどのくらい正確なのかわかりにくいので、画像中に図示してみましょう。
インスタンス数が0のものは、図示できないので、インスタンスが返されたものだけ、ボックスを描いていきます。


# %%pycodestyle
import matplotlib.pyplot as plt
from matplotlib import patches
from skimage import io

# matplotlibの可視化用に画像の読み込み
img_array = io.imread("./skateboard.jpg")

# 画像の高さと幅の取得
image_h, image_w, _ = img_array.shape

fig = plt.figure(facecolor="w", figsize=(12, 12))
ax = fig.add_subplot(111)
# 元の画像を表示する
ax.imshow(img_array)

for label in response["Labels"]:
    for instance in label["Instances"]:
        left = instance["BoundingBox"]["Left"] * image_w
        top = instance["BoundingBox"]["Top"] * image_h
        width = instance["BoundingBox"]["Width"] * image_w
        height = instance["BoundingBox"]["Height"] * image_h
        patch = patches.Rectangle(
            xy=(left, top),
            width=width,
            height=height,
            fill=False,
            ec="c",
            linewidth=2
        )
        ax.add_patch(patch)
        ax.text(x=left, y=top, s=label["Name"], fontsize=15, c="b")

出力された結果がこちらです。

人や車などの位置がしっかり検出できていますね。

Amazon Rekognitionで顔検出

AWSの機械学習・画像認識サービスであるAmazon Rekognition を試してみたので記録を残しておきます。
今回はまず、画像中の人の顔を検出するタスクをやってみました。
また、例によってPython(boto3)を使っています。

使い方はめちゃくちゃ簡単で、boto3のクライアントAPIから、detect_faces()というメソッドを呼び出すだけでした。
ドキュメントはこちらです。
参照: Rekognition — Boto3 Docs 1.17.93 documentation

対象のデータはローカルのファイルをバイト列のデータとして渡す方法と、S3にアップロードしてそのバケット名とファイル名を渡す方法の2種類があります。
ドキュメントには、
The input image as base64-encoded bytes or an S3 object.
と書いてあるので、Base64エンコードしないといけないのかと思ったのですが、これはどうやらドキュメントの誤りです。
Base64エンコードして渡すと逆にエラーになりますので、バイト列で読み込む時はファイルを読み込んだバイトデータをそのまま渡してください。
(予想ですが、boto3のライブラリが内部処理で Base64エンコードしてくれてると思います。)
せっかくこのブログでもBase64エンコーディングの方法を紹介する記事を書いて準備していたのにいらなかったですね。

引数の渡し方は少し特殊で、名前付き引数に辞書型で渡す必要があります。
ローカルのファイルを使う場合は次のようにします。

ちなみに画像は、AWSのコンソールでサンプルとして表示される家族写真を使います。(ファイル名:family.jpg)


import boto3


# 画像データを読み込む。
with open("./family.jpg", "rb") as f:
    img = f.read()

client = boto3.client('rekognition')
response = client.detect_faces(Image={'Bytes': img})

S3にアップロードしたデータを使う場合は次のようにします。


client = boto3.client('rekognition')
response = client.detect_faces(
    Image={
        'S3Object': {
            'Bucket': '{バケット名}',
            'Name': 'family.jpg',
        }
    },
)

結果は、辞書型で帰ってきます。FaceDetailsというキーの中身が、メインの検出結果で、
ResponseMetadataのほうはRequestのIdや、HTTPレスポンスのステータスコード、処理を実行した時刻などのメタデータが入ってます。


print(response.keys())
# dict_keys(['FaceDetails', 'ResponseMetadata'])

response[“FaceDetails”] の中身は配列で、検出された顔一人分ごとに辞書型で検出された情報が入っています。
今回のサンプルでは3人検出されているのですが、全部表示すると長いので一人分お見せすると次のようになります。
(整形のためにjsonライブラリ使います)


import json
print(json.dumps(response["FaceDetails"][0], indent=4))
"""
{
    "BoundingBox": {
        "Width": 0.1937681883573532,
        "Height": 0.3873019516468048,
        "Left": 0.2916979193687439,
        "Top": 0.13570082187652588
    },
    "Landmarks": [
        {
            "Type": "eyeLeft",
            "X": 0.34084638953208923,
            "Y": 0.2765427529811859
        },
        {
            "Type": "eyeRight",
            "X": 0.4209189713001251,
            "Y": 0.31195494532585144
        },
        {
            "Type": "mouthLeft",
            "X": 0.32429391145706177,
            "Y": 0.40175312757492065
        },
        {
            "Type": "mouthRight",
            "X": 0.39117804169654846,
            "Y": 0.43135175108909607
        },
        {
            "Type": "nose",
            "X": 0.3650367856025696,
            "Y": 0.3684481084346771
        }
    ],
    "Pose": {
        "Roll": 17.15113067626953,
        "Yaw": -3.947751760482788,
        "Pitch": -1.8470479249954224
    },
    "Quality": {
        "Brightness": 62.19182586669922,
        "Sharpness": 78.64350128173828
    },
    "Confidence": 99.99921417236328
}
"""

BoundingBox の中にあるのが、顔が検出された位置です。顔を囲む長方形の情報が含まれています。
Landmarksの下に、目、鼻、口の両端の位置が含まれます。
Poseは顔の向きです。
Qualityは画像の明るさなどの情報で、Confidenceは境界ボックス内に顔が含まれている信頼度になります。
詳しくはこちら
参考: イメージ内の顔を検出する – Amazon Rekognition

顔やそのパーツの位置の座標が0〜1の範囲に収まっていることから分かる通り、これらは
画像の左上を(0%,0%)、右下を(100%, 100%)とした時の相対的な位置を示しています。

上のJSON型データを見てもどのくらい正確に検出できているかわからないと思うので可視化してみましょう。
使い慣れているので、Matplotlibでやってみました。
別途、skimageで画像をNumpy配列として読み取り、画像の幅と高さを取得しています。
そして、それをRekognitionで取得した相対的な位置と掛け合わせることで、絶対値での座標に変換しています。


import matplotlib.pyplot as plt
from matplotlib import patches
from skimage import io


# matplotlib表示用に画像を配列で読み込み
img_array = io.imread("./family.jpg")
# 画像の高さと横幅を取得
image_h, image_w, _ = img_array.shape

fig = plt.figure(facecolor="w", figsize=(12, 12))
ax = fig.add_subplot(111)
# 元の画像を表示する
ax.imshow(img_array)

for face_detail in response["FaceDetails"]:
    # 検出された顔の位置を取得し、座標に変換する
    left = face_detail["BoundingBox"]["Left"] * image_w
    top = face_detail["BoundingBox"]["Top"] * image_h
    width = face_detail["BoundingBox"]["Width"] * image_w
    height = face_detail["BoundingBox"]["Height"] * image_h
    # 取得した座標の位置に長方形を描写する
    patch = patches.Rectangle(
        xy=(left, top),
        width=width,
        height=height,
        fill=False,
        ec="w",
        linewidth=2,
    )
    ax.add_patch(patch)

    # 目、鼻、口の両端に点をプロット
    ax.scatter(
        [landmark["X"] * image_w for landmark in face_detail["Landmarks"]],
        [landmark["Y"] * image_h for landmark in face_detail["Landmarks"]],
        c="w",
        s=3
    )

このコードで出力されるのが次の画像です。

3人分の顔が精度良く検出できていることがわかりますね。

さて、detect_faces()ですが、実はもう一つ引数を持っています。
それが、Attributes です。
Attributes=[“DEFAULT”] (こちらがデフォルト)
または、
Attributes=[“ALL”]
と指定します。 []も必須です。


response = client.detect_faces(
    Image={'Bytes': img},
    Attributes=["ALL"]
)

のように、 [“ALL”]を指定すると、取得できる情報が一気に増えます。
表情(笑顔かどうか)や、メガネやサングラスの有無、髭の有無や口が開いているかどうか、
大まかな年齢の推定なども行ってくれます。
また、目や口の位置情報はより詳細になり、輪郭や眉毛などに関する位置も取得されます。

結果がものすごく大きくなるのでこの記事には載せませんが、ぜひ一度試してみてください。

Matplotlibで多角形や円などの図形を描写する

Matplotlibのグラフに図形を入れる方法の紹介です。
(強調したい部分に丸や四角で目印をつけたり矢印を引いたりできます。)

Matplotlibのグラフに図形を入れるには、 matplotlib.patches の下に定義されている各種クラスを使います。
長方形は、 matplotlib.patches.Rectangle,
円は、 matplotlib.patches.Circle, など、図形に応じたクラスが用意されているので、
それぞれインスタンスを作成し、
ax.add_patch() で作ったインスタンスをグラフに挿入します。

Rectangle は左下の座標と幅と高さ、Circleは中心の座標と半径など、それぞれ固有のオプションがあり、それを指定することで思い通りの位置とサイズの図形を作れます。
詳しくはドキュメントの各クラスの説明をご参照ください。
参考: matplotlib.patches — Matplotlib 3.4.2 documentation

これらのクラスは全て、
matplotlib.patches.Patch
を継承して実装されています。
塗りつぶしの色や線のスタイルの指定、模様を入れるなど、汎用的な引数の説明は、matplotlib.patches.Patchのページに説明があるのでこちらも合わせて参照すると良いでしょう。

全てを紹介はしませんが、いくつかの図形を実際に書いてみたコードが以下です。


import matplotlib.pyplot as plt
from matplotlib import patches


fig = plt.figure(facecolor="w")
ax = fig.add_subplot(1, 1, 1, aspect="equal")

# 長方形
patch = patches.Rectangle(
    xy=(20, 10),  # 左下の頂点の座標
    width=30,  #  長方形の幅
    height=50,  # 長方形の高さ
    angle=10,  # 傾き
    facecolor="b",  # 塗りつぶしの色
    edgecolor="c",  # 辺の色
    linewidth=3,  # 辺の線幅
    hatch=".",  # 塗りつぶしの模様 {'/', '\', '|', '-', '+', 'x', 'o', 'O', '.', '*'}
    # 辺のスタイル  {'-', '--', '-.', ':', '', (offset, on-off-seq), ...}
    linestyle="-.",
    fill=True,  # 塗りつぶしあり。Falseにすると塗りつぶし無し。
)
ax.add_patch(patch)

# 円
patch = patches.Circle(
    xy=(80, 40),  #  中心
    radius=20,  # 半径
    fill=False,  # 塗りつぶし無し
)
ax.add_patch(patch)

# 矢印
patch = patches.Arrow(
    x=110,  # 始点のx座標
    y=20,  # 始点のy座標
    dx=60,  # x軸方向の長さ。 x+dx が終点のx座標。
    dy=20,  # y軸方向の長さ。 y+dy が終点のy座標。
    width=30,  # 矢印の幅
)
ax.add_patch(patch)

# 楕円
patch = patches.Ellipse(
    xy=(40, 90),  # 中心
    width=50,  # 横幅
    height=20,  # 高さ
    angle=-30,  # 傾き
)
ax.add_patch(patch)

# 多角形
patch = patches.Polygon(
    # 頂点の座標を n*2 次元配列で指定
    xy=[
        [80, 100],
        [80, 80],
        [110, 70],
        [120, 90],
        [100, 110]
    ]
)
ax.add_patch(patch)

# 正多角形
patch = patches.RegularPolygon(
    xy=(140, 70),  # 中心の座標
    numVertices=7,  # 頂点の数
    radius=20,  # 半径
    orientation=10,  # 角度
)
ax.add_patch(patch)

ax.autoscale()
plt.show()

patchたちを追加した最後に、 ax.autoscale() して、挿入した図形たちがグラフの描写範囲に収まるように調整しています。
これをしないと、グラフの描写範囲がデフォルトの x座標の区間[0, 1]、y座標の区間[0, 1] のままになってしまい、せっかく描いた図形が見えなくなってしまいます。

上記のコードで、以下の図が出力されます。

参考として長方形だけ、facecolor、edgecolorなど、色々指定しましたが、
これは継承元のPatchで定義されているので、もちろん他の図形でも指定できます。