モンティ・ホール問題で勘違いしていた話

論理パラドクスという本を読んでいたら、モンティ・ホール問題について勘違いしていたことに気付かされたのでメモ。
結論としては、普通のモンティ・ホール問題についてはこの勘違いでも(たまたま)正しいが、一般的な状況でも成り立つと思っているとまずいという話。正直、長い割にはしょうもない結論になった感が拭えない。

以前作った「意識高い」タグがとてつもなくバカっぽいが、数学関係なのでつけとく。「数理系」タグにしとけばよかった…。

勘違いの内容

モンティ・ホール問題で選択を変えることは「最初に選んだドア以外の2つのドアを選ぶ」ことと同じ意味だ。

確率を計算すると、最初の状態では当たりを選ぶ確率は1/3。
司会者が外れのドアを開けたあとは、ドアを変えた方が当たる確率は2/3になる。
たしかに最初に選んだドア以外の、2枚のドアを選んだときに当たる確率と同じになっている。

変形モンティ・ホール問題

なんとなく物語形式にしてみる。

モンティ・ホール問題で出てくるゲームをやることになった挑戦者。前日に猛勉強してコツ(a.k.a. 勘違い)をつかんだ。

「結局これは3つのドアのうち、1つを選ぶか、他の2つを選ぶかということだ。
司会者は必ず、挑戦者が選んだ1つのドア以外の2つのドアから、外れのドアを開ける。
ということは、選択を変えることは、他の2つのドアを選ぶのと同じだ。
1つより2つの方が当たる確率が高いのは当然だ。実際、確率は変えるときの方が変えないときの倍になっている。
簡単な話だ。選択を変えた方が勝ちだ!」

翌日、ゲームの説明を受ける挑戦者。
基本的にはモンティ・ホール問題なのだが、1つだけ違いがあった。

3つのドアA、B、Cのどれに賞品を入れるか、くじで決めるのだが、この確率が1/3ずつではなかった。
Aには\cfrac{9900}{10000}、Bには\cfrac{  99}{10000}、Cには\cfrac{   1}{10000}の確率で賞品を入れる。以降は同じで、司会者はどのドアに賞品が入ったか知っている。

こんなの絶対Aに決まってると思った挑戦者はあまりの興奮のため、誤ってBを選んでしまった。
だが選択を変えるチャンスが1回あるのはルールで決まっているのだ。そのときAに変えればいい…
しかし司会者が開けたのはAだった! Aは外れだったのだ。

挑戦者は思った。

「Aが外れだったのは驚きだな…奇跡的な確率じゃないか?
でもこれでBが当たりというのは確実だろう。BにはCの100倍近い確率で賞品が入るんだから。
…待てよ? 選択を変えるということは、残り2つのドアを選ぶということだったはずだ。
司会者がAを開けるまで、B以外が当たる確率は\cfrac{9900}{10000} + \cfrac{1}{10000}だった。これはBが当たる確率より100倍以上大きい。そしてそれがそっくり、Cが当たる確率になるはず…なのか?
だとすればBを選び続けるのはあまりに無謀すぎる。しかしCが当たりとはとても…」

当たる確率が高いのは選択を変える方なのだろうか、変えない方なのだろうか?

変形モンティ・ホール問題の定義

念のため、ちゃんと書いておく。この状況を変形モンティ・ホール問題と呼ぶことにする。

3つのドアA、B、Cのどれか1つに賞品が入っている(当たり)。残り2つのドアは外れ。
当たりのドアがどれか、挑戦者が当てるゲームをする。
このとき以下の条件がある。
1. 賞品をどのドアに入れるかはくじで決める。A、B、Cに賞品が入る確率をそれぞれP_A, P_B, P_C(0 \leq P_A, P_B, P_C \leq 1, P_A + P_B + P_C = 1)とする
2. 司会者はどのドアに賞品が入っているか知っているが、挑戦者は知らない
3. 挑戦者がドアを1つ選んだ後、司会者は残りの2つのドアのうち、外れのドアを(複数あれば無作為に選んで)必ず開け、選択肢から外す
4. 司会者が外れのドアを開けた後、挑戦者は必ず選択を変えるチャンスを与えられる

