ブロックごとの分布から合計点の分布

　notebookはこちら。

　次のような試験があるとする。

・試験はT(0)～T(n_tests - 1)の(n_tests)個のブロックで構成されている。

・すべての問題は、易しめのA問題と、難しめのB問題に分けられる。
・配点はすべて1点である。

・T(0)はすべてA問題のみで構成されている。

・T(i) (i≧0)は n_questions(i)個の問題で構成され、A問題 n_a(i)個とB問題 n_b(i)個で構成されている。(n_quiestions(i)=n_a(i)+n_b(i))
・受験者数は n_subjects人である。

・j番目の受験者はA問題に正答率 p_a(j)で、B問題に正答率 p_b(j)で正答する。

・正答率 p_a(j)、p_b(j)の間には、全受験者に共通のkを用いて p_b(j)=k*p_a(j) の関係が成り立つ。(0<k<1)

・試験後に受験者に返却される成績表には、各ブロックの自分の得点率(1%刻み)と、各ブロックで得点率区間(5%刻み)に何人の受験者がいるかのヒストグラムが掲載されている。

　この設定下で、合計点の分布および自分の推定順位を求めたい。

　そこで次のような処理を考えた。

(1) T(0)の成績は受験者の正答率 p_a(j) に忠実だと仮定し、p_a(j)を求める。この時、T(0)の得点率ヒストグラムにおいて各得点率区間には一様に受験者が分布するものと仮定する。(つまり55-60%の区間に4人いた場合、その4人のp_aは[55.00, 56.25, 57.50, 58.75]ということにする。)

(2) 適当なkについて、p_b(j) を求める。

(3) n_a(i)とn_b(i)を求める。

　 上で仮定したp_a(j)、p_b(j)をもつ受験者がT(i)を受けた場合の得点率の確率変数を、二項分布 B(n_a(i), p_a(j)) と B(n_b(i), p_b(j)) から求める。これを全受験者について求め合計し、実際の結果のヒストグラムとのKL情報量を求める。この作業を考えられるすべてのn_a(i), n_b(i)の組み合わせで行い、その中でKL情報量が最も小さくなるn_a(i)、n_b(i)を採用する。

(4) (2)と(3)を繰り返し、納得のいくkを求める。

(5) (4)で決定したkを用いて、受験者のp_a(j)、p_b(j)を決定する。

(6) (5)で決定したp_a(j)、p_b(j)を用いて、各受験者の合計点の確率変数を求める。

(7) 全受験者の合計点の確率変数の総和をとる。

　上記を実装する。使用するデータは、n_tests=5、n_questions=[100, 50, 30, 30, 40]、n_subjects=100としてテキトーに作成したsample_distribution.csvである。sample_distribution.csvの (x, y) には、T(x)の5*y～5*(y+1)%にいる人数が入っている。

　まずはk=0.4, 0.5, ... , 0.8 の各場合について、ブロックごとのA, B問題の配分およびKL情報量を表示する。

------------------------------------------------------------

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import math

def Combination(n, r):
    return math.factorial(n)//math.factorial(n-r)//math.factorial(r)

def Binomial(param):
    N, p=param
    return [Combination(N, k)*pow(p, k)*pow(1-p, N-k) for k in range(N+1)]

def Binomials(list_n_p):
    #list_n_p: 2-dimensioin vector (shape=(n, 2))
    #list_n_p[i]: i-th setting for Binomial distribution.
    #list_n_p[i][0]: n of i-th Binomial dist.
    #list_n_p[i][1]: p of i-th Binomial dist.
   
    list_n_p=[list_n_p[i] for i in range(len(list_n_p)) if list_n_p[i][0]!=0]
    hist=Binomial(list_n_p[0])
    for i in range(len(list_n_p)-1):
        pre_prob=hist
        add_prob=Binomial(list_n_p[i+1])
        l=np.max([len(pre_prob), len(add_prob)])
        mx=len(pre_prob)+len(add_prob)-2

        #probability map
        prob_map=np.zeros((l, l))
        prob_map[:len(pre_prob), :len(add_prob)]=np.dot(np.array(pre_prob).reshape(-1, 1), np.array(add_prob).reshape(1, -1))

        #upper-left part of matrix
        hist_former=[np.sum([prob_map[i, s-i] for i in range(s+1)]) for s in range(l)]

        #lower-right part of matrix
        hist_latter=[np.sum([prob_map[-1-i, l-(minus-i)-1] for i in range(minus+1)]) for minus in range(l-1)]
        hist=hist_former+hist_latter[::-1]
        hist=hist[:mx+1]
       
    return hist

def Sum_Prob(list_A, list_B):
    l=np.max([len(list_A), len(list_B)])
    mx=len(list_A)+len(list_B)
   
    prob_map=np.zeros((l, l))
    prob_map[:len(list_A), :len(list_B)]=np.dot(np.array(list_A).reshape(-1, 1), np.array(list_B).reshape(1, -1))
   
    hist_former=[np.sum([prob_map[i, s-i] for i in range(s+1)]) for s in range(l)]
    hist_latter=[np.sum([prob_map[-1-i, l-(minus-i)-1] for i in range(minus+1)]) for minus in range(l-1)]
   
    hist=hist_former+hist_latter[::-1]
    hist=hist[:mx+1]
   
    return hist

