2018年11月12日月曜日

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

　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問題の配分と、合計点の分布を設定した上で作っていれば、上記の処理で得られる結果とどの程度一致するかを確かめられたなァと気づいた。しかしこの先暇つぶしする時間はないので、これ以上いじるつもりはない。
投稿者 時刻:

0 件のコメント:

コメントを投稿

登録: コメントの投稿 (Atom)