以上の事項を司会者、挑戦者は共通で認識している。両者の知識に違いがあるのは、どのドアが当たりか知っているかどうかだけ。

このとき、選択を変えて当たる確率はいくらか?
そしてそれが「最初に選んだドア以外の2つのドアを選ぶ」ときに当たる確率と等しいのはどんなときか?

解く

挑戦者は常に最初にBを選び、司会者は常にAを開けるとする。

確率を計算する

ドアAに賞品が入っている事象をA、ドアBに賞品が入っている事象をB、ドアCに賞品が入っている事象をCとする。司会者がドアAを開ける事象をXとする。
選択を変えて当たる確率はベイズの定理から

  \begin{aligned}  P(C|X) &= \cfrac{P(X|C)P(C)}{P(X)} \\         &= \cfrac{P(X|C)P(C)}{P(X|A)P(A) + P(X|B)P(B) + P(X|C)P(C)} \\         &= \cfrac{1 \cdot P_C}{0 \cdot P_A + \frac{1}{2} \cdot P_B + 1 \cdot P_C} \\         &= \cfrac{2P_C}{P_B + 2P_C}  \end{aligned}

なお、特にP_B = P_Cのときは常に2/3になる。

シミュレーション

Python3.6.3で動作確認。

# Copyright (c) 2017 sankaku
# Released under the MIT license
# http://opensource.org/licenses/mit-license.php

# Simulator of the Monty Hall problem.
# The names of the three doors are 'A', 'B' and 'C'.
# The challenger always chooses 'B'.
# The host always opens 'A'.

from random import random
from functools import reduce

# Probabilities of the prize being behind the door
P_A = 1/3
P_B = 1/3
P_C = 1 - P_A - P_B
# P_A = 9900/10000
# P_B =   99/10000
# P_C =    1/10000

# number of trials
TRIAL = 100000

# Return the door the prize is hidden behind
def get_answer():
    rand = random()
    if rand < P_A:
        return 'A'
    elif rand < (P_A + P_B):
        return 'B'
    else:
        return 'C'

def get_host_opened(answer):
    if answer == 'A':
        return 'C'
    elif answer == 'B':
        return 'A' if random() < 0.5 else 'C'
    else:
        return 'A'

def simulate_no_change():
    answer = 'A'
    while True:
        answer = get_answer()
        # Challenger chooses 'B', then ...
        if get_host_opened(answer) == 'A':
            break
    # Challenger does not change the choice
    return 'B' == answer

def simulate_change():
    answer = 'A'
    while True:
        answer = get_answer()
        # Challenger chooses 'B', then ...
        if get_host_opened(answer) == 'A':
            break
    # Challenger changes the choice to the other door
    return 'B' != answer


print('If the challenger does NOT change...')
num_hit_no_change = reduce(lambda x, y: x + y, [simulate_no_change() for _ in range(0, TRIAL)], 0)
print('Hit: {0}/{1} = {2}%\n'.format(num_hit_no_change, TRIAL, 100 * num_hit_no_change/TRIAL))

print('If the challenger DOES change...')
num_hit_change = reduce(lambda x, y: x + y, [simulate_change() for _ in range(0, TRIAL)], 0)
print('Hit: {0}/{1} = {2}%'.format(num_hit_change, TRIAL, 100 * num_hit_change/TRIAL))

検証

あまり自信がないので、普通のモンティ・ホール問題のときに、確率とシミュレーションが一致していることを確かめる。

このときP_A = P_B = P_C = \cfrac{1}{3}で、
変えて当たる確率はP(C|X) = \cfrac{2P_C}{P_B + 2P_C} = \cfrac{2}{3}となる。