class test:
    def __init__(self, n_total, distribution):
        self.n_total=n_total
        self.distribution=distribution
        self.n_a=np.nan
        self.n_b=np.nan
        self.KLd=np.nan
        self.edit_distribution()
        return
   
    def edit_distribution(self):
        total=np.sum(self.distribution)
        self.distribution=[self.distribution[i]/total for i in range(len(self.distribution))]
   
    def calc_n(self, p_a, p_b, detail=False):
        temp_n_a=0
        while temp_n_a<=self.n_total:
            temp_n_b=round(self.n_total-temp_n_a)
            temp_dist=np.zeros(self.n_total+1)
            for i in range(p_a.shape[0]):
                temp_dist+=np.array(Binomials([[temp_n_a, p_a[i]], [temp_n_b, p_b[i]]]))/p_a.shape[0]
               
            temp_dist=self.convert_dist2percentage(temp_dist)
            if np.abs(mean_dist(np.arange(101), temp_dist)-mean_dist(np.arange(0, 100, 5), self.distribution))<20:
                temp_dist_5=[np.sum(temp_dist[:6])]+[np.sum(temp_dist[5*(i+1)+1:5*(i+2)+1]) for i in range(19)]
                temp_KLd=KLdivergence(temp_dist_5, self.distribution)
                if detail==True:
                    plt.plot(np.arange(len(temp_dist_5))*5+5, temp_dist_5, label='temp_dist')
                    plt.plot(np.arange(len(temp_dist_5))*5+5, self.distribution, label='self.distribution')
                    plt.title('n_a='+str(temp_n_a)+', n_b='+str(temp_n_b)+', KLd='+str(temp_KLd))
                    plt.show()

                if self.KLd>temp_KLd or np.isnan(self.KLd):
                    self.KLd=temp_KLd
                    self.n_a=temp_n_a
                   
            else:
                if detail==True:
                    print('Skip: n_a='+str(temp_n_a)+', n_b='+str(temp_n_b))
               
            temp_n_a=round(temp_n_a+1)
               
        self.n_b=self.n_total-self.n_a
        return
   
    def convert_dist2percentage(self, dist):
        percentage=[0 for i in range(101)]
        for i in range(len(dist)):
            percentage[round(100*i/(len(dist)-1))]+=dist[i]
           
        return percentage
   
def KLdivergence(p, q):
    p_range=[p[i]>1e-12 for i in range(len(p))]
    q_range=[q[i]>1e-12 for i in range(len(q))]
   
    double_true=np.sum([1 if p_range[i]==True and q_range[i]==True else 0 for i in range(len(p))])
    if double_true>0:
        out=np.sum([p[i]*np.log(p[i]/q[i]) if p_range[i]==True and q_range[i]==True else 0 for i in range(len(p))])
    else:
        out=np.inf
       
    return out

def mean_dist(value, dist):
    return np.dot(value, dist)

#data preparation
rawdata=pd.read_csv('data/sample_distribution.csv', index_col=0)
thrs=[5*i for i in range(20)]
n_questions=[100, 90, 80, 60, 70]
n_tests=len(n_questions)
n_subjects=100
subjects=np.zeros((2, n_subjects))

def dist2data(dist, thresholds):
    output=np.array([])
    i=0
    while i<len(thresholds)-1:
        if dist[i]!=0:
            interval=(thresholds[i+1]-thresholds[i])/dist[i]
            output=np.r_[output, np.arange(thresholds[i]+interval, thresholds[i+1]+interval, interval)]
           
        i+=1
    return output

subjects[0, :]=dist2data(np.array(rawdata.loc['exam0']), thrs)
subjects[0, :]/=100

for k in [0.4, 0.5, 0.6, 0.7, 0.8]:
    print('Calculation: k='+str(k))
    tests=[test(n_questions[i], np.array(rawdata.loc['exam'+str(i)])) for i in range(len(n_questions))]
    tests[0].n_a=n_questions[0]
    tests[0].n_b=0
   
    for i in range(n_tests-1):
        subjects[1]=k*subjects[0]
        tests[i+1].calc_n(subjects[0], subjects[1])
        print('tests['+str(i+1)+']: n_a='+str(tests[i+1].n_a)+', n_b='+str(tests[i+1].n_b)+', KLd='+str(tests[i+1].KLd))
       
    print()

------------------------------------------------------------

Calculation: k=0.4
tests[1]: n_a=19, n_b=31, KLd=0.07051862322854846
tests[2]: n_a=18, n_b=12, KLd=0.09610677830528232
tests[3]: n_a=19, n_b=11, KLd=0.16822053964137743
tests[4]: n_a=19, n_b=21, KLd=0.15326233787382407

