2012年01月11日

極値統計学の勉強

午前中、お金の清算とボスの書いた文章のチェック。
午後、今学会用に解析を進めているが、極値統計の理論が使えるのではないかと思い極値統計のお勉強。とりあえず最大値の分布を作って、それに従う分布のパラメータを推定するところまで。まだ使えるかどうかは分からない。とりあえず。

参考資料1(確率密度分布の式が間違っているような…)
参考資料2

################
###極値統計の勉強###
################
set.seed(1)
a <- matrix(rnorm(1e+6), 1000, 1000)

###それぞれの最大値だけ抜き取る###
a_max <- apply(a, 2, max)
ind <- seq(1, 5, length=100)
b <- matrix(NA, (length(ind)-1), 4)

###それぞれの区間に入る数をカウントする###
for(i in 1 : (length(ind)-1)){
temp <- subset(a_max, a_max >= ind[i]&
a_max < ind[i+1])

b[i, 1] <- mean(ind[i], ind[i+1])
b[i, 2] <- length(temp)
}

###平均ランク法というもので密度に直す###
b[, 3] <- b[, 2]/(sum(b[, 2])+1)

###累積和を求める###
b[, 4] <- cumsum(b[, 3])

###余計な区間は飛ばす###
b2 <- subset(b, b[, 2]>0)

#plot(b2[, 1], b2[, 3])
#plot(b2[, 1], b2[, 4])
#plot(b2[, 1], -log(-log(b2[, 4])))

###直線回帰の準備###
y <- -log(-log(b2[, 4]))
x <- b2[, 1]

###欲しかったパラメタ(平均と分散を求めた)###
alpha <- 1/coef(lm(y ~ x))[2]
lambda <- coef(lm(y ~ x))[1]/(-1*coef(lm(y ~ x))[2])

###累積分布を確認###
jpeg("Ruiseki.jpeg")
x2 <- seq(1, 5, length=100)
plot(x2, exp(-exp(-(x2-lambda)/alpha)), type="l", lwd=2, col=2)
points(b2[, 1], b2[, 4], type="l", lwd=2, col=4)
legend(4, 0.8, c("Riron", "Jissoku"), col=c(2, 4), lwd=2)
dev.off()

###密度分布で確認###
jpeg("Mitudo.jpeg")
z <- (x2-lambda)/alpha
y2 <- (1/alpha)*exp(-z)*exp(-exp(-z))
plot(x2, y2, type="l", lwd=2, col=2, xlim=c(2, 5))
points(density(a_max), type="l", lwd=2, col=4)
legend(4, 1.2, c("Riron", "Jissoku"), col=c(2, 4), lwd=2)
dev.off()
################
###極値統計の勉強###
################


累積分布関数
Ruiseki.jpeg

確率密度関数
Mitudo.jpeg
ラベル:統計学 R
posted by しばきん at 16:33| 神奈川 ☀| Comment(0) | TrackBack(0) | 統計 | このブログの読者になる | 更新情報をチェックする

2011年10月31日

n次元確率ベクトルの1次関数の平均と分散

いやはや、数理統計学でここ一週間ずうっと悩んでいたものがようやく解決した。
やはり、独立というものをきちんと理解しておかないといけない。

a,Xがそれぞれベクトルと行列だとして、cov(Xi, Xj)=0のとき(X1,...,Xnが互いに独立のとき)、
var(aX) = a1^2*var(X1) + a2^2*var(X2) +, ..., + an^2*var(Xn)
となります。これは、a*var(X)'*aと一致するはずなのですが、一致しません。なんでだろう、と思って色々やってみていたのだけれど、要は平均すれば一致する、ということのようです。シミュレーションしてみました。

