二つの母平均の差の区間推定

ふたつの母平均の区間推定

  • 正規分布に従う確率変数の和も正規分布に従う

  • ふたつの母集団から取り出した標本の平均が,xとyなら

    • xとyの差の平均は x - yになる
    • xとyの差の分散は、
      • 母分散がわかっているなら、両方の標本分散の和、つまり、(xの標本分散 = xの母分散 / xの標本数) + (yの標本分散 = yの母分散 / yの標本数)   - 母分散がわかっていないなら、xの不偏分散とyの不偏分散を合成する。
        • つまり、((xの標本数 - 1) * xの不偏分散 + (yの標本数 - 1) * yの不偏分散) / (xの標本数 + yの標本数 - 2)
        • これをプールした分散という。
  • 母分散がわかっているなら, z値で母平均を区間推定

    • z値は、((x - y) - 差の母平均) / √ (xの標本分散 = xの母分散 / xの標本数) + (yの標本分散 = yの母分散 / yの標本数)
    • 95%信頼区間なら、
    • -1.96 < ((x - y) - 差の母平均) / √ (xの標本分散 = xの母分散 / xの標本数) + (yの標本分散 = yの母分散 / yの標本数) < 1.96
    • 差の母平均 = (x - y) +- 1.96 * √ (xの標本分散 = xの母分散 / xの標本数) + (yの標本分散 = yの母分散 / yの標本数)
  • 母分散がわかっていないなら, 自由度(xの標本数 + yの標本数 - 2)のt値で母平均を区間推定

    • t値は、((x - y) - 差の母平均) / √ ((1 / xの標本数) + (1 / yの標本数) * プールした分散)
    • 95%信頼区間、 xの標本数が34,yの標本数が33なら、
    • -2 < ((x - y) - 差の母平均) / √ ((1 / xの標本数) + (1 / yの標本数) * プールした分散) < 2
  • 差の母平均 = (x - y) +- 2 * √ ((1 / xの標本数) + (1 / yの標本数) * プールした分散)

  • 対応のある標本の場合は、ひとつの標本であると想定する。自由度は標本数 - 1

母集団比率の推定/検定

  • ある大都市で1200世帯を無作為に抽出
  • 65世帯に要介護の家族がいる
  • この年全体での要介護者のいる世帯の割合はどのくらいか
  • 95%区間推定する

  • 以下の値は標準正規分布する

    • (サンプルの中での要介護者のいる世帯の割合 - 母集団の中での要介護者のいる世帯の割合) / √((サンプルの中での要介護者のいる世帯の割合)(1 -サンプルの中での要介護者のいる世帯の割合 ) / 標本数 )
  • 従って、95%区間推定なら、以下の範囲を成立させるpが信頼区間
    • -1.96 <= ( 0.054 - p ) / √((0.054 * 0.946) / 1200) <= 1.96
import math

upper = 0.054 + 1.96 * math.sqrt((0.054 * 0.946) / 1200)
lower = 0.054 - 1.96 * math.sqrt((0.054 * 0.946) / 1200)

print upper
print lower

結果は

0.0667881551445
0.0412118448555
  • 全国から無作為に2400人を抽出
  • そのうち、1250人が政策に賛成
  • 有権者全体の半数が賛成していると言えるか

  • 帰無仮説は、

    • 半数が賛成している
  • 対立仮説は

    • 半数が賛成してはいない(半数以下)
  • 母集団の比率を0.5として、1250 / 2400の割合で賛成がいるのがどの程度レアなのか調べる

  • -1.96 <= ( 0.054 - p ) / √((0.054 * 0.946) / 1200) <= 1.96
import math

sample_p = (1250 / 2400)

z = (sample_p - 0.5) / math.sqrt((sample_p * (1 - sample_p)) / 2400)

print z

結果は、

2.04
  • 95%で片側検定するなら、1.645以上なら棄却域に入る

2.04 > 1.645 なので、仮説は棄却

母分散の推定

光速の測定値は以下の通り。

850 740 900 1070 930 850 950 980 980 880 1000 980 930 650 760 810 1000 960 960

この結果から、光速の分散を区間推定する。

カイ2乗値

  • 次の値は自由度n - 1のカイ2乗分布する
    • 偏差の2乗和を分散で除算した値
  • 自由度25 -1で95%信頼区間の場合、
    • 下限は8.91、αは0.975、カイ2乗分布のグラフの面積の97.5%
    • 上限は32.85、αは0.025、カイ2乗分布のグラフの面積の2.5%
  • つまり、8.91 <= カイ2乗値 <= 32.85 になる母分散が区間になる。
data = [299850, 299740, 299900, 300070, 299930,
        299850, 299950, 299980, 299980, 299880,
        300000, 299980, 299930, 299650, 299760,
        299810, 300000, 300000, 299960, 299960]

# 標本平均
sample_mean = sum(data) / len(data)

# 偏差の2乗和
squared_deviation = sum(map(lambda x: (x - sample_mean) * (x - sample_mean), data))

# 8.91 <= カイ2乗値 <= 32.85
# カイ2乗値 = 偏差の2乗和 / 母分散
# 母分散 = 偏差の2乗和 / カイ2乗値
chisquared_upper = squared_deviation / 8.91
chisquared_lower = squared_deviation / 32.85