Calculation: k=0.5
tests[1]: n_a=13, n_b=37, KLd=0.06463637859617347
tests[2]: n_a=15, n_b=15, KLd=0.09679698263389093
tests[3]: n_a=16, n_b=14, KLd=0.16493123276583527
tests[4]: n_a=15, n_b=25, KLd=0.1427783581892589

Calculation: k=0.6
tests[1]: n_a=4, n_b=46, KLd=0.06086984932650339
tests[2]: n_a=12, n_b=18, KLd=0.0953198995596042
tests[3]: n_a=13, n_b=17, KLd=0.16350196731623778
tests[4]: n_a=8, n_b=32, KLd=0.13279440803859335

Calculation: k=0.7
tests[1]: n_a=0, n_b=50, KLd=0.2184121139905782
tests[2]: n_a=6, n_b=24, KLd=0.09600257982634174
tests[3]: n_a=7, n_b=23, KLd=0.16232317488955234
tests[4]: n_a=0, n_b=40, KLd=0.13275637941572044

Calculation: k=0.8
tests[1]: n_a=0, n_b=50, KLd=0.8637338663088325
tests[2]: n_a=0, n_b=30, KLd=0.13691841637960955
tests[3]: n_a=0, n_b=30, KLd=0.18728352476626814
tests[4]: n_a=0, n_b=40, KLd=0.45602704138127004

------------------------------------------------------------

　KLd的にはk=0.6 or 0.7が良さそう。またkが0.7以上の時、tests[1]とtests[4]の問題がすべてB問題となり不自然なので、 k=0.6 を採用する。
　続き。
------------------------------------------------------------

k=0.6

tests=[test(n_questions[i], np.array(rawdata.loc['exam'+str(i)])) for i in range(len(n_questions))]
tests[0].n_a=n_questions[0]
tests[0].n_b=0

for i in range(n_tests-1):
    subjects[1]=k*subjects[0]
    tests[i+1].calc_n(subjects[0], subjects[1])
    print('tests['+str(i+1)+']: n_a='+str(tests[i+1].n_a)+', n_b='+str(tests[i+1].n_b)+', KLd='+str(tests[i+1].KLd))

result=np.zeros(np.sum(n_questions)+1)

for i in range(n_subjects):
    As=[[tests[j].n_a, subjects[0, i]] for j in range(n_tests)]
    Bs=[[tests[j].n_b, subjects[1, i]] for j in range(n_tests)]
    params=As+Bs
   
    subj_dist=Binomials(params)
    result+=np.array(subj_dist)

result=n_subjects*result/np.sum(result)

plt.plot(np.arange(result.shape[0]), result)
plt.xlabel('score')
plt.ylabel('subject')
plt.title('score distribution')
plt.show()

------------------------------------------------------------

------------------------------------------------------------

　この分布の線下面積は n_subjectsになるようにしている。つまり自分の点数が X点だった場合に、順位は score>X の部分の線下面積　S(score>X)を求めればよい。今回の場合 X=150だと面積 S(score>150) =63.4となるため、推定順位は64位となる。

　と、ここまで書いたところで、使用するサンプルデータをテキトーに作るのではなく、正解となる全受験者の正答率p_a, p_bと、各ブロックのA問題とB問題の配分と、合計点の分布を設定した上で作っていれば、上記の処理で得られる結果とどの程度一致するかを確かめられたなァと気づいた。しかしこの先暇つぶしする時間はないので、これ以上いじるつもりはない。

二項分布の和の分布

　notebookはこちら。

　確率変数の和の分布は単純な足し算ではできない、というのを実感するために、二項分布の和の分布を求めた。要は暇つぶしです。
　まず2つの場合。内積を使って各事象の確率を求めて、斜めに和をとった。
　正しくできているかどうかは、数学的に求められる確率変数の和の平均・分散が、入手した確率変数の和の分布から計算する平均・分散と一致するかどうかで確かめた。

--------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
import math

def Combination(n, r):
    return math.factorial(n)//math.factorial(n-r)//math.factorial(r)

def Binomial(param):
    N=param[0]
    p=param[1]
    return [Combination(N, k)*pow(p, k)*pow(1-p, N-k) for k in range(N+1)]

def Sum_Binomial(na, pa, nb, pb):
    A=Binomial([na, pa])
    B=Binomial([nb, pb])
    l=np.max([na, nb])
   
    #probability map
    prob_map=np.zeros((l+1, l+1))
    prob_map[:na+1, :nb+1]=np.dot(np.array(A).reshape(-1, 1), np.array(B).reshape(1, -1))
   
    #upper-left part of matrix
    hist_former=[np.sum([prob_map[i, s-i] for i in range(s+1)]) for s in range(l+1)]
   
    #lower-right part of matrix
    hist_latter=[np.sum([prob_map[-1-i, l-(minus-i)] for i in range(minus+1)]) for minus in range(l)]
   
    hist=hist_former+hist_latter[::-1]
    return hist

#set parameters
na, pa=10, 0.8
nb, pb=8, 0.7

