sedコマンドの使い方メモ

テキストの編集は大抵vimでやっちゃうので滅多に使わないのですが、稀にコマンドラインでファイルの置換を完結させたり、一つのテンプレートファイルから中身を一部置換したファイルを複数生成する必要が発生し、それをコマンドラインで済ませたいことがあります。
そんなときに、sed (Stream Editor) コマンドを使うのですが、その度に使い方を調べているのでここメモしておきます。

基本的な使い方は、
sed -e {編集コマンド} {入力ファイル}
です。パイプラインで利用する場合は{入力ファイル}を省略できます。
ファイル中の aaa を xxx に置換する場合は次のように書きます。

$ cat input.txt
aaa,bbb,ccc
ddd,eee,fff

$ sed -e s/aaa/xxx/g input.txt
xxx,bbb,ccc
ddd,eee,fff

-e は省略可能で、省略した場合は最初の引数が編集コマンドとみなされます。
なので、大抵の場合は -e を省略して大丈夫です。

$ cat input.txt | sed s/aaa/xxx/g
xxx,bbb,ccc
ddd,eee,fff

-e オプションは複数指定することもできて、同時に複数の置換をかけることもできます。

$ sed -e s/aaa/xxx/g -e s/bbb/yyy/g input.txt
xxx,yyy,ccc
ddd,eee,fff

編集コマンドは、シングルクオーテーション、もしくはダブルクオーテーションで囲むこともできます。(置換前後の文字列のどちらかにスペースを含む場合は、確実に囲むようしましょう。そうしないとエラーになります。)

$ sed "s/aaa/X Y G/g" input.txt
X Y G,bbb,ccc
ddd,eee,fff

置換結果を別のファイルに出力する時は、他のシェルコマンド同様にリダイレクションしてあげれば大丈夫です。僕はもっぱらその形で使います。
例えば、 input.txt の aaa を xxx に置換した output.txt を生成するには次のようにします。

$ sed "s/aaa/xxx/g" input.txt > output.txt

元のファイルをそのまま書き換えることもでき、その場合は
-i {拡張子} オプションをつけます。
すると、{元のファイル名}{拡張子} というファイル名でバックアップを取った上で、入力ファイルを書き換えてくれます。
バックアップはいらないよという場合は、拡張子として長さ0の文字列を渡します。

$ sed -i '' -e 's/aaa/xxx/g' input.txt

この時、” をつけ忘れると、 input.txt-e という変な名前のファイルが残ってしまうので注意してください。

多くの編集コマンドを同時に実行したい場合などは、編集コマンドをまとめたファイルを用意しておき、それを -f オプションで渡すこともできます。

$ cat edit.sed
s/aaa/xxx/g
s/bbb/y y y/g
s/ccc/zzz/g

$ sed -f edit.sed input.txt
xxx,y y y,zzz
ddd,eee,fff

あとは、入力ファイルはスペース区切りで複数同時に渡すことも可能です。
-i オプションをつけている場合は、渡した入力ファイルたちがそれぞれ編集されます。
-i オプションがない場合は、それぞれのファイルを編集した結果が連結されて標準出力に返されます。

sed コマンドで文字列を置換することに関して主に知っておくべきことはこれくらいかなと思います。

Prestoで各行ごとに複数列の値の中から最大値/最小値を取得する

※Prestoと書いていますが、正確にはTreasure Dataで動かすことを念頭に置いた記事です。
ただ、この記事で紹介するGREATEST / LEAST という関数はMySQLにも実装されているようなので、MySQLでも同じように動作すると思います。

そんなに頻繁にあることでは無いのですが、DBのとあるテーブルのデータについて、行ごとに、複数列の最大値を取得したいことがありました。

これ、行と列が逆なら簡単です。MAX関数使うだけです。

SELECT
    MAX(col1),
    MAX(col2),
    MAX(col3)
FROM
    table_name

上のクエリで、列ごとに最大値が取得できます。
今回やりたいのはその逆で、行ごとの最大値が欲しいのです。
列が3つくらいであれば、CASE文で対応することもできなくは無いかなと思います。
こんなふうに。

SELECT
    id,
    CASE
        WHEN col1>=col2 AND col1>=col3 THEN col1
        WHEN col2>=col3 AND col2>=col1 THEN col2
        WHEN col3>=col1 AND col3>=col2 THEN col3
    END AS max_value
FROM
    table_name