シミュレーションでは以下。

If the challenger does NOT change...
Hit: 33263/100000 = 33.263%

If the challenger DOES change...
Hit: 66367/100000 = 66.367%

変えて当たる確率は66.367%という結果。だいたい合ってる。

本題

冒頭の物語での答え

P_A = \cfrac{9900}{10000}, P_B = \cfrac{99}{10000}, P_C = \cfrac{1}{10000}の場合を考える。
変えて当たる確率はP(C|X) = \cfrac{2P_C}{P_B + 2P_C} = \cfrac{2}{101}になる。
シミュレーションでも同様の結果。

If the challenger does NOT change...
Hit: 98048/100000 = 98.048%

If the challenger DOES change...
Hit: 2054/100000 = 2.054%

よって、変えない方が当たる。

…実はこれは計算しなくてもすぐわかることだった。
いまはP_Cが非常に小さいときを考えているが、極端な場合として0のときにはCは確実に外れだから、選択を変えない方が当たる。

「最初に選んだドア以外の2つのドアを選ぶ」のと等価な条件

これはP(C|X) = P(A) + P(C)ということだ。左辺は既に求めたので代入するとP_A = P_CまたはP_B = 0となる。これが答え。
2番目の場合は当然でつまらないので無視。
1番目の場合はAとCに賞品が入る確率が等しいということ。つまりAもCも同等。だから司会者がそのどちらを開けても、何も情報が増減しない。
ということは最初に選んだドアが当たる確率は前と同じ。残った1枚のドアが当たる確率は、選んでいない2枚のドアが当たる確率の合計値になるはずだ。
そう考えると納得はいく。

普通のモンティ・ホール問題のときはP_A = P_Cが成り立っていたので、たまたま冒頭のような勘違いをしても正解だった。
変形モンティ・ホール問題ではこれが成り立つとは限らないので、冒頭の物語の挑戦者のように(自分のように)考えるとおかしな結果になってしまうという話だった。

補足

よくモンティ・ホール問題のわかりやすい解説で、ドアを大量に増やすというものがある。
たとえばドアを10000枚にして、挑戦者は1番目のドアを選ぶ。司会者は3172番目以外のすべてのドアを開けて、
1番目のドアのまま変えないか、3172番目のドアに変えるかというチャンスを与える。
こう考えると、選択を変える方が当たる確率が高いのは明らかだという説明だ。最初に挑戦者が選んだドアがそのまま当たる確率は非常に低く、残ったドアに当たりがあると考えてまず間違いない。

今回の変形モンティ・ホール問題はこれに似ているように見える。
9998枚のドアをグループA、1枚のドアをグループB、残りの1枚のドアをグループCと名付けて、
どのグループに賞品が入っているかという問題にすればP_A = \cfrac{9998}{10000}, P_B = P_C = \cfrac{1}{10000}での変形モンティ・ホール問題になる。
このとき、選択を変えた方が当たる確率は2/3にしかならない…?

なぜこうなるのかというと、ドアを増やす解説で考えている状況はここでの変形モンティ・ホール問題とはまったく別物だからだ。
ドアを増やす場合、グループはあらかじめ決まっておらず、司会者がドアを開けて初めて定義できる。変形モンティ・ホール問題をグループ版にした場合は、あらかじめどのドアがどのグループかが決まっている。だからそのまま当てはめることはできない。

pythonで二重振り子のアニメーションを作る

QiitaのPythonタグをなんとなく見ているけれど、機械学習や数学関連の、意識高い投稿が多くて羨ましい。というわけで、今回は意識高く、二重振り子のアニメーションを作ってみる。
…と思って書いていたのだが、もうちょっと調べてから上げようと思っていたら2ヶ月以上経ってしまった。
このままだと永久に公開できない気がしたので、まだちゃんとできていないところがいろいろあるがとりあえず上げておく。
あとで手を加えていく…予定。

