始めに

の第4章を写経した記録をまとめます。前回はこちらPythonで体験するベイズ推論 PyMCによるMCMC入門第3章 (3.2全部) 写経。実行はJupyter Labにて実行しています。

　写経は以下のようにまとめました。

写経ではありますが、関数を作成し、コードの再実行をやりやすいようにした部分があります。
より良いと思われるコードを考えた場合は書き換えてコメントを添えるようにし、変更点をなるべく明示するようにしてあります。
個人的に気になった点のメモを残すようにしてあります。同じような疑問を持った方の助けになれば幸いです。
以前書いたコードと同じようなコード（例えばグラフの描写等）は効率化のために飛ばしているところもあります。
記事内で使用するモジュールなどは一番最初に宣言するようにしてあります。

　Pythonで体験するベイズ推論 PyMCによるMCMC入門を読んでいている際のちょっとした正誤表代わりになればと思います。

この記事で使用する主なモジュール、設定

　この記事では、主に以下のモジュールや設定を使用しています。import prawを実行する前に、pip install praw==3.6.0を実行してprawモジュールをインストールしておかないとエラーが出ます。

from os import makedirs
from urllib.request import urlretrieve
import praw
import pymc as pm
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt

4 偉大な定理、登場

4.3 サンプルサイズが小さいという災い

4.3.3 例題 : Redditコメントをソートする

makedirs("data", exist_ok=True)
urlretrieve("https://git.io/vXtX2",
            "top_showerthoughts_submissions.py")

Output:

    ('top_showerthoughts_submissions.py',
     <http.client.HTTPMessage at 0x1ce88b88f60>)

%run top_showerthoughts_submissions.py 2

print("Title of submission:")
print(top_post)  # top_postはtop_showerthoughts_submissions.pyを実行した際に代入される。型は文字列である。

Output:

    Title of submission:
    April Fools Day is proof the Purge wouldn’t work. Despite having a whole day for shitposting, there is still shitposting the rest of the year.

n_comments = len(contents)
comments = np.random.randint(n_comments, size=4)
print("Some Comments (out of %d total) \n----------"
      % n_comments)

for i in comments:
    print('"' + contents[i] + '"')
    print("upvotes/downvotes: ", votes[i, :], "\n")

Output:

    Some Comments (out of 2 total) 
    ----------
    "That guy who is always kicking robots in those Boston Dynamics videos is gonna have a really bad day when the machines take over."
    upvotes/downvotes:  [44  2] 
    
    "That guy who is always kicking robots in those Boston Dynamics videos is gonna have a really bad day when the machines take over."
    upvotes/downvotes:  [44  2] 
    
    "That guy who is always kicking robots in those Boston Dynamics videos is gonna have a really bad day when the machines take over."
    upvotes/downvotes:  [44  2] 
    
    "Shaving cream is like a progress bar of how much you've got left to shave."
    upvotes/downvotes:  [40  3]

メモ

　どうやら、まさかのコメントが2個しかないということが判明した。時間がたってしまって変化があったのかも・・・シンプルに comments = [0, 1] として進めます。

comments = [0, 1]
for i in comments:
    print('"' + contents[i] + '"')
    print("upvotes/downvotes: ", votes[i, :], "\n")

Output:

    "That guy who is always kicking robots in those Boston Dynamics videos is gonna have a really bad day when the machines take over."
    upvotes/downvotes:  [44  2] 
    
    "Shaving cream is like a progress bar of how much you've got left to shave."
    upvotes/downvotes:  [40  3]

def posterior_upvote_ratio(upvotes, downvotes, samples=20000):
    N = upvotes + downvotes
    upvote_ratio = pm.Uniform("upvote_ratio", 0, 1)
    observations = pm.Binomial("obs", N, upvote_ratio,
                               value=upvotes, observed=True)
    map_ = pm.MAP([upvote_ratio, observations]).fit()
    mcmc = pm.MCMC([upvote_ratio, observations])
    mcmc.sample(samples, samples / 4)
    return mcmc.trace("upvote_ratio")[:]

plt.rcParams['font.size'] = 17

figsize(22, 16)

colors = ["#348ABD", "#A60628", "#7A68A6", "#467821", "#CF4457"]

posteriors = []
for i in range(len(comments)):
    j = comments[i]
    label = f'({votes[j, 0]} up:{votes[j, 1]} down)\n{contents[j][:50]}...'
    posteriors.append(posterior_upvote_ratio(votes[j, 0],
                                             votes[j, 1]))
    plt.hist(posteriors[i], bins=18, density=True, alpha=.9,
             histtype="step", color=colors[i % 5], lw=3, label=label)
    plt.hist(posteriors[i], bins=18, density=True, alpha=.2,
             histtype="stepfilled", color=colors[i], lw=3)
    
plt.legend(loc="upper left")
plt.xlim(0, 1)
plt.ylabel("Density")
plt.xlabel("Probability of upvote")
plt.title("Posterior distributions of upvote ratios "
          "on differennt comments")

Output:

     [-----------------100%-----------------] 20000 of 20000 complete in 1.0 sec




    Text(0.5, 1.0, 'Posterior distributions of upvote ratios on differennt comments')

f:id:koheitsutsumi223:20190402144854p:plain

N = posteriors[0].shape[0]
lower_limits = []