#get Binomial distribution
A=Binomial([na, pa])
B=Binomial([nb, pb])

#get E(A+B) and V(A+B) by mathmatical method
print('mathmatical method: E(A+B)=E(A)+E(B)='+str(na*pa+nb*pb))
print('mathmatical method: V(A+B)=V(A)+V(B)='+str(na*pa*(1-pa)+nb*pb*(1-pb)))

#get A+B distribution
add_dist=Sum_Binomial(na, pa, nb, pb)

#get E(A+B)
mu=np.dot(np.arange(len(add_dist)), add_dist)
print('computational method: E(A+B)='+str(mu))

#get V(A+B)
rho=np.dot(np.arange(len(add_dist))**2, add_dist)-mu**2
print('computational method: V(A+B)='+str(rho))

#plot
plt.plot(np.arange(len(add_dist)), add_dist)
plt.show()
--------------------------------------------------------------
mathmatical method: E(A+B)=E(A)+E(B)=13.6
mathmatical method: V(A+B)=V(A)+V(B)=3.28
computational method: E(A+B)=13.600000000000001
computational method: V(A+B)=3.2799999999999727

--------------------------------------------------------------


　大丈夫そう。
　調子に乗って、確率変数が3つ以上になった場合もやった。

--------------------------------------------------------------
def Binomials(list_n_p):
    #list_n_p: 2-dimensioin vector (shape=(n, 2))
    #list_n_p[i]: i-th setting for Binomial distribution.
    #list_n_p[i][0]: n of i-th Binomial dist.
    #list_n_p[i][1]: p of i-th Binomial dist.

    hist=Binomial(list_n_p[0])
    for i in range(len(list_n_p)-1):
        pre_prob=hist
        add_prob=Binomial(list_n_p[i+1])
l=np.max([len(pre_prob), len(add_prob)])
        mx=len(pre_prob)+len(add_prob)-2
        #probability map
        prob_map=np.zeros((l, l))
        prob_map[:len(pre_prob), :len(add_prob)]=np.dot(np.array(pre_prob).reshape(-1, 1), np.array(add_prob).reshape(1, -1))
        #upper-left part of matrix
        hist_former=[np.sum([prob_map[i, s-i] for i in range(s+1)]) for s in range(l)]
        #lower-right part of matrix
        hist_latter=[np.sum([prob_map[-1-i, l-(minus-i)-1] for i in range(minus+1)]) for minus in range(l-1)]
        hist=hist_former+hist_latter[::-1]
        hist=hist[:mx+1]
       

    return hist

#set parameters
na, pa=10, 0.8
nb, pb=8, 0.7
nc, pc=6, 0.5

paramlist=[
    [na, pa],
    [nb, pb],
    [nc, pc]
]

#get Binomial distribution
A=Binomial([na, pa])
B=Binomial([nb, pb])
C=Binomial([nc, pc])

#get E(A+B+C) and V(A+B+C) by mathmatical method
print('mathmatical method: E(A+B+C)=E(A)+E(B)+E(C)='+str(na*pa+nb*pb+nc*pc))
print('mathmatical method: V(A+B+C)=V(A)+V(B)+V(C)='+str(na*pa*(1-pa)+nb*pb*(1-pb)+nc*pc*(1-pc)))

#get A+B+C distribution
add_dist=Binomials(paramlist)

#get E(A+B+C)
mu=np.dot(np.arange(len(add_dist)), add_dist)
print('computational method: E(A+B+C)='+str(mu))

#get V(A+B+C)
rho=np.dot(np.arange(len(add_dist))**2, add_dist)-mu**2
print('computational method: V(A+B+C)='+str(rho))

#plot
plt.plot(np.arange(len(add_dist)), add_dist)
plt.show()

--------------------------------------------------------------
mathmatical method: E(A+B+C)=E(A)+E(B)+E(C)=16.6
mathmatical method: V(A+B+C)=V(A)+V(B)+V(C)=4.779999999999999
computational method: E(A+B+C)=16.6
computational method: V(A+B+C)=4.779999999999916

--------------------------------------------------------------

　大丈夫そう。
　正直1時間もあればできると思ってたらだいぶ時間かかってしまった。一番引っかかったのが、はじめCombination()関数をC()と名付けていて、名前が分布Cと競合したところなのが完全に天然なんだが、そうはいってもなぁ。
　とはいえつい先日までリスト内包表記を全く知らなかったことを考えれば、だいぶ成長したか。

Fitzhugh-Nagumoモデルのパラメータ変化

　以下を実行したnotebookはこちら。

　Fitzhugh-Nagumoモデルのパラメータa, b, c, τを変化させた時の挙動の変化を調べた。

　まずaを変化させた。

　aが大きくなるにつれて波の振幅が狭くなる。主に上行脚よりも下降脚が急になっている印象を受ける。またわずかながらy軸方向にも平行移動しているように見えるが、振幅に対する効果に比べると無視していいレベルだと感じる。

　次にbを変化させた。

