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

学生さんに宿題で出したので自分でも書いてみた。
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()