基本の作図

例えばこんな感じのデータフレームがあったとする．

    x         y
0   2  0.954025
1   3  0.146810
2   1  0.409961
3   1  0.164558
4   3  0.782152
5   2  0.905869
6   3  0.210528
7   1  0.437970
8   1  0.801206
9   3  0.089576
10  2  0.960357
11  2  0.670732

このデータフレームから標準誤差付きの図を描こうとすると，こんな感じ．

import numpy as np
import matplotlib.pyplot as plt
 
m = df.pivot_table(index='x', values='y', aggfunc='mean')
e = df.pivot_table(index='x', values='y', aggfunc='sem')
m.plot(xlim=[0.8, 3.2], yerr=e)

このように，yerrに誤差の大きさを指定することで，誤差棒を作成することができる．

信頼区間を用いた作図

そこで，信頼区間の長さを求めるcilenを定義して利用することにする．

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt


def cilen(arr, alpha=0.95):
    if len(arr) <= 1:
        return 0
    m, e, df = np.mean(arr), stats.sem(arr), len(arr) - 1
    interval = stats.t.interval(alpha, df, loc=m, scale=e)
    cilen = np.max(interval) - np.mean(interval)
    return cilen

m = df.pivot_table(index='x', values='y', aggfunc='mean')
e = df.pivot_table(index='x', values='y', aggfunc=cilen)
m.plot(xlim=[0.8, 3.2], yerr=e)

信頼区間付きの図を作成することができた．

おまけ

信頼区間を計算する方法は「nかn-1か」問題のせいでごちゃごちゃしてます．

「車輪の再発明は避けよ」に従って，とりあえず「信頼区間 Python」とかで検索した人は，

あたりを見て，微妙に違うことに困惑しつつも両方試してみて，そして結果が違うことに辟易することになると思います (なりました)．

Rと同じ答えになるのは1の方．違うのは2の方です．ちょっと試してみます．

> x <- c(1, 1, 3, 3)
> t.test(x)

    One Sample t-test

data:  x
t = 3.4641, df = 3, p-value = 0.04052
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 0.1626138 3.8373862
sample estimates:
mean of x
        2

1の方法

import numpy as np
from scipy import stats

alpha     = 0.95
data = [1, 1, 3, 3]

mean_val = np.mean(data)
sem_val  = stats.sem(data)  # standared error of the mean
ci       = stats.t.interval(alpha, len(data)-1, loc=mean_val, scale=sem_val)

print(ci)

>> (0.16261376896260105, 3.837386231037399)

2の方法

import math
import numpy as np
from scipy import stats

alpha     = 0.95
data = [1, 1, 3, 3]

mean_val = np.mean(data)
std_val = np.std(data)
ci = stats.t.interval(alpha,len(data)-1,loc=mean_val,scale=std_val/math.sqrt(len(data)))
print(ci)

>> (0.40877684735786857, 3.5912231526421312)

今回は天下のRに従うことにしたので1を採用しています．

なお，2は何が違うのかというと，ciを計算する箇所の最後の方

math.sqrt(len(data))

これです．nで割ってます．しかし推測統計を行うのであればn-1で割った方が良いです．t分布を仮定するからですね．実際，

math.sqrt(len(data) - 1)

とすると，方法2の答えもRと完全に一致します．