　大きくなるにつれて立ち上がりの閾値が上がり、それによって波と波の間隔が広がっているようにみえる。その考えならば、b=1で発火活動がみられないのは閾値を超えなくなったからとすると納得がいく。

　またその変化幅がbが大きいほど大きいように見えるので、分岐点に近い値で再度調べた。

　分岐点の近くだと少しの変化でかなり効果が出る。

　次にcを変化させた。

　cを大きくするに従って(1)閾値が下がり、(2)波の幅が短くなった ようにみえる。

　ここでbとcがいずれも閾値に影響したことから、二つの合わせ技を考えてみる。つまり先の結果では、b=1では発火活動がみられなかった。そこで、b=1のままcを大きくした場合に、発火活動がみられるようになるかを調べた。

　いけた。

　次にτを変化させた。

　τの変化によって時間軸方向に伸び縮みしているように見える。ただしある値より小さいと発火しなくなるようで、τ=0.4では発火活動がみられない。
　τの変化によって閾値が変化したのかはこの結果からはっきりとはわからないが、先ほどのように他の値との合わせ技を考えてみる。つまりτ=0.4としたまま、b, cを変化させて発火するようにすることはできるかを試す。
　まずはτ=0.4のままbを変化させる。

いけた。
　同様にτ=0.4のままcを変化させる。

これもいけた。
　これらの結果から一昨日の記事に対するアンサーは、「bを小さく」「cを大きく」すると閾値が下がり発火するようになり、「bを大きく」「cを小さく」すると閾値が上がって発火しにくくなりそうである。τに関しては、元の式みるとdw/dt全体にかかる係数だし、時間軸の伸び縮みくらいしか役割がない気がしているので一旦固定でいいかななどと思っている。

Fitzhugh-Nagumoモデルをつなげた

　以下を実行したnotebookはこちら。

　棘徐波複合の再現に挑戦してみようという機運になった。手始めに、以前1ノードで再現したFitzhugh-Nagumoモデルを連結するところから始めた。
　内部処理はオリジナルと同じものを用いた。外部入力は、他ノードからの入力と、「中枢」的なものからの入力の和とした。すなわち

　I=(他ノードからの入力)+(中枢からの入力)

とした。
　他ノードからの入力は、連結するノードのuの総和とした。(将来的には、これを発火時1、静止時0のようにしたい)
　中枢からの入力は、今後「視床由来の3Hz徐波」のようなものを考える可能性を考えて加えたものである。これはひとまず1+正弦波とした。1を足した理由は、マイナスになることのないようにしたかったため。
　Iの大きさが色々と変わるので、最終的な入力としては係数kとの積を取った。
　ノード数はひとまず5個。ネットワーク構成は、隣り合ったノードと結合して1周するような5角形。

　そのようにした上で散々色々なパラメータで実行した結果たどり着いたあるパラメータでの結果。

　最上部画像：全5ノードのuプロット
　下の5画像：各ノードのuプロット

　うまい具合にバラバラに発火するような感じになった。また上の結果はたまたまこうなったわけではなく、何度も実行した場合にもちゃんとバラバラに発火するような感じが再現された。なおkを大きくした場合には全ノードが同期してひたすら発火、小さくした場合には全ノードが静止する。
　本来は中枢からの入力を3Hz 波にしたかったので、"1+sin 6πt"とするべきなのだが、早とちりして"1+sin 3t"としてしまっていた。そこを正した場合で色々とパラメータ探索を試みたが、うまい具合にバラバラに発火するようなパターンを見つけられなかった。ある時の結果がこちら。

　見ればわかるように、3Hz波がモロに残ってしまっている。この結果をもってさっきの結果をもう一度みてみたら、同じように中枢由来の入力の周波数3/2π Hz程度で揺れていることに気付いた。
　Fitzhugh-Nagumoモデルで再現される波(以後、FN波)の周波数が2回のシミュレーションで変化しないこと、および先のシミュレーションではFN波と中枢由来の入力の周波数が比較的近い(または整数倍に近い)ことから、この中枢由来の入力の周波数がFN波の周波数と何らかの条件を満たした場合にのみ、上で見られたような適度に発火する系を再現できるような気がする。
　すなわち、このような「適度に発火する」「同期しない」の状況になるにはかなりパラメータの調整が必要そうである。そのため、(おそらく)同じ機構で働く神経細胞が様々なネットワーク・入力の下で意味のある発火を行うためには、各神経細胞がそれぞれ適度に発火するよう自分を調整するようにする仕組みが必要かなという気がする。それはパラメータ調節かもしれないし、入力がなくても発火するようなものかもしれないし、そこは神経生理学の教科書などを読まないとだめそうである。

今後の課題
・上のモデルに色んな変化いれて棘徐波ができないか試す
・全く入力が無い神経細胞は発火するか調べる
・FNモデルのパラメータをいじったときの波形の変化
・色んな自動調節機構を創作する 調べながら考える

何もしてないのに仮想環境が壊れた(機械弱者)

