「データ解析のための統計モデリング入門」の2章のグラフをPythonで描画する

chap2.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from collections import Counter
from scipy.stats import poisson
import rpy2.robjects as robjects

if __name__=='__main__':

    rdata = r('load("data.RData")')
    data = np.array(cdata)
    cdata = r('data')
    data = np.array(cdata)

    ## Min.
    print np.min(data)
    ## 1st Qu.
    print np.percentile(data, 25)
    ## Median
    print np.median(data)
    ## Mean
    print np.mean(data)
    ## 3rd Qu.
    print np.percentile(data, 75)
    ## Max.
    print np.max(data)

    ## table?
    print Counter(data)

    ## Fig.2.2
    plt.hist(data,bins=arange(-0.5,8.5,1))
    plt.xlabel('data')
    plt.ylabel('Frequency')
    plt.title('Histogram of data')

    ## Sample Variability
    print np.var(data)

    ## Sample Standard Deviation
    print np.std(data)
    print np.sqrt(np.var(data))

    x = arange(0,10)
    prob = poisson.pmf(x,3.56)

    ## Fig.2.4
    plt.plot(x,prob,'bo--')
    plt.xlabel('y')
    plt.ylabel('prob')

    ## Fig.2.5
    plt.hist(data,bins=arange(-0.5,8.5,1),color='white')
    plt.plot(x,prob*50,'bo--')
    plt.xlabel('data')
    plt.ylabel('Frequency')
    plt.title('Histogram of data')

    ## Fig.2.6
    frame= pd.DataFrame()
    x = arange(0,20)
    mlist = [3.5,7.7,15.1]

    for m in mlist:
        prob = poisson.pmf(x,m)
        frame[str(m)] = prob

    frame.plot(kind='line',style=['bo-','bd-','b^-'])
    plt.xlabel('y')
    plt.ylabel('prob')
    plt.legend(title='lambda')

    ## Fig.2.7
    mlist = [2.0, 2.4, 2.8, 3.2, 3.6, 4.0, 4.4, 4.8, 5.2]
    plt.figure(figsize=(12,12))
    for i in arange(1,10):
        m = mlist[i-1]
        subplot(330+i)
        x = arange(0,10)
        logl = sum([d*log(m) - m - sum(log(arange(1,d+1))) for d in data])
        prob = poisson.pmf(x,m)
        plt.hist(data,bins=arange(-0.5,8.5,1),color='white')
        plot(x,50*prob,'bo--')
        plt.xlim(-0.5,8)
        plt.ylim(0,15)
        plt.text(5,13,'$\lambda='+str(m)+'$\n $\log L='+str(round(logl,1))+'$',fontsize=10)