怪しいダイス(プロットしてみる)

ダイスを振った記録があるとして、その一部で、極端に偏りのあるダイス(この例では奇遇の偏り)が使われているのではという疑念を持ったときに、どうして調べればいいか。

前回 怪しいダイス(hmmlearn) は、隠れマルコフモデルで調べてるんだけど…というお話だったので、hmmlearnを使って、怪しい痕跡が見えるなーということを確かめました。

ただ、思ったようには収束しない場合があり、もうちょっとわかりやすい判別方法がありそうなので、今回は単純に二項分布を使って調べてみます。

ただし、ここで注意点として、このダイスを振った記録は、不正でありつつ、かつ全体では偏りを調節してあるということです。(単に偏りがあるだけなら普通に検定を行えばいい)

やっていることは単純で、

  1. ダイスを振った記録について、ダイスのそれぞれの面ごとに、ある面が出現していない繋がりを求める。

  2. それぞれの繋がりについて、確率を、二項分布 Bi(n, p) n=繋がりの長さ, p=5/6、確率密度関数に x=n として求める。 要するに、ダイスのある面が出る確率は1/6なので、ダイスのある面が出ない確率は5/6になり、それがn回続く確率を求めます。

  3. 必要に応じて、それぞれの面について計算した値を、ダイスを振った記録の位置ごとに足したりする。

ランダムに生成したデータと怪しいデータで若干違いが見えるようになりました。

import itertools
import numpy as np


np.random.seed(1)

nf = 6
s = np.random.randint(0, nf, 200)
ns = len(s)
state = [0 for _ in range(ns)]
result = [[0 for _ in range(ns)] for _ in range(nf)]

i = 0
while(True):
    v = s[i]
    if i > 0:
        for k in range(nf):
            if k == v:
                d = i - state[k]
                score = ss.binom(d, (nf-1.)/nf).pmf(d)
                for j in range(state[k], i):
                    result[k][j] += np.log(score)
                state[k] = i
    i += 1
    if ns == i:
        for k in range(nf):
            d = i - state[k]
            score = ss.binom(d, (nf-1.)/nf).pmf(d)
            for j in range(state[k], i):
                result[k][j] += np.log(score)
        break

%matplotlib inline
import matplotlib.pyplot as plt
plt.title("logP of faces (random generated data)")
for res in result:
    plt.plot(res)

f:id:ruby-U:20171001182048p:plain

import itertools
import numpy as np


np.random.seed(1)

nf = 6
s = np.random.randint(0, nf, 200)
ns = len(s)
state = [0 for _ in range(ns)]
result = [0 for _ in range(ns)]

i = 0
while(True):
    v = s[i]
    if i > 0:
        for k in range(nf):
            if k == v:
                d = i - state[k]
                score = ss.binom(d, (nf-1.)/nf).pmf(d)
                for j in range(state[k], i):
                    result[j] += np.log(score)
                state[k] = i
    i += 1
    if ns == i:
        for k in range(nf):
            d = i - state[k]
            score = ss.binom(d, (nf-1.)/nf).pmf(d)
            for j in range(state[k], i):
                result[j] += np.log(score)
        break
        
%matplotlib inline
import matplotlib.pyplot as plt
plt.title("logP (random generated data)")
plt.plot(result)

f:id:ruby-U:20171001182105p:plain

import itertools
import numpy as np


np.random.seed(1)

nf = 6
s = np.random.randint(0, nf, 200)
ns = len(s)
state = [0 for _ in range(ns)]
result = [[0 for _ in range(ns)] for _ in range(nf)]

i = 0
while(True):
    v = s[i]
    if i > 0:
        for k in range(nf):
            if k == v:
                d = i - state[k]
                score = ss.binom(d, (nf-1.)/nf).pmf(d)
                for j in range(state[k], i):
                    result[k][j] += np.log(score)
                state[k] = i
    i += 1
    if ns == i:
        for k in range(nf):
            d = i - state[k]
            score = ss.binom(d, (nf-1.)/nf).pmf(d)
            for j in range(state[k], i):
                result[k][j] += np.log(score)
        break
        