　諸々の解析ソフトがWindows未対応なため、解析時にはVirtualBoxでUbuntu16.04を立ち上げて解析を行っている。しかしその仮想環境がある日突如立ち上がらなくなるということが過去たびたびあり、そのたびに仮想環境ごと削除して一から環境構築をすることにより消極的解決を行ってきた。

　そして本日、再び仮想環境が立ち上がらなくなった。また一から環境構築しようかとも思ったのだが、環境構築が面倒な解析ソフトが以前より増えていたため、再度一から環境構築をするのはあまりにも面倒に感じた。そのため、仮想環境が完全に破壊されたらそれはそれでいいや(鼻ホジ)くらいの気概で、仮想環境の再生に取り組んだ次第。バックアップさまさまですな。

　起動時ターミナルがこちら

　"Welcome to emergency mode!" がパーティの司会風に脳内再生されるので「なにがウェルカムじゃこの」という気分になった。

　'journalctl -xb'でログを見たところ、上の画像の頭に表示されているエラー2点に加えて、共有フォルダがマウントできないエラーが吐かれていた。そこでVirtualBoxの設定から共有フォルダの共有設定を全て削除したり戻したりしたが、該当エラーは消えなかったので一旦保留。

　"Failed to start load kernel modules"でググったところ、多くの人が"apt-get update"と"dpkg --configure -a"あたりのコマンドをしろと言っていたので、ひとまず実行して再起動したが変化なし。そのまま色々ググって見回っていたら、このサイトにこんな記述がされていた。

'Remember it's quite common to have to repeat the above commands in no specific order'

繰り返せと(笑)

　他にやれることがあるわけでもないので、上記サイトに書かれた通りひたすら4コマンドと再起動を繰り返した。(5, 6周くらい？) 途中、文字化けにより真偽は不明だが、" 'sudo apt-get autoremove'しろ"と言われてそうな時には'sudo apt-get autoremove'を実行し、"/var/lib/dpkg/lockをopenできない"と言ってそうな時には'sudo rm /var/lib/dpkg/lock'を実行した。
　また途中わけわからんタイミングでGUI起動して完全に文字化けした何かが表示されたことがあった。エスケープ押したら<♦♦>と<♦♦♦>の2択になったので、無難にNOっぽい2文字の方選択したらターミナルに戻ったけど、あれなんだったんやろ(ヤケになりすぎてスクショ撮り損ねた)

　さすがに飽きてきたななどと思った矢先の起動画面がこちら。

1行目消えた！まじで繰り返して消えんのかよ！

　次のエラーについてググってたらまたさっきの4コマンドのうちの一つが解決策として書かれていたので実行するなどしたところ

増えた(笑)

　その後、ググリ着いたこのサイトの記載通りにコマンド実行し、再起動した起動画面がこちら。

消えた～！しかし起動しない。「頭の2行のエラーによってemergency modeになっている」と思いこんでいたが、どうやらそうではないらしいということが判明する。
　残すは共有フォルダの問題しかない。"Welcome to emergency mode"でググってこのサイトに行きつき、vimで/etc/fstabを開いてみた。

　この時、共有フォルダ 'share' はVirtualBoxマネージャーから共有を外していたので、「shareについての記述があるのはおかしい！(決めつけ)」と判断し、行ごと削除して保存。そして一旦仮想環境をシャットダウンし、VirtualBoxマネージャーから再度共有フォルダに'share'を設定して、そして起動！

　ログイン画面出た！中身もそっくりそのまま残っていました。めでたしめでたし。


　さて、振り返ってみると、実は最後の共有フォルダの設定だけが問題だったのではないかという気がしてくる。そういわれてみると、正常にログインできていた頃も起動時に数行エラーを吐かれていた気がしないでもない。もう復活してしまったので試せないのだが、同じような使い方していれば定期的に仮想環境が壊れる()ので、次の機会には共有フォルダの件だけで解決するかどうかを試そうと思う。またその場合には何故共有フォルダの設定が狂ったのかが不明だが、これは先日のWindows自動更新によるものだと信じたい。

　決して！俺は！何もしていない！(機械弱者)

Sample Entropyで遊ぶ

(Sample Entropyの日本語ソースがあまりに少ないので最初に断っておくと、ただの学生の遊びブログなのであまりあてにしないでいただければと。)
(むしろ詳しい人いたら教えてください。)

　時系列データのエントロピーを示すSample Entropyという指標に出会ったので遊んでみる。長さNの時系列データ U={u(1), u(2), ... , u(N)}と、定数m(通常m=2)とした時、Sample Entropyは次のようにして求める。

(1) x(i)={u(i), ... , u(i+m-1)}とし、X={x(1), x(2), ... , x(N-m+1)}を得る。
(2) Ci(r)を求める。Ci(r)=(d[x(i), x(j)<r] となるjの数(i≠j))/(N-m+1) (ただしd[a, b]はa, bのChebyshev distance)(r=0.2*データの標準偏差)
(3) Φm(r)を求める。Φm(r)=Σ[i=1, N-m]Ci(r)/(N-m)
(4) 同様にΦm+1(r)を求める。
(5) SampEn=ln[Φm(r)/Φm+1(r)]