ただ、列数が増えるとこの方法で対応するのはなかなか厄介です。(あまりやりたく無い)
それに、上のクエリでは対象の列にNULLが含まれていた場合に正常に動作しないので、NULLも考慮する必要がある場合はもっと複雑なクエリを書く必要があります。

もう少しスマートな方法としては、以前紹介した横縦変換の方法で値を縦持ちに変換して、
GROUP BY と MAX を使うこともできるかと思います。
参考: PrestoのUNNESTを利用した横縦変換

WITH
    unpivot_table AS (
        SELECT
            id,
              t.key,
              t.value
        FROM
            table_name
        CROSS JOIN UNNEST (
            array['col1', 'col2', 'col3'],
            array[col1, col2, col3]
        ) AS t (key, value)
    )
SELECT
    id,
    MAX(value) AS max_value
FROM
    unpivot_table
GROUP BY
    id

ただ、これはこれで仰々しくてちょっと嫌だなと思っていました。

それでドキュメントを調べてみると、どうやらGREATEST というメソッドが用意されていたようです。
参考: GREATEST and LEAST

これを使うと非常に話は単純で、次のクエリで行ごとに3列(co1, col2, col3)の最大値が取得できます。

SELECT
    id,
    GREATEST(
        col1,
        col2,
        col3
    ) AS max_value
FROM
    table_name

同様に、最小値を求めるLEAST も容易されています。(使用例略)

ちなみに、どの列が最大だったのかを取得できる GREATEST_BY みたいなのもあるといいなと思って探してみたのですが、流石にそれはなさそうでした。最大値と合わせてどの列が最大だったのかも欲しい場合は、上の UNNEST を使うクエリで縦持ちに変換して、MAX_BYするのが現実的かなと思います。(もしくはSQLで実行するのを諦めてPythonなどで書くか)

pandasのメソッドで、上位n件や下位n件のデータを取得する

先日紹介したbar chart raceのライブラリのドキュメントやソースコードを読んでいて、その中で nlargest というメソッドを見つけたのでその紹介です。その対となる nsmallest というメソッドもあります。

これが何をするメソッドとかというと、DataFrameやSeriesのデータの値が大きい方からn件(nlargest)や小さい方からn件(nsmallest)を取得してくれるものです。
え、sort_values() して、 head(n)やtail(n)すればいいじゃん、という声も聞こえてきそうですし、実際僕もそう思ってるのですが、多少の利点がちゃんとあるので読んでいただければ幸いです。

公式ドキュメントはこちらになります。
pandas.DataFrame.nlargest
pandas.DataFrame.nsmallest
pandas.Series.nlargest
pandas.Series.nsmallest

使い方は簡単で、Seriesの方であれば、取得したい件数を最初の引数nに渡してあげるだけ、DataFrameの方は、取得したい件数と合わせて、どの列の上位/下位を取得したのかを2つ目の引数columnsに渡してあげればOKです。

とりあえず、適当に作ったDataFrameに対して適当に列を指定して5項目ほど取得してみましょう。

import pandas as pd
import numpy as np


# 50行3列の乱数データを生成する
data = np.random.randint(1, 50, size=(50, 3))
df = pd.DataFrame(data, columns=["col1", "col2", "col3"])
print(df.shape)
# (50, 3)

print(df.nlargest(5, "col2"))
"""
    col1  col2  col3
46    30    48    28
17    47    47    31
33     9    45    30
16     3    44    33
26    16    44     2
"""

見ての通り、指定した”col2″でソートした上でその値が大きいものから順番に、5項目選択されています。

nlargest/ nsmallest にはもう一つ、keepという引数があります。これは、値が等し鋳物が複数あって、n位にランクインするものが一意に決められないときにその取り扱いを指定するものです。
“first”(デフォルト)を指定すると、元のデータで先に登場指定したものが優先され、”last”を指定すると、最後に登場したものが優先されます。また、”all”にすると、同率だったものが全部含まれます。

print(df.nlargest(5, "col2", keep="first"))
"""
    col1  col2  col3
46    30    48    28
17    47    47    31
33     9    45    30
16     3    44    33
26    16    44     2
"""

print(df.nlargest(5, "col2", keep="last"))
"""
    col1  col2  col3
46    30    48    28
17    47    47    31
33     9    45    30
40    48    44     1
26    16    44     2
"""

