×
[PR]上記の広告は3ヶ月以上新規記事投稿のないブログに表示されています。新しい記事を書く事で広告が消えます。
確率分布は基本的にRjpwikiのまとめを見ればわかるが、少しだけ覚書しておく
これも便利。感謝
標準正規分布で特定の値の有意水準を調べる
z得点が1.96のとき
pnorm(1.96)
[1] 0.9750021
有意水準としては (zが2だけど、これは棄却閾5%に含まれる？)
1-pnorm(2)
[1] 0.02275013 # OK
逆に、有意水準5%以下となるzの値を調べるには
qnorm(0.975)
[1] 1.959964
正規分布曲線を描く
x <- seq(-4, 4, 0.01) # 最小値-4, 最大値4、0.01刻みのベクトル
curve(dnorm(x, 0, 1), -4, 4) # ベクトルxについて、平均0、標準偏差1の曲線を描く
確率分布関数の覚書
d... distribution。分布関数曲線の数式。確率密度関数のこと
p... percentile 。特定のポイントでの面積 (積分) 。よく斜線が引いてあるアレ
q... quantile。ある確率密度のときの横軸 (標準得点の値) 。よく1.96とか使われるアレ
r... random。特定の分布に基づいて乱数を発生させる (例えば正規分布に基づく乱数だと、平均0、標準偏差1になる)
これも便利。感謝
標準正規分布で特定の値の有意水準を調べる
z得点が1.96のとき
pnorm(1.96)
[1] 0.9750021
有意水準としては (zが2だけど、これは棄却閾5%に含まれる？)
1-pnorm(2)
[1] 0.02275013 # OK
逆に、有意水準5%以下となるzの値を調べるには
qnorm(0.975)
[1] 1.959964
カイ二乗では
あるカイ二乗値の有意水準が知りたい (カイ二乗値 = 18.30704, df = 10)
pchisq(18.30704, 10)
[1] 0.95
自由度と有意水準がわかっていて、カイ二乗値を知りたいとき
qchisq(0.95, 10)
[1] 18.30704
あるカイ二乗値の有意水準が知りたい (カイ二乗値 = 18.30704, df = 10)
pchisq(18.30704, 10)
[1] 0.95
自由度と有意水準がわかっていて、カイ二乗値を知りたいとき
qchisq(0.95, 10)
[1] 18.30704
正規分布曲線を描く
x <- seq(-4, 4, 0.01) # 最小値-4, 最大値4、0.01刻みのベクトル
curve(dnorm(x, 0, 1), -4, 4) # ベクトルxについて、平均0、標準偏差1の曲線を描く
確率分布関数の覚書
d... distribution。分布関数曲線の数式。確率密度関数のこと
p... percentile 。特定のポイントでの面積 (積分) 。よく斜線が引いてあるアレ
q... quantile。ある確率密度のときの横軸 (標準得点の値) 。よく1.96とか使われるアレ
r... random。特定の分布に基づいて乱数を発生させる (例えば正規分布に基づく乱数だと、平均0、標準偏差1になる)
PR
母平均の信頼区間を求める。母平均0との検定を行う
t.test関数で両方わかって超便利
##データを用意
x <- c(1,2,3,4,5,6,7,8,9)
t.test(x, m=0)
t.test(x, m=0)
One Sample t-test
data: x
t = 5.4772, df = 8, p-value = 0.0005894
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
2.894916 7.105084
sample estimates:
mean of x
5
data: x
t = 5.4772, df = 8, p-value = 0.0005894
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
2.894916 7.105084
sample estimates:
mean of x
5
##これでいいのだが、信頼区間の覚書
95%信頼区間の加減と上限の絶対値は
x.qt <- abs(qt(0.025, 8))
[1] 2.306004
平均値の標準誤差
x.se <- sqrt(var(x)/length(x))
[1] 0.912871
mean(x)-(x.qt*x.se)
[1] 2.894916
以下、理屈の覚書
ある母集団からサンプリングした標本平均は、平均μ、分散σ2/nの正規分布に従う (これはすなわち標本平均のサンプリング分布のことであり、この分散の平方根は標本平均の標準誤差と同じ)
母集団分散σ2を標本分散 (本当は不偏分散) に置き換えて標準化する
(mean(x) - μ )/ se(x)
とこれが自由度n-1のt分布に従う。
これでqt関数を使用すればt分布の下限と上限がわかり、式を展開すると
下限値はmean(x)-(x.qt*x.se)
上限値はmean(x)+(x.qt*x.se)
母比率の区間推定
100サンプルで40のヒットが出た場合
n <- 100
p <- 40/100
比率の標準誤差
p.se <- sqrt((p*(1-p))/n))
これにqnorm(0.975)をかけ、標本比率にプラスマイナスすれば信頼区間が出る、と教科書には書いてある
p - qnorm(0.975)*p.se
[1] 0.3039818
Rでbinom.test関数を用いると信頼区間が算出され、特定の比率との検定ができる
binom.test(40, 100, 0.5)
Exact binomial test
data: 40 and 100
number of successes = 40, number of trials = 100, p-value = 0.05689
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
0.3032948 0.5027908
sample estimates:
probability of success
0.4
### この信頼区間は上の計算と一致しない。binom.test関数でどんな計算をしてるかわからないが、たぶんqnorm(0.975) をかけるやりかたはあくまでも正規分布への近似だからじゃないかと思う。
参考: http://oku.edu.mie-u.ac.jp/~okumura/stat/tests_and_CI.html
相関係数の標準誤差
x <- c(1,2,3,4,5,6)
y <- c(2,1,4,3,6,5)
r <- cor(x,y)
cor.se <- (1-r^2)/sqrt(6)
cor.test(x,y)
Pearson's product-moment correlation
data: x and y
t = 2.9598, df = 4, p-value = 0.04156
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.05192932 0.98068460
sample estimates:
cor
0.8285714
相関係数の信頼区間はフィッシャーのz変換とかして面倒だからパスしよう
最尤法、最尤推定
母数の分布はわかっており、母数が未知であるときに、標本データから母数を推定する手法
L = px(1-p)(n-x)
pは母集団比率、nはサンプリング数、xはデータ数、Lは確率
pxはそのデータが得られる確率を示し、これに得られない確率(1-p)(n-x)をかける
pは母集団比率であるが、これを確率変数とみなすと
pが0.001であるとき標本データが得られる確率、0.002であるときの…というふうに変動させることができる
この中で、標本データが得られる確率の最も高いものを、母数 (母集団比率) とみなす。
のが最尤法。ホンマか嘘かはわからんが、たぶん一番もっともらしい値を確率的に定義する
自然対数をとって微分するとpは容易に求められる
すなわち、px (および(1-p)(n-x)) のように、あるデータが得られる確率を規定できれば最尤推定が可能になる。
これが確率モデルと呼ばれるもので、事象が得られる (変動する) 確率を尤度 (もっともらしさ) といい、
標本データより算出された尤度の最も高いものを最尤推定値という
因子分析では、標本共分散行列の確率密度がウィッシャート分布に従うので、これと
因子分析モデル (係数x因子+独自因子) の確率変数との差を最小化することにより
特定の因子数モデルの尤度を最大化する
母数の分布はわかっており、母数が未知であるときに、標本データから母数を推定する手法
L = px(1-p)(n-x)
pは母集団比率、nはサンプリング数、xはデータ数、Lは確率
pxはそのデータが得られる確率を示し、これに得られない確率(1-p)(n-x)をかける
pは母集団比率であるが、これを確率変数とみなすと
pが0.001であるとき標本データが得られる確率、0.002であるときの…というふうに変動させることができる
この中で、標本データが得られる確率の最も高いものを、母数 (母集団比率) とみなす。
のが最尤法。ホンマか嘘かはわからんが、たぶん一番もっともらしい値を確率的に定義する
自然対数をとって微分するとpは容易に求められる
すなわち、px (および(1-p)(n-x)) のように、あるデータが得られる確率を規定できれば最尤推定が可能になる。
これが確率モデルと呼ばれるもので、事象が得られる (変動する) 確率を尤度 (もっともらしさ) といい、
標本データより算出された尤度の最も高いものを最尤推定値という
因子分析では、標本共分散行列の確率密度がウィッシャート分布に従うので、これと
因子分析モデル (係数x因子+独自因子) の確率変数との差を最小化することにより
特定の因子数モデルの尤度を最大化する
確率変数
ある確率で特定の値をとる変数
確率的に変動する変数。逆に言えば、確率的に既定される変数。
あるバッターがある特定の打席でヒットを打つかどうかは打率によって既定される (ように表現されている)
心理学は (大体) 推測統計をするから、どのデータも確率変数
サイコロの1-6はそれぞれ1/6の確率でとる確率変数
確率分布
確率変数がそれぞれどういう値をとるか決める分布。
サイコロの例ではp(x) = P(X=x) = 1/6。1/6という棒グラフが横から6つ並んでいるイメージ
※Σp(x) は常に1になる
確率密度関数 probability density function
サイコロの場合、1-6、それぞれの生起確率は1/6であるが、連続変量の場合、その取りうる値は無限である。この場合、確率密度関数を定義し、ある確率変数が任意の区間に含められる確率を求められるようにしておく
Rで標準正規分布
x <- seq(-3,3,0.01)
curve(dnorm(x,0,1), -3, 3, type="l")
調整平均
(x2+x3+x4...xn-1)/(n-2)
最小値と最大値を省いて平均をとったものが100%調整平均
母数
parameter. 母集団分布の特徴を表す数字。母集団平均、母集団相関とか
標本分布
sampling distribution. 標本の度数分布ではない。間違えそうなのでサンプリング分布と読んでみよう (南風原, 2003)
標本統計量 (観測データの平均とか) の変動を既定する分布
大学生、という母集団から研究のため20人サンプリングして知能検査をする。このときの平均が100だった。もう一回同じ実験をしたら110だった。もう一回やったら95だった...ということを何度もやったとき、平均の変動を決める分布。真ん中に集まってたら変動しないし、広かったら変動しやすいだろう
標準誤差
サンプリング分布の標準偏差。上述の真ん中に集まってたら、とはサンプリング分布の標準偏差が小さい、すなわち標準誤差が小さいという意味で、広かったらとは同様に大きいという意味。
Rでは sqrt(var(as.vector(x))/length(x))
推定量estimatorと推定値estimate
母数の点推定 (母集団の平均はこれだ！) に使う標本統計量。標本平均、標本比率…等
推定量は概念であり、これを適用して実際のデータを収集し求めた値を推定値という
大学生の平均知能指数 (母数) は100人くらい標本をとって平均をとれば推定できる (推定量) だろう。で、実際取ったら平均は108だった (推定値)
標準化
粗点から平均を引いて標準偏差で割る。このときの標準偏差は不偏標準偏差。
あくまでも母集団標準偏差の代わりとして使うから
Rではscale(x)
ある確率で特定の値をとる変数
確率的に変動する変数。逆に言えば、確率的に既定される変数。
あるバッターがある特定の打席でヒットを打つかどうかは打率によって既定される (ように表現されている)
心理学は (大体) 推測統計をするから、どのデータも確率変数
サイコロの1-6はそれぞれ1/6の確率でとる確率変数
確率分布
確率変数がそれぞれどういう値をとるか決める分布。
サイコロの例ではp(x) = P(X=x) = 1/6。1/6という棒グラフが横から6つ並んでいるイメージ
※Σp(x) は常に1になる
確率密度関数 probability density function
サイコロの場合、1-6、それぞれの生起確率は1/6であるが、連続変量の場合、その取りうる値は無限である。この場合、確率密度関数を定義し、ある確率変数が任意の区間に含められる確率を求められるようにしておく
Rで標準正規分布
x <- seq(-3,3,0.01)
curve(dnorm(x,0,1), -3, 3, type="l")
調整平均
(x2+x3+x4...xn-1)/(n-2)
最小値と最大値を省いて平均をとったものが100%調整平均
母数
parameter. 母集団分布の特徴を表す数字。母集団平均、母集団相関とか
標本分布
sampling distribution. 標本の度数分布ではない。間違えそうなのでサンプリング分布と読んでみよう (南風原, 2003)
標本統計量 (観測データの平均とか) の変動を既定する分布
大学生、という母集団から研究のため20人サンプリングして知能検査をする。このときの平均が100だった。もう一回同じ実験をしたら110だった。もう一回やったら95だった...ということを何度もやったとき、平均の変動を決める分布。真ん中に集まってたら変動しないし、広かったら変動しやすいだろう
標準誤差
サンプリング分布の標準偏差。上述の真ん中に集まってたら、とはサンプリング分布の標準偏差が小さい、すなわち標準誤差が小さいという意味で、広かったらとは同様に大きいという意味。
Rでは sqrt(var(as.vector(x))/length(x))
推定量estimatorと推定値estimate
母数の点推定 (母集団の平均はこれだ！) に使う標本統計量。標本平均、標本比率…等
推定量は概念であり、これを適用して実際のデータを収集し求めた値を推定値という
大学生の平均知能指数 (母数) は100人くらい標本をとって平均をとれば推定できる (推定量) だろう。で、実際取ったら平均は108だった (推定値)
標準化
粗点から平均を引いて標準偏差で割る。このときの標準偏差は不偏標準偏差。
あくまでも母集団標準偏差の代わりとして使うから
Rではscale(x)
Q1
イチローは現在のところ、打率.365である
これにより、1打数でヒットを打つ確率が.365とする
イチローがある試合で5打数3安打となる確率はどのくらいか
5打数3安打となるケース
1，2，3打席でヒットを打ち、残り2打席で凡退…ケース1
1，2, 4打席でヒットをうち、3打席目と4打席目で凡退…ケース2
…
ケースとしては5C3 = 10とおりある
ケース1の確率を計算すると
.365x.365x.365x(1-.365)x(1-.365) = 0.01960...となる
ケース2は
.365x.365x(1-.365)x.365x(1-.365) = 0.01960...で同じ
全てのケースを総和する (= 5C3 = 10をかけると)
.196となる
A. イチローがある試合で5打数3安打となる確率は.196である
※ 4安打や5安打の場合は考えない
逆に考えてみる
Q2
1打数でヒットを打つ確率がどのくらいあれば、最も5打数3安打となりやすいか (5打数3安打の結果となるバッターは打率何割のバッターか)
Q1の計算を色々な打率について行ってみる
打率.133 -> .017
打率.365 -> .196
打率.900 -> .073
式はこうなる: (打率^3)x((1-打率)^2)x10
このデータをプロットする
# Rのコード optim関数の例
A.並べてみれば明らかだが、最も5打数3安打となりやすいのは1打数のヒットの確率が.600のとき (打率.600のバッター)
このバッターが5打数3安打になる確率は.346
この.346も含めた、グラフの縦軸を尤度と呼ぶ
定義的にはある確率モデルがあるとき、ある事象が得られる度合い
上の例で言えば、確率モデルとは(打率^3)x((1-打率)^2)x10、ある事象とは「ある試合での5打数3安打」
別の例で言えば、
コインを投げて、9回表、9回裏がでる度合いを調べるのも尤度推定
尤度が最大 (.346) になるパラメータ (打率.600) を求めるのが最尤法で、微分すると簡単らしい
わかったようなわからんような話だな…
イチローは現在のところ、打率.365である
これにより、1打数でヒットを打つ確率が.365とする
イチローがある試合で5打数3安打となる確率はどのくらいか
5打数3安打となるケース
1，2，3打席でヒットを打ち、残り2打席で凡退…ケース1
1，2, 4打席でヒットをうち、3打席目と4打席目で凡退…ケース2
…
ケースとしては5C3 = 10とおりある
ケース1の確率を計算すると
.365x.365x.365x(1-.365)x(1-.365) = 0.01960...となる
ケース2は
.365x.365x(1-.365)x.365x(1-.365) = 0.01960...で同じ
全てのケースを総和する (= 5C3 = 10をかけると)
.196となる
A. イチローがある試合で5打数3安打となる確率は.196である
※ 4安打や5安打の場合は考えない
逆に考えてみる
Q2
1打数でヒットを打つ確率がどのくらいあれば、最も5打数3安打となりやすいか (5打数3安打の結果となるバッターは打率何割のバッターか)
Q1の計算を色々な打率について行ってみる
打率.133 -> .017
打率.365 -> .196
打率.900 -> .073
式はこうなる: (打率^3)x((1-打率)^2)x10
このデータをプロットする
# Rのコード optim関数の例
xax <- seq(0.001, 0.999, length=999)
yax <- (xax^3)*((1-xax)^2)*10
par(family="Japan1GothicBBB")
plot(xax, yax, xlab="打率", ylab="5打数3安打の確率")
yax <- (xax^3)*((1-xax)^2)*10
par(family="Japan1GothicBBB")
plot(xax, yax, xlab="打率", ylab="5打数3安打の確率")
A.並べてみれば明らかだが、最も5打数3安打となりやすいのは1打数のヒットの確率が.600のとき (打率.600のバッター)
このバッターが5打数3安打になる確率は.346
この.346も含めた、グラフの縦軸を尤度と呼ぶ
定義的にはある確率モデルがあるとき、ある事象が得られる度合い
上の例で言えば、確率モデルとは(打率^3)x((1-打率)^2)x10、ある事象とは「ある試合での5打数3安打」
別の例で言えば、
コインを投げて、9回表、9回裏がでる度合いを調べるのも尤度推定
尤度が最大 (.346) になるパラメータ (打率.600) を求めるのが最尤法で、微分すると簡単らしい
わかったようなわからんような話だな…