　現時点で不明点は2つ。
　なぜmを1足したものとの関係を考えることでエントロピーと扱えると考えたのかがわからない。これについてはSample Entropyの元となったApproximate Entropy以前へ遡らなければならなさそうなので今後の課題ということで。
　2点目は、近年の文献では上に示したようなΦmとΦm+1の比の対数を取っているが、Richman(2000)では分母が揃うようにした二つの微妙に異なるAi, Biの比の対数を取っていた。Nがある程度大きければ「あんま変わんないんじゃね？」という話にでもなったのだろうか。当然誤読の可能性も十分にある。こちらも今後の課題。
　とりあえず以下ではΦmを用いた方法でSample Entropyを計算した。色んなデータを与えてどうSample Entropyが変化するかをみてみる。

(1)正弦波の振幅変化
　信号源を x=k*sin(T) とし、このkを変化させた時に、0.1秒間隔でサンプリングしたデータのSample Entropyを調べた。k=[1, 2, ... , 10]

　一定。rが標準偏差に合わせて変化するので、振幅の影響は基本的にはないということですな。

(2)正弦波の周波数変化
　信号源を x=sin(k*T)とし、このkを変化させた時に、0.1秒間隔でサンプリングしたデータのSample Entropyを調べた。k=[1, 2, ... , 10]

　ある値を境に大きく変化するというのは違和感があるので、もう少し多くのkで調べてみる。k=[1, 2, ... , 100]としてみる。

　信号源が正弦波であることにより「たまたまSample Entropyが大きくなるような信号値の並びになった」やその逆が起きることで乱高下していると思われる。こんなキッチリ周期的な信号源というのは現実的ではないのであまり気にしなくてもいいと思うが、一応周波数変化により乱高下することが「ありうる」ということは知っておいていいかもしれない。また周波数変化は「時間軸の伸び縮み」と同義なので、異なる時系列データのSample Entropyを比較するには、サンプリング間隔をそろえる必要はありそう。

(3) サンプリング数の変化
　信号源を x=sin(T)とし、0.1秒間隔でサンプリングするデータ数Nを変化させた時のSample Entropyを調べた。N=[200, 300, ... , 1000]

　Sample Entropyはその元となったApproximate Entropyに比べて「データ数の影響を受けにくい」と記述されていた気がするが、このデータでは漸減している。現実データではあまり影響しないのかもしれないが、時系列データを比較する際にはデータ数もそろえた方が無難そう。Nは出来れば1000以上がよいとしていた文献があったのだが、確かに1000を超えてこればだいぶ安定しそうではある。

(4) 正弦波を重ね合わせた信号源
　信号源を x=Σ[i=1, k]sin(i(T+δi))/k とし 0.1秒間隔でサンプリングしたデータのSample Entropyを調べた。重ね合わせる正弦波の始点をそろえると信号のバリエーションが少ないので、要素sin(iT)/kをそれぞれ[0, 2π/i]からランダムに選択したδiだけ平行移動させた。乱数を用いるため、各kで20回ずつ実行して求めたSample Entropyを平均したものを用いた。k=[1, 2, ... , 30]

かなりきれいに上昇してくれてとてもうれしい気持ちになった。
(5) 階段関数
　正弦波ばかりで飽きたので、階段関数を解析してみる。信号源x=T%k (T=[1, 2, ... , N])、サンプリング間隔は1、N=200、k=[2, 3, ... , 30]

　ファッ！？と思ったが、rが1を超えたことで一気に増加しただけでした。これも現実データではあまりないと思われるが、解析に際してrを多少上下させてみる必要はあるかもしれない。


　このSample Entropyを身体動作や心電図、脳波で計算して、「加齢でエントロピーが～」みたいな研究がいくつか見られるが、上で述べたように「mとm+1の場合を比較することで何故エントロピーといえるのかがわからない」ので、現時点ではあまりパッとしないなぁという印象。ここがすっきり解決してこれがほんまにエントロピーだと言えそうということになったら、解釈についてもすっきりしそうではある。

上記を実行したjupyter notebookはコチラ

ボーダーライン再考

　昨年末の総合試験前、「点数を伸ばすことによって"マイナス k SD"のボーダーラインを下げる人の条件」を、平均M、標準偏差SDとして「点数でM+SD/k以上 (偏差値で50+10/k以上) の人」と考えたら、次のような指摘が来ていた。

その後ウーンと考えていたまま忘れ去っていたのをふと思い出し、どこが間違っていたのか結局わかっていないので、導出過程を貼っておいてまたコメントでももらおうという記事です。前回のコメ主でも通りすがりの人でも誰でもいいのでご指摘ください。

相関行列からネットワークへ

　(昨日の記事から色々追加・修正したので、新しく書き直しました。)
　FisherのZ変換を施された相関行列でネットワーク解析を行う。

お品書き
1. FisherのZ変換の逆作業をして生の相関行列を作る
2. 有意な相関係数のみ残す
3. 相関行列からsemi-metricな距離行列を作成する
4. semi-metricな距離行列から重み付き隣接行列を作成する