print(df.nlargest(5, "col2", keep="all"))
"""
    col1  col2  col3
46    30    48    28
17    47    47    31
33     9    45    30
16     3    44    33
26    16    44     2
40    48    44     1
"""

“col2″に値が44のレコードが3つ存在するのですが、”first”と”last”で選択されたレコードが違うのがわかりますね。そして”all”を指定すると3レコードとも返され、結果が6行になっています。

このkeep引数が存在することのほか、sort_values/head に比べると、速度面でも優れているそうです。

This method is equivalent to df.sort_values(columns, ascending=False).head(n), but more performant.

とドキュメントにもあります。
ソースを読んで無いので予想ですが、sort_values/headの方は最終的な結果に必要ない行まで全部ソートを完了させるに対して、nlargest/nsmallestの方は必要なデータだけ並べ替えてソートを打ち切ってるのではないかと思っています。

コードの実行例は載せませんでしたが、nsmallestもnlargestと同じように使うことができ、こちらは結果が小さい順に取得されます。

pandasのデータの順位を取得する

稀にではあるのですが、Pandasのデータ(DataFrame/Series)のデータの順位を取得したくなることがあります。
これまでは、DataFrameの列内の順位であれば、sort_valuesで並べ替えて、インデックスを振り直して、といった手順で対応することが多かったです。しかし、この方法では、値が等しい項目の扱いが少々厄介になります。また、最近、列内の順位ではなく、各行ごとに行内での順位を取得したいことがあり、ちょっと面倒だなと感じることがありました。

そこで、改めて調べてみたのですが、DataFrameもSeriesもそれぞれ、rankというメソッドを持っていて、これを使えば簡単に順位が取得できることがわかりました。
参考:
pandas.DataFrame.rank
pandas.Series.rank

使い方非常に簡単で、rank()を呼び出すだけです。適当なDataFrameでやってみます。

import pandas as pd


# 適当にデータを生成する
df = pd.DataFrame(
    {
        "col1": [20, 30, None, 20, 10, 20],
        "col2": [10, 50, 20, 20, 30, 60],
        "col3": [30, None, 60, None, 20, 80]
    }
)
print(df)
"""
   col1  col2  col3
0  20.0    10  30.0
1  30.0    50   NaN
2   NaN    20  60.0
3  20.0    20   NaN
4  10.0    30  20.0
5  20.0    60  80.0
"""

# 列内の順位を取得する
print(df.rank())
"""
   col1  col2  col3
0   3.0   1.0   2.0
1   5.0   5.0   NaN
2   NaN   2.5   3.0
3   3.0   2.5   NaN
4   1.0   4.0   1.0
5   3.0   6.0   4.0
"""

結果を見てわかる通り、順序は昇順で、値が小さいほど高順位(数値が小さい)ですね。

さて、このrank()メソッドはとても気が利いていて、多くの引数で細かく結果を制御できます。
まず、列ごとではなく、行ごとの順位が欲しい場合は、axis引数に1を渡します。
ちなみに、Seriesの方のドキュメントにも、axis引数があって、1を渡せるような記載があるのですがこれはおそらくドキュメントの誤りです。(普通にエラーになります。)

# 行内の順位を取得する
print(df.rank(axis=1))
"""
   col1  col2  col3
0   2.0   1.0   3.0
1   1.0   2.0   NaN
2   NaN   1.0   2.0
3   1.5   1.5   NaN
4   1.0   3.0   2.0
5   1.0   2.0   3.0
"""

昇順ではなく降順の順位が欲しい、という場合は、ascending にFalse を渡します。(デフォルトはTrueです。)

# 降順の順位を取得する
print(df.rank(ascending=False))
"""
   col1  col2  col3
0   3.0   6.0   3.0
1   1.0   2.0   NaN
2   NaN   4.5   2.0
3   3.0   4.5   NaN
4   5.0   3.0   4.0
5   3.0   1.0   1.0
"""

na_option という引数で、NaN値に対応する順位を指定できます。
“keep”(デフォルト) であれば、NaNのままです。
“top”にすると、最も高い順位(要するに1)がNaN値に振り分けられます。
“bottom”にすると、逆にもっとも低い順位が割り振られます。
それぞれ実行した結果が以下です。

df = pd.DataFrame(
    {"data":  [20, 30, None, 20, 10, 20]}
)
df["na_keep"] = df.data.rank(na_option="keep")
df["na_top"] = df.data.rank(na_option="top")
df["na_bottom"] = df.data.rank(na_option="bottom")

