信頼区間の計算で使用するために t 統計を取得するための Python 関数 (または、存在しない場合は独自の関数を作成する) を探しています。
このようなさまざまな確率/自由度の答えを示す表を見つけましたが、任意の確率に対してこれを計算できるようにしたいと考えています。この自由度にまだ慣れていない人にとっては、サンプルのデータ ポイントの数 (n) は -1 で、上部の列見出しの数字は確率 (p) です。たとえば、両側の有意水準 0.05 が使用される場合n回のテストを繰り返した場合、結果が平均+/-信頼区間内に収まるという95％の信頼度の計算に使用するtスコアを調べています。
scipy.stats 内でさまざまな関数を使用することを調べましたが、上記の単純な入力を許可しているように見えるものはありません。
=TINV(0.05,999)Excel には、たとえば 1000 のサンプルの tスコアを取得するための簡単な実装があります。ここでは、95% の信頼が必要です。
これまでのところ、信頼区間を実装するために使用したコードを次に示します。現在、t スコアを取得するのに非常に大雑把な方法を使用していることがわかります (perc_conf にいくつかの値を許可し、それが正確ではないことを警告するだけです)。サンプル < 1000):
# -*- coding: utf-8 -*-
from __future__ import division
import math
def mean(lst):
# μ = 1/N Σ(xi)
return sum(lst) / float(len(lst))
def variance(lst):
"""
Uses standard variance formula (sum of each (data point - mean) squared)
all divided by number of data points
"""
# σ² = 1/N Σ((xi-μ)²)
mu = mean(lst)
return 1.0/len(lst) * sum([(i-mu)**2 for i in lst])
def conf_int(lst, perc_conf=95):
"""
Confidence interval - given a list of values compute the square root of
the variance of the list (v) divided by the number of entries (n)
multiplied by a constant factor of (c). This means that I can
be confident of a result +/- this amount from the mean.
The constant factor can be looked up from a table, for 95% confidence
on a reasonable size sample (>=500) 1.96 is used.
"""
if perc_conf == 95:
c = 1.96
elif perc_conf == 90:
c = 1.64
elif perc_conf == 99:
c = 2.58
else:
c = 1.96
print 'Only 90, 95 or 99 % are allowed for, using default 95%'
n, v = len(lst), variance(lst)
if n < 1000:
print 'WARNING: constant factor may not be accurate for n < ~1000'
return math.sqrt(v/n) * c
上記のコードの呼び出し例を次に示します。
# Example: 1000 coin tosses on a fair coin. What is the range that I can be 95%
# confident the result will f all within.
# list of 1000 perfectly distributed...
perc_conf_req = 95
n, p = 1000, 0.5 # sample_size, probability of heads for each coin
l = [0 for i in range(int(n*(1-p)))] + [1 for j in range(int(n*p))]
exp_heads = mean(l) * len(l)
c_int = conf_int(l, perc_conf_req)
print 'I can be '+str(perc_conf_req)+'% confident that the result of '+str(n)+ \
' coin flips will be within +/- '+str(round(c_int*100,2))+'% of '+\
str(int(exp_heads))
x = round(n*c_int,0)
print 'i.e. between '+str(int(exp_heads-x))+' and '+str(int(exp_heads+x))+\
' heads (assuming a probability of '+str(p)+' for each flip).'
この出力は次のとおりです。
1000 回のコイン投げの結果が 500 の +/- 3.1% 以内、つまり 469 ～ 531 回の表になることを 95% 確信できます (各コイン投げの確率を 0.5 と仮定)。
また、範囲のt 分布を計算し、必要な確率に最も近い t スコアを返すことも検討しましたが、式の実装に問題がありました。これが関連しており、コードを見たいかどうか教えてください。
前もって感謝します。