#!/usr/bin/env python
# coding: utf-8

# 簡単な具体例について考えることにする。

# ## 例1

# 意見は有限でなく、$[0,1]$の間の値を一様な確率で取るとする。このとき、そうして得られた確率変数$x$に関して、
# 確率密度関数$f(x)$は
# 
# $$f(x) = 1\ \ \ (0\le x \le 1)\ \ \ \text{otherwise}\ \ 0$$
# 
# である。累積分布関数$F(x)$は
# 
# $$F(x) = x\ \ \ (0\le x \le 1)$$
# 
# となる。

# $i$が$k$番目に発言したあと、$j$が$k+1$番目に発言する確率が、発言者の人数を$n$人として、沈黙を含めない時、
# 
# $$p_{i}(k,j) = p(n) = \frac{1}{n}$$
# 
# のように、誰が発言する確率も等しい場合を考えることとする。

# $j$が発言するとき、その発言$x^{j}$は確率変数であり、その期待値と分散は
# 
# $$ E(x^{j}) = \frac{1}{2}$$
# $$ V(x^{j}) = \frac{1}{12}$$
# である。
# 

# 以上の議論とは関係なく、誰が発言しようと、その発言が一様な確率で$[0,1]$の間の値を取るとすると、この問題は結局ある閾値$r(0<r\le1)$を定めた時に、領域$[x^{k}-r, x^{k}+r]$の中に入るそれまでに出た点の数と、その個数がある値を超えるときまでに必要なステップ数はいくつか、という問題に帰着される。
# 
# 図にすると以下のような場面を考えていることになる。
# 
# <img src=simple001_1.jpg width=400 /> 
# 

# 一様な確率で$[0,1]$の間の数が選ばれるとき、その確率変数が$[\max(0,x-r), \min(x+r,1)]$の範囲に入っている確率は、確率密度関数を用いて、
# 
# 
# \begin{align}
# p(\max(0, x-r), \min(x+r, 1)) &= \int ^{\min(x+r,1)}_{\max(0, x-r)} 1 dx \\
# &= \left[ x\right]^{\min(x+r,1)}_{\max(0, x-r)}\\
# \end{align}
# 
# よって
# 
# $$p(x,r)= \left\{ \begin{array}{ll}x+r & 0\le x< \min(r,1-r) \\
# p(r) = \min(2r, 1) & \min(r, 1-r)\le x \le \max(1-r, r) \\
# 1 - x+r & \max(1-r, r) < x \le 1
# \end{array}\right.$$
# 
# を得る。これはグラフにすると下図のようになる。
# 

# In[1]:


get_ipython().run_line_magic('matplotlib', 'inline')
import matplotlib.pyplot as plt
from IPython.html.widgets import interactive
from IPython.display import display
import numpy as np

def plot_p(r=0.1):
    plt.figure(figsize=(10,8))
    x = np.array([0, min(r, 1-r), max(1-r,r), 1])
    y = np.array([r, min(2*r, 1), min(2*r, 1), r])
    plt.plot(x,y)
    _x = [-0.1,1.1]
    _y = [1,1]
    plt.plot(_x, _y, '--d')
    plt.xlabel(r'$x$')
    plt.ylabel(r'$p(x)$')
    plt.title(r'graph of $x$ - $p(x)$ when $r=%1.2f$' % r)
    plt.xlim(0,1)
    plt.ylim(0, 1.2)
    plt.show()

w = interactive(plot_p, r=(0.01, 1., 0.01))
display(w)


# これまで考えてきたように、会議の合意が得られたとする指標を、ある発言がなされた時に、それまでの発言のうちのいくつかに関連があることが条件であった。このリンクの結び方を、単純にその間の長さが一定値より小さければつなぐことにすると、これまでの考察をそのまま使うことができて、
# 
# $k+1$番目の発言$x_{k}$がなされた時、$k$番目までの発言のうち$y$個の発言が領域$[\max(0,x-r), \min(x+r,1)]$の中に存在する確率は、
# 
# $$_{k}C_{y}[p(r)]^{y}[p(r)]^{k-y}$$
# 
# で表せる。

# これは、
# 
# $$P[X=y] =\ _{k}C_{y}[p(r)]^{y}[1-p(r)]^{k-y}$$
# 
# のように書けば明らかなように、確率変数$X$に対するパラメータ$k,p$の二項分布$B(k,p)$を表している。ここで注意すべき点は、$p(x_{k}, r)$は$x_{k}$によって変わるものであったから、$x$に対して期待値を取ったものを考えなければならないことである。$p(r)$はそのようにして得られた期待値であることを意味している。
# 
# 
# この時、確率変数$X$に対する期待値と分散は、
# 
# $$E(X) = kp(r)$$
# 
# $$V(X) = kp(r)(1-p(r))$$
# 
# $X$の最頻値は、$(k+1)p$以下の最大の整数。ただし、$m = (k+1)p$において$m$が整数である場合、$m − 1$と$m$の双方が最頻値となる。