1. FisherのZ変換の逆作業をして生の相関行列を作る
　某解析ソフトは相関行列を計算した結果にFisherのZ変換を施すために値が[-1, 1]から溢れる。邪魔くさいので元の相関係数に戻したい。
　入力のcsvファイルを逆関数にぶちこむだけのシンプルなもの。

matrix_Z2R.py

--------------------------------------------------

#Usage

#python3 matrix_Z2R.py [input filename] [output filename]

#

#Input

#[input filename]: csv file (correlation coefficient matrix converted Z-score by Fisher's method)

#[output filename]: output filename

#

#Output

#[output filename]: raw correlation coefficient matrix

import sys

import numpy as np

filename=sys.argv[1]

output=sys.argv[2]

matrix=np.loadtxt(filename, delimiter=',')

matrix=(np.exp(2*matrix)-1)/(np.exp(2*matrix)+1)

np.savetxt(output, matrix, delimiter=',')

--------------------------------------------------

2. 有意な相関係数のみ残す

　プログラムどうこうというより、無相関検定の復習になった。

Rmat_significnat.py

--------------------------------------------------

#Usage

#python3 Rmat_significant.py [input filename] [output filename] [number of sample] [significance level]

#

#Input

#[input filename]: Input csv filename (correlation coefficient matrix)

#[output filename]: Output csv filename

#[number of sample]: Number of samples to calculate each correlation coefficient

#[significance level]: significance level (usually 0.05)

#

#Output

#[output filename]: adjacency matrix filled with 0 or 1

import sys

import numpy as np

import scipy.stats

filename=sys.argv[1]

output_filename=sys.argv[2]

n_sample=np.int(sys.argv[3])

alpha=np.float(sys.argv[4])

#load inputfile

matrix=np.loadtxt(filename, delimiter=',')

matrix[np.isnan(matrix)]=0

#get T-value matrix by correlation coefficient matrix

Tmat=np.abs(matrix)*pow(n_sample-2, 1/2)/((1-matrix**2)**(1/2))

#get T threshold for significance

T_sig=-1*scipy.stats.t.ppf(q=[alpha/2], df=n_sample-2)

#set unsignificant values 0

matrix[Tmat<T_sig]=0

#set minus R 0

matrix[matrix<0]=0

#Output

np.savetxt(output_filename, matrix, delimiter=',')

--------------------------------------------------

3. 相関行列からsemi-metricな距離行列を作成する

　話が前後するのだが、4で距離を重みにするときに単に「距離の逆数を重みとする」のではなく、「距離に1を足した数の逆数を重みとする」のが個人的にすごいグッときたんすよね。以前遊んでいた時に、ゼロの逆数が無限大になるのがとても扱いにくかったので、なるほどなぁと。

　そんなわけで距離行列作成では「重み(相関係数)の逆数から1を引く」作業をしている。

tnorm_r2len.py

--------------------------------------------------

#Usage

#python3 tnorm_r2len.py [input filename] [output filename]

#Input

#

#[input filename]: input weighted matrix csv filename

#[output filename]: output filename

#

#Output

#

#[output filename]: t-norm length matrix

import sys

import numpy as np

inf=np.inf

input_filename=sys.argv[1]

output_filename=sys.argv[2]

matrix=np.loadtxt(input_filename, delimiter=',')

matrix[matrix==0]=inf

matrix=matrix**(-1)

matrix-=1

matrix[matrix<0]=inf

np.savetxt(output_filename, matrix, delimiter=',')

--------------------------------------------------

4. semi-metricな距離行列から重み付き隣接行列を作成する

　現時点では有意とならない相関となったノード間には結合が存在しないが、結合が存在しないノード間にも重みを見出したいとする。そこで、距離行列からダイクストラ法でノード間の最短距離を求め、Dombi t-normで変換したものを重みとしているものがあって、なるほどなぁと思いそれを採用。なお対角成分はゼロにセットしなおすので、2分グラフは未対応。

tnorm_len2wei.py

--------------------------------------------------

#Usage

#python3 tnorm_len2wei.py [input filename] [output filename]

#

#Input

#[input filename]: length matrix csv file

#[output filename]: save weighted matrix as the argv

import sys

import numpy as np

import networktools

input_filename=sys.argv[1]

output_filename=sys.argv[2]

matrix=np.loadtxt(input_filename, delimiter=',')

matrix=networktools.shortest_path_length(matrix)

matrix=(matrix+1)**(-1)

for i in range(matrix.shape[0]):

 matrix[i, i]=0

np.savetxt(output_filename, matrix, delimiter=',')

--------------------------------------------------

　そんなわけで上4つを順に呼べば、某ソフトで出力される相関行列からネットワークの距離行列と重み付き隣接行列を作成することができるようになった。

PS
　久々に文字サイズがそろわない。大量のHTMLが突っ込まれているのはわかるんだが、どの作業の時に突っ込まれるのか全くわからん…