print(df)
"""
   data  na_keep  na_top  na_bottom
0  20.0      3.0     4.0        3.0
1  30.0      5.0     6.0        5.0
2   NaN      NaN     1.0        6.0
3  20.0      3.0     4.0        3.0
4  10.0      1.0     2.0        1.0
5  20.0      3.0     4.0        3.0
"""

さて、最初の方のコードの実行例で、2.5など小数点の順位のものがあるのがわかると思います。これは同率順位の項目に対して、デフォルトではその平均順位を返す設定になっているからです。
この設定は、 method 引数で制御できます。値はデフォルトの’average’の他、最小値(もっとも高順位)を採用する’min’、その逆に最大値を採用する’max’、元の配列に表示されていた順に順位がつく’first’、’min’と同じように、最小値が採用されるが、その次の順位の項目の順位が数が飛ばないように採番される’dense’の5種類の値が指定できます。
ちょっとわかりにくいと思うので実例でやってみます。

df = pd.DataFrame(
    {"data":  [20, 30, 40, 20, 10, 20, 40]}
)
df["m_average"] = df.data.rank(method="average")
df["m_min"] = df.data.rank(method="min")
df["m_max"] = df.data.rank(method="max")
df["m_first"] = df.data.rank(method="first")
df["m_dense"] = df.data.rank(method="dense")

print(df)
"""
   data  m_average  m_min  m_max  m_first  m_dense
0    20        3.0    2.0    4.0      2.0      2.0
1    30        5.0    5.0    5.0      5.0      3.0
2    40        6.5    6.0    7.0      6.0      4.0
3    20        3.0    2.0    4.0      3.0      2.0
4    10        1.0    1.0    1.0      1.0      1.0
5    20        3.0    2.0    4.0      4.0      2.0
6    40        6.5    6.0    7.0      7.0      4.0
"""

値が20の項目が3つあって順位的には、2位,3位,4位に相当するのですが、
averageであれば3、minであれば2、maxであれば4が割り振られているのが確認できましたね。firstであれば元の配列に出てきた通り、2,3,4位が当てられています。
そして、denseの結果を見ると、minと同様に20は2位になっているのですが、その次の30が、minの時は5位だったのに、denseでは欠番がなくこれが3位になっています。

あとは、あまり使わないと思うのですが、 pct という引数をTrueにすると、順位の数値ではなくパーセンタイルで結果が受け取れます。

df = pd.DataFrame(
    {"data":  [20, 30, 10, 20, 40]}
)
df["pct_false"] = df.data.rank(pct=False)
df["pct_true"] = df.data.rank(pct=True)
print(df)
"""
   data  pct_false  pct_true
0    20        2.5       0.5
1    30        4.0       0.8
2    10        1.0       0.2
3    20        2.5       0.5
4    40        5.0       1.0
"""

順位が一番低い項目が1になるのは想像通りですが、最高順位の項目は0では無いんですね。

WordPress 5.X系のブロックエディタで Prism.js を使う方法

WordPress 5系で新しくなったエディタ(ブロックエディタ、もしくはGutenbergというらしいですね)を使いにくいと感じていたので、つい最近まで4系のまま使い続けていたのですが、サポート期間終了の警告が出るようになってしまったので、諦めて5系にバージョンアップしました。

この新しいエディタには慣れるしか無いので諦めて使っていこうと思います。

実は先日のBar Chart Raceの記事はブロックエディタで書いたのですが、ソースコードのシンタックスハイライトをやってくれているPrism.js を動作させる方法がなかなかわからず苦戦したので、使い方を記録しておこうと思います。

Prismjs のページで公式な対処法を探したのですが、そこでは記載を見つけられなかったのであくまで僕はこうやって解決したという非公式な方法になります。

具体的には、次の手順でprismjsが動作してくれます。

  1. ブロックを追加するときに、「コード」のブロックを選択して追加する。
    もしくはブロック追加後に + ボタンを押して「コード」に変換する。
  2. 右ペインのメニューの「ブロック」の「高度な設定」タブを開き、追加 CSS クラスに「language-python」など、有効化したい言語のクラスを設定する。

例えば、ブロックを「コード」にしただけで、追加CSSを設定しないと次のような表示になります。

print("Hello World!")

追加CSSにlanguage-pythonを入れるとこうなります。

print("Hello World!")

公式ドキュメント等で確認できてないので少々不安ではありますが、ちゃんと動作してるように見えますね。