環境

  • ArchLinux
    • Python 3.6.1
    • pip 9.0.1
    • matplotlib 2.0.0
    • numpy 1.12.1
    • PyQt 5.8.2
    • scipy 0.19.0
  • Ubuntu 16.04
    • Python 3.5.2
    • 他は同じ

インストール方法

$ sudo pacman -S python-pip
$ sudo pip install matplotlib scipy numpy pyqt5

Ubuntuなら

$ sudo apt-get install python3-pip
$ pip3 install matplotlib numpy scipy sympy

でいけたはず。Ubuntuだとpipで落とすパッケージは ~/.local/lib/python3.5/site-packages に入れてるっぽいのでsudoが不要なのか。あといちいちpip3とかpython3とか打たないとダメなのが面倒。デフォルトがpython2のため。

matplotlibの設定

ArchLinuxだとそのままでは描画できなかった。

/usr/lib/python3.6/site-packages/matplotlib/mpl-data/matplotlibrc のbackendを qt5aggに変更。

Python で matplotlib を試そうとしたら ‘cairo.Context’ ってナニソレ – Qiita参照。

matplotlibのテスト: sin

sin曲線を描画してみる。

import numpy as np
import matplotlib.pyplot as plt

PI = np.pi
x = np.arange(-PI, PI, 0.1)
y = np.sin(x)
plt.grid()
plt.plot(x, y)
plt.show()

numpyパッケージをnpという名前でimportしている。
numpy.arange (arrangeではない)はNumPy配列を返し、 numpy.sin は配列の各要素に対してsinの値を返す。

sin

  • NumPy配列

    numpy.ndarray のことをNumPy配列と呼ぶらしい。これはN-Dimension Arrayだからndarrayのよう。
    面白いのは y = np.sin(x)xy もNumPy配列になっていること。つまりxの各要素のsinを取ったものがy。
    なおPythonでは普通の配列はlist。型は type でわかる。

    >>> import numpy as np
    >>> a = np.array([1, 2, 3])
    >>> a
    array([1, 2, 3])
    >>> type(a)
    <class 'numpy.ndarray'>
    >>> b = [1, 2, 3]
    >>> b
    [1, 2, 3]
    >>> type(b)
    <class 'list'>
    
  • numpy.arange([start, ]stop[, step])

    startからstopをstepずつ増やしたNumPy配列を返す。numpy.arange — NumPy v1.12 Manual参照。

    >>> import numpy as np
    >>> np.arange(5)
    array([0, 1, 2, 3, 4])
    >>> np.arange(5, 9)
    array([5, 6, 7, 8])
    >>> np.arange(5, 9, 2)
    array([5, 7])
    
  • pyplot.grid()

    グリッドを描く。pyplot — Matplotlib 2.0.2 documentation参照。
    色とか線の形などを変えるときは以下のようにする。

    plt.grid(color='r', linestyle='dotted', linewidth=1)
    
  • pyplot.plot()

    グラフを描く。pyplot — Matplotlib 2.0.2 documentation参照。
    色とか線の形などを変えるときはgridと同様。

    plt.plot(x, y, color='r', linestyle='dotted', linewidth=3)
    

    これと同じことを略記して

    plt.plot(x, y, 'r--', linewidth=3)
    

    とも書ける。

  • numpy.sin, numpy.cos, numpy.pi

    はじめに

    from numpy import sin, cos, pi
    

    と書いておけば、いちいち np.sin とか書かなくても sin だけで書ける。

比率の固定

上で描画したsin曲線のグラフは、ウィンドウのサイズを変えると縦長になったり横長になったりして比率が変わる。それを固定化する。
図の中にもう一つ図(サブプロット)を入れて、サブプロットにいろいろ設定して描画する感じ。

import numpy as np
import matplotlib.pyplot as plt

PI = np.pi
x = np.arange(-PI, PI, 0.1)
y = np.sin(x)