> set.seed(1)
> n <- 100
> S <- 40
> a <- seq(1, n)
> X <- X2 <- X3 <- NULL
> temp <- temp2 <- numeric(s)
> for (s in 1 : S) {
+ X <- matrix(runif(n^2), nrow=n, ncol=n)
+ X2 <- matrix(NA, nrow=n, ncol=n)
+ X3 <- matrix(0, nrow=n, ncol=n)
+
+ ###期待値ベクトル###
+ ###ほぼ0だが、微妙に差が出てしまうようだ###
+ ###aが横ベクトルであることに注意###
+ round(apply(a*X, 1, mean)-
+ a*apply(X, 1, mean))
+
+ ###共分散行列###
+ for(j in 1 : ncol(X)){
+ for (i in 1 : nrow(X)) {
+ X2[j, i] <- cov(X[j, ], X[i, ])
+ }}
+ ###対称行列となっている###
+ #X2==t(X2)
+ #var(X)==t(var(X))
+ #X2==var(t(X))
+
+ ###独立とは言っても一回一回はばらつく###
+ temp[s] <- a%*%X2%*%a-
+ (a^2)%*%apply(X, 1, var)
+
+ ###対角だけ取ってきて無理矢理独立にすれば(?)いつでも差は0###
+ diag(X3) <- diag(X2)
+ temp2[s] <- a%*%X3%*%a-
+ (a^2)%*%diag(X3)
+ }
> jpeg("variance_independent.jpeg")
> hist(temp)
> dev.off()
windows
2
> temp;temp2
[1] -173.4482 1113.3236 -7706.5664 502.3490 -7383.8686 -993.2250
[7] 1724.8721 2393.0135 5236.6691 3939.5105 827.2226 -1489.3155
[13] -2922.6024 447.2575 -1754.0834 -1780.8433 1820.3705 -10123.5030
[19] -2082.8554 2079.5110 -3439.7196 -2314.1255 -5338.1769 6276.1852
[25] -2316.7343 1050.0900 -2247.3146 -3779.2307 8175.3117 1395.5665
[31] -3659.9512 -2296.1468 -7303.6068 116.1688 -1407.6821 -4242.2000
[37] -1149.0380 6796.3341 952.2544 797.3517
[1] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 -3.637979e-12
[6] -3.637979e-12 0.000000e+00 0.000000e+00 -3.637979e-12 0.000000e+00
[11] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[16] 0.000000e+00 3.637979e-12 0.000000e+00 0.000000e+00 0.000000e+00
[21] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 -3.637979e-12
[26] 0.000000e+00 -3.637979e-12 0.000000e+00 3.637979e-12 -3.637979e-12
[31] 0.000000e+00 3.637979e-12 0.000000e+00 7.275958e-12 0.000000e+00
[36] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00

出来た図。
variance_independent.jpeg
確かに、何度も繰り返してやれば一応イコールで結ばれるようです。なるほどなぁ。独立って難しい。
posted by しばきん at 00:17| 神奈川 ☔| Comment(0) | TrackBack(0) | 統計 | このブログの読者になる | 更新情報をチェックする

2011年07月26日

二項分布からのポアソン分布の導出と幾何分布のお勉強と確率母関数からの平均と分散の導出

以前も二項分布ポアソン周りは計算したのだけれど、もう一度おさらい。それから幾何分布の勉強と、確率母関数の微分を使った平均値と分散の導出の練習。微分で平均値とか分散が出るなんて不思議。すごいな、数理統計。感動する。

> ###強度の決定###
> lambda <- 3
>
> ###適当な数値列と階乗計算###
> x <- seq(0, 10)
> x2 <- NULL
> for (i in 1 : length(x)){
+ x2[i] <- gamma(x[i]+1)
+ }
>
> ###ポアソン分布###
> fx_lambda <- ((lambda^x)/x2)*exp(-1*lambda)
>
> ###積分して1になることを確認(まぁ10までなのできちんと1にならない)###
> sum(fx_lambda)
[1] 0.9997077
>
> ###プロットして確認###
> #pdf("Poissons.pdf")
> plot(x, fx_lambda, main="Poisson distributions")
>
> ###既成関数とも一致###
> points(x, dpois(lambda=lambda, x), cex=1.3, pch=0)
>
> ###二項分布からのポアソン分布導出###
> prob <- 0.000001
> size <- lambda/prob
> points(x, dbinom(x=x, size=size, prob=prob), col=2, cex=2)
> legend(7, 0.2, c("Po_hand", "Po_R", "Po_bin"), col=c(1, 1, 2), pch=c(1, 0, 1))
> #dev.off()

Poissons.jpeg

ここから幾何分布と確率母関数を微分してからの平均値と分散の導出。