print chisquared_lower
print chisquared_upper

結果は、

6367.73211568
23476.9921437

母平均の推定/検定

区間推定

光速の測定値は以下の通り。

850
740
900
1070
930
850
950
980
980
880
1000
980
930
650
760
810
1000
960
960

この結果から、光速を区間推定する。信頼区間は95%。

  • 母集団の分散がわからない場合は、
    • 標本の平均
    • 標本の普遍分散 を使って、推定する。
import math

data = [850, 740, 900, 1070, 930, 850, 950, 980, 980, 880, 1000, 980,
        930, 650, 760, 810, 1000, 960, 960]

# 平均
sample_mean = sum(data) / len(data)

# 普遍分散
deviation =sum(map(lambda n: (n - sample_mean) * (n - sample_mean), data)) / (len(data) - 1)

# t = 標本平均 - 母平均 / (√不変分散/n) 
# 標本平均 = sample_mean, 不変分散 = deviation, n = len(data)  - 1 
# -2.09 <= t <= 2.09 となる 母平均の範囲を求める

# まず下限
# t = (sample_mean - 母平均) / (√deviation/len(data)) = 2.09 だから
mean_lower = sample_mean - 2.09 * (math.sqrt(deviation / len(data) - 1))

# まず上限
# t = (sample_mean - 母平均) / (√deviation/len(data)) = -2.09 だから
mean_upper = sample_mean + 2.09 * (math.sqrt(deviation / len(data) - 1))

print 'lower ' + str(mean_lower)
print 'upper ' + str(mean_upper)

結果は

lower 853.449643918
upper 954.550356082

検定

  • 新入生 3000人
  • 毎年試験
  • 今年の新入生から36人をサンプリングしたら、平均点は480点
  • 今年は例年よりも賢いのか

 z検定

  • これまでの平均は450点
  • これまでの標準偏差は80点

帰無仮説を「今年も例年通りで学力に変化なし」 = 平均は450 対立仮説を「今年は変化があった」 = 平均は450ではない

とする。

  • 変化があったかどうかを検定するので、両側検定
  • 有意水準5パーセント
  • 従来の試験得点は平均450 分散80の正規分布
  • この正規分布を前提にした場合、36人の平均値である480はどの程度レアなのか
  • あまりにもレアなら、変化があったと判断できる
z = (480 - 450) / (80 / 6)

結果は、2.25なので、1.96よりも大きい。棄却域に入るので、帰無仮説を棄却する。つまり、「変化があった」ということ。

 t検定

  • パン屋チェーン店のクリームパンの重さは105g
  • しかし、A支店はもっと重くしている疑いがある
  • A支店のクリームパン26個の重さを調べたところ

  • 帰無仮説

    • A支店のクリームパンの重さは平均したて105gである。
  • t値は、2.54

import math
t = (106 - 105) / (2 / math.sqrt(26))
  • 「重くしている可能性がある」ので右側片側検定
  • 自由度25で5%に対するt値は、1.708
  • 2.54の方が大きいので、仮説は棄却する。

回帰分析について

以下のデータの公転周期=orbital-periodを応答変数、説明変数をsemi-major-axisにして分析する。

planet,inclination,eccentricity,semi-major-axis,orbital-period
Mercury,7.004,0.2056,0.3871,0.241
Venus,3.39471,0.0068,0.72333,0.615
Earth,0.00005,0.0167,1,1
Mars,1.85061,0.0934,1.52366,1.881
Jupiter,1.3053,0.0484,5.20336,11.86
Saturn,2.48446,0.0542,9.53707,29.46
Uranus,0.774,0.0461,19.19138,84.01
Neptune,1.76917,0.0086,30.06896,164.79
Pluto,17.089,0.25025,39.445,247.74
import math
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

data = pd.read_csv("/planet.data")

axismean = data['semi-major-axis'].mean()

periodmean = data['orbital-period'].mean()

axismean_var = data.apply(lambda x: (x["semi-major-axis"] - axismean) * (x["semi-major-axis"] - axismean) , axis=1).sum() / data['semi-major-axis'].count() - 1

periodmean_var = data.apply(lambda x: (x["orbital-period"] - periodmean) * (x["orbital-period"] - periodmean) , axis=1).sum() / data['orbital-period'].count() - 1

covariance = data.apply(lambda x: (x["semi-major-axis"] - axismean) * (x["orbital-period"] - periodmean), axis=1).sum() / data['orbital-period'].count() - 1

correlation_coefficient = covariance / (math.sqrt(axismean_var) * math.sqrt(periodmean_var))

regression_coefficient = covariance / axismean_var

alpha = periodmean - regression_coefficient * axismean

plt.scatter(data['semi-major-axis'], data['orbital-period'])

r = np.corrcoef(data['semi-major-axis'], data['orbital-period'])

x = np.linspace(0, 50)
y = alpha + regression_coefficient * x
plt.plot(x, y, "r-")
plt.title(" alpha = " + str(round(alpha,2)) + " regression_coefficient = " + str(round(regression_coefficient,2)) + " correlation_coefficient = " + str(round(correlation_coefficient,2)))
plt.show()

結果

f:id:stokutake:20151025130147p:plain