スクラッチでBar Chart Raceを実装(コード供養)

前回の記事で、Bar Chart Raceを作るライブラリを紹介しましたが、実は僕はこのライブラリが登場するよりも前、スクラッチでBar Chart Raceを実装したことがあります。
便利なライブラリが登場したので、今後スクラッチで作ることはおそらく無いのですが、せっかく作ったコードが勿体無いので供養も兼ねて紹介させていただこうと思います。

棒の伸びもライブラリのように滑らかな動きでは無いですし、順位の入れ替わりなどもバーが上下に滑らかに移動して入れ替わるのではなく、パッと切り替わるなど、全体的にパラパラ漫画感が強く出てる出来栄えなのであまり期待せずによろしくお願いします。

データだけではライブラリ付属のコロナウィルス感染者のデータを拝借します。僕のコードはNaNに対応できないので、NaNは0埋めしておきます。

# データだけはライブラリから拝借
import bar_chart_race as bcr
# サンプルデータ読み込み
df = bcr.load_dataset('covid19')
# NaNに対応できてないので0埋めしておく
df.fillna(0, inplace=True)

では、早速作っていきます。実装としては、matplotlibのアニメーション機能を使います。
FuncAnimation を使うので、実装としては次の記事と似ています。
参考: matplotlibの3次元プロットを回転するアニメーションで保存する

まず、パラパラ漫画の各コマを生成する関数を実装します。

import matplotlib.pyplot as plt
import matplotlib.ticker as ticker


def draw_barchart(date):
    target_row = df.loc[date]
    target_data = target_row.T.sort_values(ascending=True).tail(10)

    ax.clear()
    # 棒グラフを描写
    ax.barh(target_data.index, target_data.values)
    dx = target_data.max() / 200

    for i, (name, value) in enumerate(target_data.items()):
        # 棒の先端部に項目名を出力
        ax.text(value-dx, i, name, size=14, ha='right',
                va='bottom', color="white", weight=600)
        # 棒の先に値を出力
        ax.text(value+dx, i, f'{value:,.0f}', size=14, ha='left', va='center')

    # 日付を出力
    ax.text(1, 0.4, date.strftime("%Y-%m-%d"), transform=ax.transAxes,
            color='#777777', size=23, ha='right', weight=800)
    # x軸のメモリの設定
    ax.xaxis.set_major_formatter(ticker.StrMethodFormatter('{x:,.0f}'))
    ax.xaxis.set_ticks_position('top')
    ax.tick_params(axis='x', colors='#777777', labelsize=12)
    # y軸のメモリ(項目名)を消す
    ax.set_yticks([])

    ax.margins(0, 0.01)
    ax.grid(which='major', axis='x', linestyle='-')
    ax.set_axisbelow(True)

    # 外枠を消す
    plt.box(False)

コメントを多めに付けましたが、関数の中で順に、棒グラフを書いたり文字を挿入したりメモリを調整したりとコマを組み立てています。
上記の関数でパラパラ漫画のコマが描写できるので、それを使って、アニメーションにします。

import matplotlib.animation as animation

fig = plt.figure(figsize=(10, 6), facecolor="w")
ax = fig.add_subplot(111)

animator = animation.FuncAnimation(
    fig, draw_barchart, frames=df.index, interval=400)
animator.save('bar-chart-race.mp4', writer="ffmpeg")

これで出力されるのが次の動画です。

やっぱり全体的にカクカクなりますね。
データとデータの間を補完してコマ数をもっと増やすなどしないとなめらなかなアニメーションにならないようです。

PythonのライブラリでBar Chart Raceを作ってみた

皆さんもどこかでご覧になったことがあると思うのですが、項目ごとに増加し続けるデータの面白い可視化方法として、Bar Chart Race というものがあります。
棒グラフがグイグイ伸びて順位を争っているようなアニメーションですね。

その、Bar Chart Race をPythonで手軽に作れるライブラリを見つけたので今回の記事ではそれを紹介します。
その名も、 bar_chart_race です。 そのままですね。
公式ドキュメントはこちらになります。
Bar Chart Race

インストールはpipでもcondaでも可能です。(僕はcondaで入れました)
次の2行のコードのどちらかを実行してください。

pip install bar_chart_race
conda install -c conda-forge bar_chart_race

使い方は非常に簡単で、ライブラリをインポートしたら、bar_chart_raceというメソッドにデータと保存するファイル名を渡すだけです。