fig = plt.figure()
ax = fig.add_subplot(111, aspect='equal', autoscale_on=False, xlim=(-PI, PI), ylim=(-PI, PI))
ax.grid()

plt.plot(x, y)
plt.show()
  • figure.Figure.add_subplot

    figure — Matplotlib 2.0.2 documentation参照。

    これを使うと複数のグラフを並べて描ける。
    最初の’111’というのはグラフの位置指定。たとえばこれが’234’なら、ウィンドウを縦2横3に分割して、その(左上から右、そして下に数えて)4番目を使うということ。
    以下のような感じ。

    import numpy as np
    import matplotlib.pyplot as plt
    
    PI = np.pi
    
    fig = plt.figure()
    
    ax = fig.add_subplot(234, aspect='equal', autoscale_on=False, xlim=(-PI, PI), ylim=(-PI, PI))
    x1 = np.arange(-PI, PI, 0.1)
    y1 = np.sin(x1)
    ax.grid()
    plt.plot(x1, y1)
    
    ax = fig.add_subplot(232, aspect='equal', autoscale_on=False, xlim=(-PI, PI), ylim=(-PI, PI))
    x2 = np.arange(-PI, PI, 0.1)
    y2 = np.cos(x2)
    ax.grid(color='b', linestyle='dotted', linewidth=0.1)
    plt.plot(x2, y2, 'r--')
    
    plt.show()
    

    sin_cos

微分方程式の数値計算: 自由落下

dy/dt = f(t, y)という形の常微分方程式を解くのはscipy.integrate.odeintでできる。scipy.integrate.odeint — SciPy v0.18.1 Reference Guide 参照。
FORTRANのodepackにあるlsodaという関数(?)を使って数値計算しているそう。

