Demo entry 6686497

TastingTea

   

Submitted by anonymous on Dec 23, 2017 at 15:22
Language: Python. Code size: 4.1 kB.

def massivetest(n, rate, testnum, threshold):
    import numpy as np
    from scipy.special import comb

    def judge(n, rate):                   #代表被试做判断的函数。默认前n杯茶为茶,后n杯为奶
        import random
        tearandom = 0
        for i in range(0, n):             #先根据rate生成一个数,代表多少杯茶需要乱蒙
            if random.random() < rate:
                tearandom += 1

        index = [1] * (2 * n)
        tea = [1] * n
        for i in range(0, 2 * n):
            index[i] = i

        for i in range(0, n - tearandom):                   #先选出确定的茶
            tea[i] = i
            index[i] = index[2 * n - 1 - i]
        for i in range(0, tearandom):                       #再选出乱蒙的茶
            r = int((tearandom + n - i) * random.random())
            tea[n - tearandom + i] = index[r]
            index[r] = index[tearandom + n - i - 1]
        #for i in range(0, n):
        #    print(tea[i])
        return tea                     #返回判断是茶的样本的编号数组(0-7)

    def runtest(n, rate):              #做实验的函数,返回选对的茶数
        tea = [1] * n
        tea = judge(n, rate)
        correctcount = 0
        for i in range(0, n):
            if tea[i] < n:
                correctcount += 1
        #print(correctcount)
        return correctcount

    def getprob(k, n):                 #返回2n杯茶时,在零假设条件下选对k杯的概率
        a = comb(n, k)
        b = comb(n, n - k)
        c = comb(2 * n, n)
        prob = a * b / c
        return prob

    for i in range(0, n):
        probsum = 0
        for j in range(0, i + 1):
            probsum += getprob(j, n)
        if probsum > threshold:
            cutoff = n + 1 - i          #cutoff值即是我们根据拒绝零假设的阈值判断得到的,表示大于等于多少样本正确时可以拒绝零假设
            break

    result = [1] * testnum
    judgeability = [0] * testnum
    ratio_of_admit = 0

    for i in range(0, testnum):
        result[i] = runtest(n, rate)            #做实验并将判断正确数分别储存在数组result中,将最终判断储存在数组judgeability中
        if result[i] >= cutoff:
            judgeability[i] = 1
        ratio_of_admit += judgeability[i]
    ratio_of_admit /= testnum                   #计算做多次实验后承认被试判断能力情形的频率

    return ratio_of_admit                       #输出做多次实验后承认被试判断能力情形的频率

    
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif']=['SimHei']
plt.figure(figsize=(15,10))
plt.xlabel("实际瞎蒙率")
plt.ylabel("拥有判断能力的显著性水平")
plt.title("被试拥有判断能力的显著性水平-被试实际瞎蒙率的相关图")

rate = [1] * 101
admit = [1] * 101
for i in range(0, 101):
    rate[i] = i / 100
    admit[i] = massivetest(4, rate[i], 1000, 0.05)
plt.plot(rate, admit, label = 'n = 4', color = 'r')

rate = [1] * 101
admit = [1] * 101
for i in range(0, 101):
    rate[i] = i / 100
    admit[i] = massivetest(6, rate[i], 1000, 0.05)
plt.plot(rate, admit, label = 'n = 6', color = 'g')

rate = [1] * 101
admit = [1] * 101
for i in range(0, 101):
    rate[i] = i / 100
    admit[i] = massivetest(10, rate[i], 1000, 0.05)
plt.plot(rate, admit, label = 'n = 10', color = 'b')

rate = [1] * 101
admit = [1] * 101
for i in range(0, 101):
    rate[i] = i / 100
    admit[i] = massivetest(30, rate[i], 1000, 0.05)
plt.plot(rate, admit, label = 'n = 30', color = 'brown')

rate = [1] * 101
admit = [1] * 101
for i in range(0, 101):
    rate[i] = i / 100
    admit[i] = massivetest(100, rate[i], 1000, 0.05)
plt.plot(rate, admit, label = 'n = 100', color = 'orange')

rate = [1] * 101                                      #这部分是第三张图的代码
admit = [1] * 101
for i in range(0, 101):
    rate[i] = i / 100
    admit[i] = massivetest(2, rate[i], 1000, 0.2)
plt.plot(rate, admit, label = 'n = 2, threshold = 0.2', color = 'r')

rate = [1] * 101
admit = [1] * 101
for i in range(0, 101):
    rate[i] = i / 100
    admit[i] = massivetest(3, rate[i], 1000, 0.055)
plt.plot(rate, admit, label = 'n = 3, threshold = 0.055', color = 'b')

rate = [1] * 101
admit = [1] * 101
for i in range(0, 101):
    rate[i] = i / 100
    admit[i] = massivetest(4, rate[i], 1000, 0.05)
plt.plot(rate, admit, label = 'n = 4, threshold = 0.05', color = 'g')

plt.legend()                             #这个是共有的代码
plt.show()

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).