実際にデータを渡す前に、どんなデータを渡せばいいのか確認しておきましょう。
このライブラリ自体にサンプルデータとして国別のコロナウィルス感染者数のデータが同梱されているので、それをみてみます。

import bar_chart_race as bcr


# サンプルデータ読み込み
df = bcr.load_dataset('covid19')

# データのサイズ
print(df.shape)
# (57, 20)

print(df.index[: 5])
"""
DatetimeIndex(['2020-02-26', '2020-02-27', '2020-02-28', '2020-02-29',
               '2020-03-01'],
              dtype='datetime64[ns]', name='date', freq=None)
"""
print(df.columns[: 5])
"""
Index(['Belgium', 'Brazil', 'Canada', 'China', 'France'], dtype='object')
"""

上記の通り、インデックスに日付、カラムに国名(比較するアイテム)を持ったデータフレームがこのライブラリが想定しているデータ形式のようです。

早速、Bar Chart Raceを作ってみましょう。まずはメソッドに単純に渡して初期設定の出来を見てみます。

bcr.bar_chart_race(
    df=df,
    filename='covid19_default.mp4',
)

初期設定でも十分な見栄えのBar Chart Race ができましたね。

さらにこのライブラリは、見た目を変えたい場合に備えて、非常に多くのオプションが用意されています。
公式ドキュメントに色々設定を加えたバージョンも出てるのでそれも見ておきましょう。

bcr.bar_chart_race(
    df=df,
    filename='covid19_horiz.mp4',
    orientation='h',
    sort='desc',
    n_bars=6,
    fixed_order=False,
    fixed_max=True,
    steps_per_period=10,
    interpolate_period=False,
    label_bars=True,
    bar_size=.95,
    period_label={'x': .99, 'y': .25, 'ha': 'right', 'va': 'center'},
    period_fmt='%B %d, %Y',
    period_summary_func=lambda v, r: {'x': .99, 'y': .18,
                                      's': f'Total deaths: {v.nlargest(6).sum():,.0f}',
                                      'ha': 'right', 'size': 8, 'family': 'Courier New'},
    perpendicular_bar_func='median',
    period_length=500,
    figsize=(5, 3),
    dpi=144,
    cmap='dark12',
    title='COVID-19 Deaths by Country',
    title_size='',
    bar_label_size=7,
    tick_label_size=7,
    shared_fontdict={'family' : 'Helvetica', 'color' : '.1'},
    scale='linear',
    writer=None,
    fig=None,
    bar_kwargs={'alpha': .7},
    filter_column_colors=False)  

これの出力がこちらです。

引数名を見れば各引数が何を設定しているのかは大まかにわかるのではないかと思います。
また、公式ドキュメントのAPIリファレンスにも解説が充実しています。

指数分布について

前回の記事で名前だけ登場した指数分布についてついでに整理しておきます。
参考: 幾何分布の無記憶性について

指数分布は幾何分布の連続分布版のような確率分布です。
(古さに関係なく一定確率で壊れる機械について)機械が故障するまでの時間や、
(単位時間あたり一定確率で発生する災害について)災害が発生するまでの時間など、
一定確率で発生する何かしらの事象が、次に発生するまでの時間が従う分布です。