> x <- seq(0, 20)
> p <- 0.1
> fx_p <- p*((1-p)^x)
>
> ###確認と既成関数###
> plot(x, fx_p)
> points(x, dgeom(x=x, prob=p), cex=2, col=2)
> legend(15, 0.1, c("G_hand", "G_R"), col=c(1, 2), pch=1)
>
> ###幾何分布の確率母関数###
> func_g <- function(t, p) {
+ p/(1-t*(1-p))
+ }
>
> ###一階微分###
> func_g_d <- deriv(~ p/(1-t*(1-p)), c("t"))
> func_g_d1 <- function(t, p){}
> body(func_g_d1) <- func_g_d
>
> ###微分からの平均値###
> attr(func_g_d1(t=1, p=p),"gradient")
t
[1,] 9
>
> ###教科書に書かれている平均値の定義からの計算と一致###
> (1-p)/p
[1] 9
>
> ###高階微分のための関数(R-tipsより)###
> DD <- function(expr, name, order = 1) {
+ if(order < 1) stop("'order' must be >= 1")
+ if(order == 1) D(expr, name)
+ else DD(D(expr, name), name, order - 1)
+ }
> func_g_d_e <- expression(p/(1-t*(1-p)))
> func <- function(t, p) {}
> body(func) <- DD(func_g_d_e, "t", 2)
>
> ###二階微分による分散###
> func(t=1, p=p) + (1-p)/p -((1-p)/p)^2
[1] 90
>
> ###教科書に書かれている分散の定義からの計算と一致###
> (1-p)/(p^2)
[1] 90
ラベル:統計学 数学 微分 R
posted by しばきん at 09:15| 神奈川 ☁| Comment(0) | TrackBack(0) | 統計 | このブログの読者になる | 更新情報をチェックする

2011年07月17日

Fisherの正確確率検定とオッズ比と超幾何分布と二項分布

超幾何分布を勉強していたのに、いつのまにかFisherの正確確率検定の勉強になっていました。そういえば自力で計算したことがなかったなぁと思ったので、理解できているか確認してみます。
> ###########################
> ###フィッシャーの正確確率検定###
> ###########################
>
> ###2*2分割表の各値の定義###
> a <- 13
> b <- 5
> c <- 6
> d <- 8
>
> ###ひとまとめ###
> x <- matrix(c(a, b, c, d), byrow=TRUE, nc=2)
>
> ###各段の合計を計算###
> e <- sum(x[1, ])
> f <- sum(x[2, ])
> g <- sum(x[, 1])
> h <- sum(x[, 2])
> n <- sum(e,f)
>
> ###正確確率を求めるために全通りを計算してみる###
> a2 <- 0:e
> b2 <- e-a2
> c2 <- g-a2
> d2 <- h-b2 #f-c2でも良い
>
> ###これが各段の合計を固定した時のセル内の値が取りうる全通り###
> com <- data.frame(a2, b2, c2, d2)
>
> ###マイナスの値は省く###
> temp <- {com[, 1:4]<0}+0
> com2 <- data.frame(com[apply(temp, 1, sum)==0,])
>
> com2[, 5] <-
+
+ ###一段目に関して、a2が起こる確率###
+ choose(e,com2[, 1])*
+
+ ###二段目に関して、c2が起こる確率###
+ choose(f,com2[, 3])/
+
+ ###合計の数から、gが起こる確率###
+ choose(n,g)
>
> ###合計すると1になることを確認###
> sum(com2[, 5])
[1] 1
>
> ###確率を求めるために、今回得られた値より小さなp値を全て合計する###
> p <- sum(com2[which(com2[,1]==a),5], com2[com2[which(com2[,1]==a),5] > com2[,5], 5])
>
> ###得られたp値###
> p
[1] 0.1491833
>
> ###Rの関数とも出力が一致###
> fisher.test(x, B = 500)

Fisher's Exact Test for Count Data

data: x
p-value = 0.1492
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
0.6373828 19.6279452
sample estimates:
odds ratio
3.324929

>
> ###オッズも計算###
> ###簡略化###
> a*d/(b*c)
[1] 3.466667
>
> ###行方向でも列方向でも一緒の値になる###
> (a/b)/(c/d)
[1] 3.466667
> (a/c)/(b/d)
[1] 3.466667
>
> ###オッズ比の信頼区間###
> (a*d/(b*c))*exp(c(1, -1)*qnorm(0.025)*sqrt(1/a+1/b+1/c+1/d))
[1] 0.790796 15.197063
>
> ###########################
> ###フィッシャーの正確確率検定###
> ###########################

