ヘルスライスでバーサーク

「『カザンの闘技場』でもらえるヘルスライス(42D)でバーサークするとどの位ヒットが出るのか?」という子供の頃の疑問に答えるスクリプトを誘惑に負けて作成。関数を再帰的に使うことの練習ということで(汗)。

100万回トライして
平均 :  357
標準偏差 :  23
最大 :  501
最小 :  268

非バーサーク時(期待値 147)の2.4倍。

分布でなんか変なヒゲが生えているのは何なのだろう...?

ハイパーバーサークだと1万回トライして
平均 :  499
標準偏差 :  51
最大 :  765
最小 :  348

...仕事します。

import random
import numpy
import matplotlib.pyplot as plt

def normal_dice(dice_num) :
    vals = []
    total_val = 0
    
    random.seed()
    for num in range(0,dice_num):
        vals.append(random.randint(1, 6))
    
    #print vals
    
    total_val = sum(vals)
    
    return total_val


def berserk_dice(dice_num) :
    vals = []
    total_val = 0
    
    random.seed()
    for num in range(0,dice_num):
        vals.append(random.randint(1, 6))
    
    #print vals
    
    total_val = total_val + sum(vals)
   
    for num in range(1,6):
        num_of_same = vals.count(num)
        if num_of_same > 1 :
            #print "Bersrk by %d dices of %d marks" % (num_of_same,num) 
            total_val=total_val+berserk_dice(num_of_same)
    
    return total_val

def hyper_berserk_dice(dice_num) :
    vals = []
    total_val = 0
    
    random.seed()
    for num in range(0,dice_num):
        vals.append(random.randint(1, 6))
    
    #print vals
    
    total_val = total_val + sum(vals)
   
    for num in range(1,6):
        num_of_same = vals.count(num)
        if num_of_same > 1 :
           # print "Bersrk by %d dices of %d marks" % (num_of_same,num) 
            total_val=total_val+hyper_berserk_dice(num_of_same+1)
    
    return total_val
    
def main():
    dice_num = input("Input dice_num :")
    trial = input("Input Num. of Trials :")
    vals=[]

    for num in range(0, trial) :
        vals.append(berserk_dice(dice_num)) 
        if num%10000==0:
            print num

    vals=numpy.array(vals)
    print "Mean : ", numpy.mean(vals)
    print "SD : ", numpy.std(vals)
    print "Max : ", vals.max()
    print "Min : ", vals.min()
   
    plt.rc('font', **{'family': 'serif'})
    fig = plt.figure()
    
    ax = fig.add_subplot(111)
    ax.hist(vals, bins=100)
    ax.set_title('Berserk with %d dices! (%d trials)' % (dice_num, trial))
    ax.set_xlabel('Hits!')
    ax.set_ylabel('Frequency')

    plt.show()    
    return None 
    
main()

中心極限定理をテストするスクリプト

学生さんに宿題で出したので自分でも書いてみた。
Python (matplotlib, numpy) 便利。
Enthought Canopyも便利。

# -*- coding: utf-8 -*-
import random
import numpy
import matplotlib.pyplot as plt

'''
あまり格好は良くないけど、中心極限定理をテストするスクリプト。
歪んだベータ分布から正規分布が出てくる。
以下を参考にした。
http://akiyoko.hatenablog.jp/entry/2013/06/07/213448
http://docs.python.jp/2/library/random.html
'''

def main():
    random_val=[]
    
    print("This is a test of central limit theorem by beta functions.")
    alpha = input("Input alpha :")
     
    beta = input("Input beta :")
    
    N = input("Input N :")
    trial = input("Input Num. of Trials :")

    random.seed()

    for num in range(0, N):
        random_val.append(random.betavariate(alpha,beta))
        
    random_val=numpy.array(random_val)

        
    mean=numpy.mean(random_val)    
    print "mean=%.3f" % mean

    sigma=numpy.std(random_val)
    print "sigma=%.3f" % sigma
    
    vals=[]
    
    for num0 in range(0,trial):
        val=0.0
        for num1 in range(0, N):
            val=val+random.betavariate(alpha,beta)
        vals.append(val) 
                 
    vals=numpy.array(vals)

    plt.rc('font', **{'family': 'serif'})
    fig = plt.figure()
    
    ax0 = fig.add_subplot(311)
    ax0.hist(random_val, bins=100, range=(0, 1), normed=False, facecolor='r', alpha=0.8)
    ax0.set_xlim(0, 1)

    y,x=numpy.histogram(random_val, bins=100, range=(0, 1), normed=False)

    ax0.set_title('Beta distribution')
    ax0.set_xlabel('Random variable')
    ax0.set_ylabel('Frequency')
    ax0.text(0.8, y.max()*0.8,r'alpha=%.2f' % alpha)
    ax0.text(0.8, y.max()*0.7,r'beta=%.2f' % beta)
    ax0.text(0.8, y.max()*0.6,r'N=%d' % N)

    ax1 = fig.add_subplot(313)
    ax1.hist(vals, bins=100, normed=False, facecolor='g', alpha=0.8)
    ax1.set_title('Distribution of sums Trial=%d' % trial)
    ax1.set_xlabel('Sum of random variables')
    ax1.set_ylabel('Frequency')
  
    plt.show()
    
    return None

main()