数学的には次のように定義されます。
パラメーター$\lambda > 0$に対して、次の確率密度関数を持つ分布を指数分布と呼び、$Exp(\lambda)$と書きます。
$$
f(x;\lambda) = \left\{
\begin{align}
&\lambda e^{-\lambda x} \quad & (x \geq 0)\\
&0 \quad & (x < 0)
\end{align}
\right.
$$
期待値は$\frac{1}{\lambda}$、分散は$\frac{1}{\lambda^2}$です。

モーメント母関数を使うと簡単に導出できますので見ておきましょう。
まず、モーメント母関数は$t<\lambda$の範囲で次のように定義されます。
$$
\begin{align}
M_X(t) &= E(e^{tX})\\
&= \int_{0}^{\infty} e^{tx}\lambda e^{-\lambda x} dx\\
&= \lambda \int_{0}^{\infty} e^{(t-\lambda)x} dx\\
&= \frac{\lambda}{\lambda -t}.
\end{align}
$$

これの微分は簡単ですね。1回微分と2回微分はそれぞれ次のようになります。
$$
\begin{align}
\frac{d}{dt}M_X(t) &= \frac{\lambda}{(\lambda-t)^2}\\
\frac{d^2}{dt^2}M_X(t) &= \frac{2\lambda}{(\lambda-t)^3}.
\end{align}
$$

これを使うと期待値と分散は次のように計算できます。
$$
\begin{align}
E(X) &= \left.\frac{d}{dt}M_X(t)\right|_{t=0}\\
&= \frac{\lambda}{(\lambda-0)^2}\\
&= \frac{1}{\lambda}.
\end{align}
$$
$$
\begin{align}
E(X^2) &= \left.\frac{d^2}{dt^2}M_X(t)\right|_{t=0}\\
&= \frac{2\lambda}{(\lambda-0)^3}\\
&= \frac{2}{\lambda^2}
\end{align}
$$
より、
$$
\begin{align}
V(X) &= E(X^2) – E(X)^2\\
&= \frac{2}{(\lambda-0)^2} – \left(\frac{1}{\lambda}\right)^2\\
&= \frac{1}{\lambda^2}.
\end{align}
$$

前回の記事でも触れました通り、指数分布は無記憶性を持つ連続分布です。
$x_1, x_2 \geq 0$に対して、$P(X\geq x_1+x_2|X\geq x_1) = P(X\geq x_2)$が成り立ちます。
実際、$x>0$とすると、 $P(X\geq x)=e^{-\lambda x}$ ですから、
$$
\begin{align}
P(X\geq x_1+x_2|X\geq x_1) &= \frac{P(X\geq x_1+x_2)}{P(X\geq x_1)}\\
&= \frac{e^{-\lambda (x_1+x_2)}}{e^{-\lambda x_1}}\\
&= e^{-\lambda x_2}\\
&= P(X\geq x_2)
\end{align}
$$
となります。

幾何分布の無記憶性について

ここ数回の記事で幾何分布に関連する話を取り上げているので、ついでに幾何分布が持つ無記憶性という性質について紹介します。
これは条件付き確率を用いて、次の数式で表される性質です。

$$
P(X > m+n|X > m) = P(X > ​n) \quad \text{ただし}m, n\geq 0.
$$

まず、幾何分布について上の数式が成り立つとを確認しておきましょう。
$P(X=k) = p(1-p)^{k-1}$ ですから、
$$
\begin{align}
P(X>n) &= \sum_{k=n+1}^{\infty}p(1-p)^{k-1}\\
&= p\cdot\frac{(1-p)^{n}}{1-(1-p)}\\
&= (1-p)^{n}
\end{align}
$$
となります。

よって、
$$
\begin{align}
P(X> m+n|X > m) &= \frac{P(X> m+n \land X > m)}{P(X > m)}\\
&= \frac{P(X > m+n)}{P(X > m)}\\
&= \frac{(1-p)^{m+n}}{(1-p)^{m}}\\
&= (1-p)^{n}\\
&= P(X>n)
\end{align}
$$
となり、幾何分布が冒頭の数式を満たすことが示されました。

これはどういうことか説明します。
幾何分布は確率$p$で成功する独立な試行を、初めて成功するまで繰り返すときに要した回数の分布ですから、
$P(X>n)$というのは、初めて成功するまでに$n+1$回以上かかる確率、言い換えると初めて成功するまでに$n$回以上失敗する確率になります。
これに対して、$P(X > m+n|X > m)$はどういうことかというと、成功するまでに$m+1$回以上かかる、つまりすでに$m$回失敗したという条件のもとで、
成功するのに$m+n+1$回以上かかる、つまり追加で$n$回以上失敗し成功するまでに$n+1$回以上かかる確率を意味します。

この二つが等しいということはどういうことかというと、
成功するまでに$n$回以上失敗する確率は、今の時点で何回失敗しているかという事実に全く影響を受けないということです。

例えば、$1/20$の確率で当たりが出るクジで、連続して20回ハズレを引くと、
そろそろ当たるんじゃないかなという気がしてくる人も多いと思うのですが、
そんなことは全くなく、この先あたりを引くまでにかかる回数の期待値は全く変わってないということを示しています。

この無記憶性は、離散分布の中では幾何分布だけが持つ性質です。
(逆にいうと、離散分布で、無記憶性を持っていたらそれは幾何分布だと言えます。)
このほか、連続分布まで範囲を広げると、指数分布が幾何分布同様に無記憶性を持ちます。

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

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

おさらいしておくと、$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

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