本当はオッズ比を超幾何分布を最尤推定して求めたかったのですけどよくわからなかったので、数理統計の理論と応用という参考文献を使って推定してみました。

それから超幾何分布は二項分布で近似できるらしいので、両方をグラフ化して確認してみます。確かにnが多いと線が重なります。図はnが少ないときです。
超幾何分布と二項分布の図
 ####################
> ###超幾何分布と二項分布###
> ####################
>
> ###故障なし製品の数を例に###
> ###ここの数が多いと、二項分布に近似されていく###
> n <- 80
>
> ####不良品の数##
> m <- 10
>
> ###サンプル数###
> k <- 5
>
> ###サンプル数に含まれる不良品の数###
> x <- 0:k
>
> ###超幾何分布を手計算###
> temp <- (choose(m, x)*choose(n, k-x))/choose(m+n, k)
>
> ###確認###
> #pdf("hyper_geometric.pdf")
> plot(dhyper(x=x, n=n, m=m, k=k), type="l", col=2, lwd=2)
> points(dbinom(x=x, size=k, prob=m/n), type="l",col=3, lwd=2)
> points(temp, pch=16)
> legend(4, 0.5, c("hypergeo", "binom"), col=c(2, 3), lwd=2)
> #dev.off()
> ####################
> ###超幾何分布と二項分布###
> ####################
ラベル:統計学 数学 R
posted by しばきん at 18:37| 神奈川 ☀| Comment(0) | TrackBack(0) | 統計 | このブログの読者になる | 更新情報をチェックする

2011年07月12日

log()にゲタをはかせることについて

前のエントリーで中学時代の友人からlog(x+1)でもp-valueはあんまりかわんねーぞ!と言われたので試しにやってみた。

> set.seed(1)
> x <- c(rep(0, 100), seq(1, 62)) #疑似共変量の生成
> y <- rpois(162, exp(0.01*x + 0.005)) #疑似応答変数の生成
> res1 <- lm(log(y+1) ~ x) #まずは+1というゲタをはかせて解析
> anova(res1)
Analysis of Variance Table

Response: log(y + 1)
Df Sum Sq Mean Sq F value Pr(>F)
x 1 1.027 1.02663 4.6504 0.03254 * #ゲタ+1で有意。
Residuals 160 35.322 0.22076
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> res2 <- lm(log(y+0.1) ~ x) #次にゲタ+0.1をはかせてみる
> anova(res2)
Analysis of Variance Table

Response: log(y + 0.1)
Df Sum Sq Mean Sq F value Pr(>F)
x 1 6.593 6.5930 4.0792 0.04508 * #まだ有意だけどp値は変わった。
Residuals 160 258.600 1.6162
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> res3 <- lm(log(y+0.01) ~ x) #+0.01
> anova(res3)
Analysis of Variance Table

Response: log(y + 0.01)
Df Sum Sq Mean Sq F value Pr(>F)
x 1 18.78 18.7803 3.6879 0.05659 . #有意じゃなくなった。
Residuals 160 814.78 5.0924
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> res6 <- glm(y ~ x, family="poisson") #gaussianで解析しないでpoissonで。
> summary(res6)

Call:
glm(formula = y ~ x, family = "poisson")

Deviance Residuals:
Min 1Q Median 3Q Max
-1.78259 -1.43247 -0.02576 0.47585 2.58039

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.025653 0.091217 0.281 0.7785
x 0.007540 0.003581 2.105 0.0353 * #poissonでやればゲタなくても有意
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 169.05 on 161 degrees of freedom
Residual deviance: 164.85 on 160 degrees of freedom
AIC: 436.92

Number of Fisher Scoring iterations: 5

上記のサンプルからわかるように、ゲタのはかせかたで結果は変わってしまうので、それなら最初からそのデータにあった解析を見つけてあげるのが良いと思うわけです。なんか勘違いしてたらご指摘ください。
ラベル:統計学 R
posted by しばきん at 08:47| 神奈川 ☀| Comment(0) | TrackBack(0) | 統計 | このブログの読者になる | 更新情報をチェックする
×

この広告は1年以上新しい記事の投稿がないブログに表示されております。