# また、時刻$k$が大きい時には二項分布は正規分布に近似できるので、期待値は中央値とほぼ等しくなる。

# 以上から、リンクの張り方に関する閾値$y$の選び方として
# 
# $$y = kp(x_{k+1}, r),$$
# $$\ ^{\forall} k\ \text{which satisfy}\ kp(x_{k+1}, r) > 5\\
# \text{and}\ kp(x_{k+1}, r)(1-p(x_{k+1},r)) > 5$$
# 
# というのが考えられる。このようにすれば、この条件を満たすどんな時刻$k$においても、またどのような閾値$r$を設定したとしても、時刻$k+1$で選ばれた$x_{k+1}$の$r$によって定まる領域において、$y$個以上の他の点を見出す確率を$1/2$とできる。

# すなわち、時刻の依存性はないため、$x_{k+1}$のまわりの領域に$y$個以上の他の点が存在する場合を1、存在しない場合を0とすれば、これらの起こる確率は、それぞれ$1/2$であるから、$k$回目に試行が終了する確率は$p=1/2$の[幾何分布](http://ja.wikipedia.org/wiki/%E5%B9%BE%E4%BD%95%E5%88%86%E5%B8%83)に従う:

# $$P[X=k] = \left( \frac{1}{2} \right)^{k}$$
# 
# このときの期待値と分散は
# 
# $$E(X) = \frac{1}{p} = 2$$
# $$V(X) = \frac{1-p}{p^{2}} = 2$$
# 
# すなわち、これまでのような場合を考えたとすると、平均として2回で試行が終了し、ほとんどの試行は1~3回で終了することになる。しかし、これは先程述べたような二項分布の正規分布への近似ができない領域であることに注意が必要である。したがって、実際には$y$の選び方をうまく選ぶことによって(たとえば試行が終了する回数の期待値が、二項分布を正規分布に近似できるような$k$となるように$p$を決めるような$y$)、実際に$k$の期待値はこの確率分布のものと等しくなるだろう。

# 以上の試行はあくまで閾値$y$を$y=kp(r)$とした場合のものであったが、実際の会議の場面を想定するなら、人数が多くなるほど結論に必要な時間は増大すると見るべきで、したがって今のような完全にランダムな発言の場合において$y$の選び方は、少なくとも今の選び方よりも大きい値となるようにしなければならない。

# また、二項分布$B(k,p)$のことを思い出せば、$y$の選び方をうまくとれば、任意の時刻$k$と閾値$r$について$p$の値が固定されるように$y$の値を決めることができ、その時の試行回数は先ほどの幾何分布に従う。

# 最後にこの例の条件と結果を整理する。
# 
# 条件:
# 
# - 沈黙なし
# - 発言者は等確率で選ぶ
# - 発言は[0,1]の一様乱数
# 
# 結果:
# 
# - 発言がなされたときにリンクが張られる数は二項分布$B(k,p(x_{k+1},r))$に従う
# - $y$を$k,r$にしたがって変化する適当な関数におけば、リンクの個数が$y$個を超えるのに必要な時間は幾何分布に従う
# 
# 

# ## 単純に試したもの

# In[58]:


import scipy.misc as sc

k = np.arange(2, 100)

def prob(r, y):
    # p(x,r) = p(r) = 2*rとする
    # 時系列に対して確率を求める
    p = []
    for _k in k:
        _p = 0.
        _y = int(y)
        while _y < _k+1:
            _p += sc.comb(_k, _y, exact=True)*((2*r)**_y)*((1-2*r)**(_k-_y))
            _y += 1
        p.append(_p)
    return p


# In[70]:


y = np.arange(1, 30)
r = 0.1
for _y in y:
    p = prob(r, int(_y))
    plt.plot(k, np.array(p), label=r"$y=%1.1f$" % _y)
plt.xlabel(r'$k$')
plt.ylabel(r'$p_{y}(r=%1.1f,k)$' % r)

plt.show()
    


# In[69]:


r = np.linspace(0.01, 0.49, 10)
y = 10
for _r in r:
    p = prob(_r, y)
    plt.plot(k, np.array(p), label=r"$r=%1.1f$" % _r)
plt.xlabel(r'$k$')
plt.ylabel(r'$p_{r}(y=%d, k)$' % y)

plt.show()


# In[64]:


y = 19
r = 0.1
p = prob(r, y)
a = np.array(p)<0
print np.where(a)
plt.plot(k, p)
plt.show()


# In[66]:


_p = 0.
_y = 19
_k = 77
r = 0.1
while _y < _k+1:
    _p += sc.comb(_k, _y, exact=True)*((2*r)**_y)*((1-2*r)**(_k-_y))
    _y += 1
print _p


# In[18]:


import itertools
a = itertools.combinations('ABCD', 2)
len(list(a))


# In[72]:


import scipy.misc as sc
sc.comb(77, 19)


# 組み合わせの大きさは膨大になりすぎ、小さい数に対して何十乗の計算を行っているために、誤差は無視できないほど大きくなる。大体$k=60$を超えると、信頼できる結果は得られそうにない。