for i in range(len(comments)):
    j = comments[i]
    label = f'({votes[j, 0]} up:{votes[j, 1]} down)\n{contents[j][:50]}...'
    plt.hist(posteriors[i], bins=20, density=True, alpha=.9,
             histtype="step", color=colors[i], lw=3, label=label)
    plt.hist(posteriors[i], bins=20, density=True, alpha=.2,
             histtype="stepfilled", color=colors[i], lw=3)
    v = np.sort(posteriors[i])[int(0.05 * N)]
    plt.vlines(v, 0, 10,
               color=colors[i], linestyles="--", linewidths=3)
    lower_limits.append(v)
    
plt.legend(loc="upper left")
plt.xlabel("Probability of upvote")
plt.ylabel("Density")
plt.title("Posterior distributions of upvote ratios "
          "on different comments")

order = np.argsort(-np.array(lower_limits))
print(order, lower_limits)

Output:

    [0 1] [0.87004917925665515, 0.82966075938026218]

f:id:koheitsutsumi223:20190402144858p:plain

4.3.5 でも計算が遅すぎる！

def intervals(u, d):
    a = 1. + u
    b = 1. + d
    mu = a / (a + b)
    std_err = 1.65 * np.sqrt((a * b) / ((a + b)**2 * (a + b + 1.)))
    return (mu, std_err)

posterior_mean, std_err = intervals(votes[:, 0], votes[:, 1])
lb = posterior_mean - std_err

print(f"Approximate lower bounds:\n{lb}, \n")

order = np.argsort(-lb)
ordered_contents = []
for i in order:
    ordered_contents.append(contents[i])
    print(votes[i, 0], votes[i, 1], contents[i])
    print("-----------------")

Output:

    Approximate lower bounds:
    [ 0.88044266  0.84187791], 
    
    44 2 That guy who is always kicking robots in those Boston Dynamics videos is gonna have a really bad day when the machines take over.
    -----------------
    40 3 Shaving cream is like a progress bar of how much you've got left to shave.
    -----------------

r_order = order[::-1]
plt.errorbar(posterior_mean[r_order],
             np.arange(len(r_order)),
             xerr=std_err[r_order],
             capsize=0, fmt="o", color="#7A68A6")
plt.xlim(0.3, 1)
plt.yticks(np.arange(len(r_order) - 1, -1, -1),
           map(lambda x: x[:30].replace("\n", ""),
               ordered_contents))

Output:

    ([<matplotlib.axis.YTick at 0x2730d10c0b8>,
      <matplotlib.axis.YTick at 0x2730caaaf28>],
     <a list of 2 Text yticklabel objects>)

f:id:koheitsutsumi223:20190402144903p:plain

メモ

map(lambda x: x[:30].replace("\n", ""),ordered_contents)はordered_contentsの各文字列の先頭30文字だけ抜き出し、改行を置換した要素を返すイテレータを生成します。

for i in map(lambda x: x[:30].replace("\n", ""),ordered_contents):
    print(i)

Output:

    That guy who is always kicking
    Shaving cream is like a progre

　一般的に、mapは第一引数の関数を第二引数の各値に適用させたイテレータを返す関数です。

for i in map(lambda x: 2*x, range(10)):
    print(i)

Output:

メモ

　今更だが、np.arangeの説明を読んでみた。rangeとほぼ同じような引数なんですね。

np.arange??（クリックすると展開されます）

Docstring:
arange([start,] stop[, step,], dtype=None)

Return evenly spaced values within a given interval.

Values are generated within the half-open interval ``[start, stop)``
(in other words, the interval including `start` but excluding `stop`).
For integer arguments the function is equivalent to the Python built-in
`range <http://docs.python.org/lib/built-in-funcs.html>`_ function,
but returns an ndarray rather than a list.

When using a non-integer step, such as 0.1, the results will often not
be consistent.  It is better to use ``linspace`` for these cases.

Parameters
----------
start : number, optional
    Start of interval.  The interval includes this value.  The default
    start value is 0.
stop : number
    End of interval.  The interval does not include this value, except
    in some cases where `step` is not an integer and floating point
    round-off affects the length of `out`.
step : number, optional
    Spacing between values.  For any output `out`, this is the distance
    between two adjacent values, ``out[i+1] - out[i]``.  The default
    step size is 1.  If `step` is specified, `start` must also be given.
dtype : dtype
    The type of the output array.  If `dtype` is not given, infer the data
    type from the other input arguments.

Returns
-------
arange : ndarray
    Array of evenly spaced values.

    For floating point arguments, the length of the result is
    ``ceil((stop - start)/step)``.  Because of floating point overflow,
    this rule may result in the last element of `out` being greater
    than `stop`.

See Also
--------
linspace : Evenly spaced numbers with careful handling of endpoints.
ogrid: Arrays of evenly spaced numbers in N-dimensions.
mgrid: Grid-shaped arrays of evenly spaced numbers in N-dimensions.

Examples
--------
>>> np.arange(3)
array([0, 1, 2])
>>> np.arange(3.0)
array([ 0.,  1.,  2.])
>>> np.arange(3,7)
array([3, 4, 5, 6])
>>> np.arange(3,7,2)
array([3, 5])
Type:      builtin_function_or_method

　今回の写経は以上です。ここまで読んでいただきありがとうございました。

一日一万字の感謝の写経

一日一万字の感謝の写経をして強くなります。そのうち本を置き去りにして何かを作り始める予定です。

Pythonで体験するベイズ推論 PyMCによるMCMC入門第4章写経

目次

始めに

この記事で使用する主なモジュール、設定

4 偉大な定理、登場

4.3 サンプルサイズが小さいという災い

4.3.3 例題 : Redditコメントをソートする

メモ

4.3.5 でも計算が遅すぎる！

メモ

メモ