import numpy as np
result = [np.array(result[0])+np.array(result[2])+np.array(result[4]), 
          np.array(result[1])+np.array(result[3])+np.array(result[5])]

%matplotlib inline
import matplotlib.pyplot as plt
plt.title("logP of odd and even (random generated data)")
for res in result:
    plt.plot(res)

f:id:ruby-U:20171001182117p:plain

import itertools
import numpy as np


np.random.seed(1)

nf = 6
x = "12442552156235244611513665441521334522422516131461324351224626256625324354313241416144163515363423135135313531535155116446122436632362242642626642654256461421553232365354361461451264546335353161532646"
s = [int(s)-1 for s in x]
ns = len(s)
state = [0 for _ in range(ns)]
result = [[0 for _ in range(ns)] for _ in range(nf)]

i = 0
while(True):
    v = s[i]
    if i > 0:
        for k in range(nf):
            if k == v:
                d = i - state[k]
                score = ss.binom(d, (nf-1.)/nf).pmf(d)
                for j in range(state[k], i):
                    result[k][j] += np.log(score)
                state[k] = i
    i += 1
    if ns == i:
        for k in range(nf):
            d = i - state[k]
            score = ss.binom(d, (nf-1.)/nf).pmf(d)
            for j in range(state[k], i):
                result[k][j] += np.log(score)
        break

%matplotlib inline
import matplotlib.pyplot as plt
plt.title("logP of faces (fishy data)")
for res in result:
    plt.plot(res)

f:id:ruby-U:20171001182132p:plain

import itertools
import numpy as np


np.random.seed(1)

nf = 6
x = "12442552156235244611513665441521334522422516131461324351224626256625324354313241416144163515363423135135313531535155116446122436632362242642626642654256461421553232365354361461451264546335353161532646"
s = [int(s)-1 for s in x]
ns = len(s)
state = [0 for _ in range(ns)]
result = [0 for _ in range(ns)]

i = 0
while(True):
    v = s[i]
    if i > 0:
        for k in range(nf):
            if k == v:
                d = i - state[k]
                score = ss.binom(d, (nf-1.)/nf).pmf(d)
                for j in range(state[k], i):
                    result[j] += np.log(score)
                state[k] = i
    i += 1
    if ns == i:
        for k in range(nf):
            d = i - state[k]
            score = ss.binom(d, (nf-1.)/nf).pmf(d)
            for j in range(state[k], i):
                result[j] += np.log(score)
        break

%matplotlib inline
import matplotlib.pyplot as plt
plt.title("logP (fishy data)")
plt.plot(result)

f:id:ruby-U:20171001182142p:plain

import itertools
import numpy as np


np.random.seed(1)

nf = 6
x = "12442552156235244611513665441521334522422516131461324351224626256625324354313241416144163515363423135135313531535155116446122436632362242642626642654256461421553232365354361461451264546335353161532646"
s = [int(s)-1 for s in x]
ns = len(s)
state = [0 for _ in range(ns)]
result = [[0 for _ in range(ns)] for _ in range(nf)]

i = 0
while(True):
    v = s[i]
    if i > 0:
        for k in range(nf):
            if k == v:
                d = i - state[k]
                score = ss.binom(d, (nf-1.)/nf).pmf(d)
                for j in range(state[k], i):
                    result[k][j] += np.log(score)
                state[k] = i
    i += 1
    if ns == i:
        for k in range(nf):
            d = i - state[k]
            score = ss.binom(d, (nf-1.)/nf).pmf(d)
            for j in range(state[k], i):
                result[k][j] += np.log(score)
        break
        
import numpy as np
result = [np.array(result[0])+np.array(result[2])+np.array(result[4]), 
          np.array(result[1])+np.array(result[3])+np.array(result[5])]

%matplotlib inline
import matplotlib.pyplot as plt
plt.title("logP of odd and even (fishy data)")
for res in result:
    plt.plot(res)

f:id:ruby-U:20171001182154p:plain