たとえば自由落下はy軸を鉛直上向きに取って \cfrac{d^{2}y}{dt^2} = -g なので、これを1階に分解して
  \left\{  \begin{array}{l}  \cfrac{dy}{dt} = v\\  \cfrac{dv}{dt} = -g  \end{array}  \right.

import numpy as np
from scipy.integrate import odeint

def ode(f, t):
    y, v = f
    dfdt = [v, -G]
    return dfdt

G  = 9.8        # [m/s^2] gravitational acceleration
V0 = 0          # [m/s]   initial velocity
DURATION = 10   # [s] duration time
INTERVAL = 1    # [s] interval time

f0 = [0, V0]    # [y, v] at t = 0
t = np.linspace(0, DURATION, DURATION + INTERVAL)

sol = odeint(ode, f0, t)
print(sol)

f0には[y, v]の初期値を入れている。ode(f, t)はただの連立微分方程式で、ここに上に書いた自由落下の式を入れている。

結果は以下。
[[y, v] at t = 0,
[y, v] at t = 1,
…,
[y, v] at t = 10]
ってこと。

[[   0.     0. ]
 [  -4.9   -9.8]
 [ -19.6  -19.6]
 [ -44.1  -29.4]
 [ -78.4  -39.2]
 [-122.5  -49. ]
 [-176.4  -58.8]
 [-240.1  -68.6]
 [-313.6  -78.4]
 [-396.9  -88.2]
 [-490.   -98. ]]
  • numpy.linspace(start, stop, num)

    startからstopまでをnum等分した配列を返す。 numpy.arange とよく似てる。こんなメソッドもあるということで使ってみた。
    numpy.linspace — NumPy v1.12 Manualを参照。

    >>> import numpy as np
    >>> np.linspace(1, 3, 5)
    array([ 1. ,  1.5,  2. ,  2.5,  3. ])
    
  • y, v = f

    Pythonではこんな書き方ができるらしい。いらない要素は _ に入れればいいっぽい。

    >>> x = [1, 2, 3]
    >>> a, _, c = x
    >>> a
    1
    >>> c
    3
    

数値計算の結果をアニメーションにする

上の結果をアニメーションにする。自由落下だと退屈なので、ちょっと上向きに放ってみる。
なお、最後の2行はアニメーションをmp4/gifで保存するためのもの。いらなければコメントアウトしていい。

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.animation import FuncAnimation
from scipy.integrate import odeint

def ode(f, t):
    y, v = f
    dfdt = [v, -G]
    return dfdt

G  = 9.8        # [m/s^2] gravitational acceleration
V0 = 20         # [m/s]   initial velocity
DURATION = 10   # [s] duration time
INTERVAL = 0.1  # [s] interval time

f0 = [0, V0]    # [y, v] at t = 0
t = np.arange(0, DURATION + INTERVAL, INTERVAL)

sol = odeint(ode, f0, t)



fig = plt.figure()
ax = fig.add_subplot(111, aspect='equal', autoscale_on=False, xlim=(-100, 100), ylim=(-500, 100))
ax.grid()

ball, = plt.plot([], [], 'ro', animated=True)

def init():
    return ball,

def update(y):
    ball.set_data(0, y)
    return ball,

FRAME_INTERVAL = 1000 * INTERVAL # [msec] interval between frames
FPS = 1000/FRAME_INTERVAL # frames per second

ani = FuncAnimation(fig, update, frames=sol[:, 0],
                    init_func=init, blit=True, interval=FRAME_INTERVAL, repeat=True)
plt.show()
ani.save('fall.mp4', fps=FPS, extra_args=['-vcodec', 'libx264'])
ani.save('fall.gif', writer='imagemagick', fps=FPS)
  • set_data()

    lines — Matplotlib 2.0.2 documentation参照。
    座標を指定すればその通りに描画される。

  • ball, = plt.plot([], [], ‘ro’, animated=True)

    この , はどう見ても間違ってると思ったのだが…。
    これは plt.plot([], [], 'ro', animated=True) がlist(要素数1)で、 ball にそのlistの中身を入れるためにやるよう。
    animation example code: simple_anim.py — Matplotlib 2.0.2 documentationでもやってるので同じようにやった。
    plt.plot([], [], 'ro', animated=True) で2つ[]を入れているのは自信なし。

    >>> a, = [999]
    >>> a
    999
    
  • sol[:, 0]

    NumPy多重配列のslice。solは多重配列だが、要素の配列の0番目を集めた配列ということ。
    Python (numpy) で出てくるコロン,カンマ [:,]参照。

    >>> import numpy as np
    >>> a = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
    >>> a
    array([[1, 2, 3],
           [4, 5, 6],
           [7, 8, 9]])
    >>> a[:, 0]
    array([1, 4, 7])
    

単振り子のアニメーション

二重振り子ではどうせ解析力学を使うので、これも同じようにやってやる。

  \mathcal{L} = \cfrac{1}{2}ml^2\dot{\theta}^2 + mgl\cos\theta
なのでEuler-Lagrange eq.(意識高いので英語で書く)は
  \ddot{\theta} = - \cfrac{g}{l} \sin\theta

これを数値計算するために1階に分解してやると
  \left\{  \begin{array}{l}  \cfrac{d\theta}{dt} = \dot{\theta}\\  \cfrac{d\dot{\theta}}{dt} = -\cfrac{g}{l}\sin\theta  \end{array}  \right.
のような形になる。あとは落下のときと同様に描画してやればいい。
経過時間を表示するようにしてみた。

import matplotlib.pyplot as plt
import numpy as np
from numpy import sin, cos, pi
from matplotlib.animation import FuncAnimation
from scipy.integrate import odeint

def ode(f, t):
    theta, dtheta = f
    dfdt = [dtheta, -(G/L) * sin(theta)]
    return dfdt

G  = 9.8         # [m/s^2] gravitational acceleration
THETA0 = pi/4    # [rad]   initial angle
V0 = 1           # [m/s]   initial velocity
L = 1            # [m]     length of the pendulum
DURATION = 10    # [s]     duration time
INTERVAL = 0.05  # [s]     interval time

f0 = [THETA0, V0/L]    # [theta, v] at t = 0
t = np.arange(0, DURATION + INTERVAL, INTERVAL)

sol = odeint(ode, f0, t)

fig = plt.figure()
ax = fig.add_subplot(111, aspect='equal', autoscale_on=False, xlim=(-L, L), ylim=(-L, L))
ax.grid()

theta = sol[:, 0]
x = L * sin(theta)
y = - L * cos(theta)
line, = plt.plot([], [], 'ro-', animated = True)

time_template = 'time = %.1fs'
time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)

def init():
    time_text.set_text('')
    return line, time_text

def update(i):
    next_x = [0, x[i]]
    next_y = [0, y[i]]
    line.set_data(next_x, next_y)

    time_text.set_text(time_template % (i*INTERVAL))
    return line, time_text

FRAME_INTERVAL = 1000 * INTERVAL # [msec] interval between frames
FPS = 1000/FRAME_INTERVAL # frames per second

ani = FuncAnimation(fig, update, frames=np.arange(0, len(t)),
                    interval=FRAME_INTERVAL, init_func=init, blit=True)
plt.show()
ani.save('single_pendulum.mp4', fps=FPS, extra_args=['-vcodec', 'libx264'])
ani.save('single_pendulum.gif', writer='imagemagick', fps=FPS)

二重振り子のアニメーション

数式を書いてるが、正しさは一切責任持たない。typoあると思う。
角度や長さなどの記号は二重振り子 – Wikipediaと同じ。

Lagrangianは
  \mathcal{L} = \cfrac{1}{2} (m_{1} + m_{2}) l_{1}^{2} \dot{\theta_{1}}^{2} + \cfrac{1}{2} m_{2} l_{2}^{2} \dot{\theta_{2}}^{2} + m_{2} l_{1} l_{2} \dot{\theta_{1}} \dot{\theta_{2}} \cos(\theta_{2} - \theta_{1}) + (m_{1} + m_{2}) gl_{1}\cos\theta_{1} + m_{2}gl_{2}\cos\theta_{2}

あとはE-L eq.を求めて、それを1階の連立微分方程式に分解すればいい。
面倒な計算をすると
  \left\{  \begin{array}{l}  \cfrac{d\theta_{1}}{dt} = \dot{\theta_{1}}\\  \cfrac{d\theta_{2}}{dt} = \dot{\theta_{2}}\\  \cfrac{d\dot{\theta_{1}}}{dt} = \cfrac{1}{l_{1}l_{2} \left( m_{1} + m_{2} \sin^2 (\theta_{2} - \theta_{1})\right)} l_{2} \left( m_{2} g \sin\theta_{2} \cos(\theta_{2} - \theta_{1}) + m_{2} l_{1} \dot{\theta_{1}}^2 \sin(\theta_{2} - \theta_{1}) \cos(\theta_{2} - \theta_{1}) + m_{2} l_{2} \dot{\theta_{2}}^2 \sin(\theta_{2} - \theta_{1}) - (m_{1} + m_{2}) g \sin\theta_{1}\right)\\  \cfrac{d\dot{\theta_{2}}}{dt} = \cfrac{1}{l_{1}l_{2} \left( m_{1} + m_{2} \sin^2 (\theta_{2} - \theta_{1})\right)} l_{1} \left( -(m_{1} + m_{2}) l_{1} \dot{\theta_{1}}^2 \sin(\theta_{2} - \theta_{1}) - (m_{1} + m_{2}) g \sin\theta_{2} - m_{2}l_{2}\dot{\theta_{2}}^2\sin(\theta_{2} - \theta_{1}) \cos(\theta_{2} - \theta_{1}) + (m_{1} + m_{2})g\sin\theta_{1} \cos(\theta_{2} - \theta_{1}) \right)  \end{array}  \right.
を(たぶん)得る。

あとは同じようにやればいい。
単振り子と違うのは、質点が2個あること。これは単振り子のときに f = [\theta, \dot{\theta}] としていたのを f = [\theta_{1}, \dot{\theta_{1}}, \theta_{2}, \dot{\theta_{2}}] とすればいい。

from numpy import sin, cos, pi
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import numpy as np
from scipy.integrate import odeint

def ode(f, t):
    theta1, dtheta1, theta2, dtheta2 = f
    dtheta1dt = (M2 * G * sin(theta2) * cos(theta2 - theta1)
            + M2 * L1 * (dtheta1 ** 2) * sin(theta2 - theta1) * cos(theta2 - theta1)
            + M2 * L2 * (dtheta2 ** 2) * sin(theta2 - theta1)
            - (M1 + M2) * G * sin(theta1)) / (M1 * L1 + M2 * L1 * (sin(theta2 - theta1) ** 2))
    dtheta2dt = (-(M1 + M2) * L1 * (dtheta1 ** 2) * sin(theta2 - theta1)
            - (M1 + M2) * G * sin(theta2)
            - M2 * L2 * (dtheta2 ** 2) * sin(theta2 - theta1) * cos(theta2 - theta1)
            + (M1 + M2) * G * sin(theta1) * cos(theta2 - theta1)) / (M1 * L2 + M2 * L2 * (sin(theta2 - theta1) ** 2))
    return [dtheta1, dtheta1dt, dtheta2, dtheta2dt]

G  = 9.8            # [m/s^2]  gravitational acceleration
THETA1_0 =  pi      # [radian] initial theta
THETA2_0 =  pi/6    # [radian] initial theta
V1_0 = 0            # [m/s]    initial velocity
V2_0 = 0            # [m/s]    initial velocity
L1 =  1             # [m]      length of pendulum
L2 =  1             # [m]      length of pendulum
M1 =  1             # [kg]     mass
M2 =  1             # [kg]     mass
DURATION = 10       # [s]      duration time
INTERVAL = 0.05     # [s]      interval time

f_0 = [THETA1_0, V1_0/L1, THETA2_0, V2_0/L2]
t = np.arange(0, DURATION + INTERVAL, INTERVAL)

sol = odeint(ode, f_0, t)
theta1, theta2 = sol[:, 0], sol[:, 2]

x1 = L1 * sin(theta1)
y1 = - L1 * cos(theta1)
x2 = x1 + L2 * sin(theta2)
y2 = y1 - L2 * cos(theta2)

fig = plt.figure()
ax = fig.add_subplot(111, aspect='equal', autoscale_on=False, xlim = [-L1 - L2, L1 + L2], ylim = [-L1 - L2, L1 + L2])
ax.grid()

line, = plt.plot([], [], 'ro-', animated = True)

time_template = 'time = %.1fs'
time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)

def init():
    time_text.set_text('')
    return line, time_text

def update(i):
    next_x = [0, x1[i], x2[i]]
    next_y = [0, y1[i], y2[i]]
    line.set_data(next_x, next_y)

    time_text.set_text(time_template % (i * INTERVAL))
    return line, time_text

FRAME_INTERVAL = 1000 * INTERVAL # [msec] interval time between frames
FPS = 1000/FRAME_INTERVAL        # frames per second
ani = FuncAnimation(fig, update, frames=np.arange(0, len(t)),
                    interval=FRAME_INTERVAL, init_func=init, blit=True)

plt.show()
ani.save('double_pendulum.mp4', fps=FPS, extra_args=['-vcodec', 'libx264'])
ani.save('double_pendulum.gif', writer='imagemagick', fps=FPS)

感想

Qiitaとかに数式ゴリゴリ書いている人たちの情熱はすごい。もう数式とか書きたくない。
あと複雑系はけっこう好きなので、こうやって数値計算できるのは楽しい。

参考