2010年10月31日

標準誤差と標準偏差

中心極限定理から、xが平均μ,標準偏差σのある分布にしたがうとき,大きさnの無作為標本に基づく標本平均 は,nが無限に大きくなるとき,平均μ,標準偏差σ/sqrt(n)の正規分布に近づく。というお話が色んな教科書でされています。そういえば一度も自分で確認していなかったので、ちょっとシミュレーションをしてみました。


###母集団から何回抽出するか(本当は無限回試行にしたい)###
iter <- 3000

###母集団の標準偏差###
sigma <- 10

###母集団からいくつ抽出するか###
u <- c(seq(10, 1000, length=50))

###入れ物の用意###
a <- matrix(0, iter, length(u))

for(j in 1 : length(u)) {
for (i in 1 : iter) {
x <- rnorm(u[j], 0, sigma)
a[i, j] <- mean(x)
}
}
###中心極限定理からの期待値と合っているか確認###
#pdf("SE_SD_sqrt_sigma_.pdf")
plot(sigma/sqrt(u), apply(a, 2, sd), xlab="expected", ylab="Observation")
segments(0, 0, 4, 4, col="red")
#dev.off()

出来た図。ものすごい当てはまりよう。
SE_SD_sqrt_sigma_.pdf
ラベル:R 統計学
posted by しばきん at 17:25| 神奈川 ☁| Comment(0) | TrackBack(0) | 統計 | このブログの読者になる | 更新情報をチェックする

2010年10月29日

サンプリング方式の違いによる種の分布の推定結果に与える影響(Mohler, 1983)を検証してみた

サンプリングする場所によって多少なりとも結果が変わるのは、ランダムサンプリングならばよくあること(っていうかそれを許容するしかないし、無限回試行の後には不偏)。ならば、バイアスサンプリングによってならばどうか?まんべんなくサンプリングするよりも良い場合があるのではないか?と思ったかどうかは知らないが、そういうことをやっている著者がいたので、追認試験をやってみた。全部やるのは面倒だったので1パターンだけ(下記コードでbの中身を変えれば、実は他のパターンはできる)。汚いコードだけど、まぁよし。結果としては、ランダムよりもバイアスをかけてサンプリングするほうが、推定されたパラメータはより真の値と近かった。

参照論文
Effect of sampling pattern on estimation of species distributions along gradients


set.seed(1)
###真の価の設定。ここを推定する###
a <- 2
b <- c(0)
c <- 40

###繰り返し数の設定###
iter <- 2000

###入れ物の用意###
mat_rand <- mat_con <- mat_edge1 <- mat_edge2 <- matrix(0, iter*length(b), 3)
x <- seq(0, 100, by=0.1)

###真の関係の用意###
###xが水温だとすると、どの水温にどのくらい生物がいるのかわかるような関係###
#for (j in 0:1){
j <- 0
for(i in 1 : iter) {
f <- function(x, n, s){
a*exp(-(x-b[j+1])^2/(2*c^2)) + rnorm(n, 0, s)
}
#plot(x, f(x=x, n=length(x), s=0.1))

###サンプリングパターンの決定###
###ランダムに水温を抽出###
pat_random <- runif(30, 1, 100)

###中心付近を重点的にサンプリング###
pat_con <- c(runif(10, 40, 60), runif(20, 1, 100))

###端っこを重点的に###
pat_edge1 <- c(runif(10, 0, 20), runif(10, 20, 80), runif(10, 80, 100))

###少しゆるめに端っこをサンプリング###
pat_edge2 <- c(runif(10, 0, 30), runif(10, 30, 70), runif(10, 70, 100))

###誤差を伴ってサンプルが得られるはず###
D_rand <- data.frame(x=pat_random, y=f(x=pat_random, n=30, s=0.1))
D_con <- data.frame(x=pat_con, y=f(x=pat_con, n=30, s=0.1))
D_edge1 <- data.frame(x=pat_edge1, y=f(x=pat_edge1, n=30, s=0.1))
D_edge2 <- data.frame(x=pat_edge2, y=f(x=pat_edge2, n=30, s=0.1))

###最適化する関数の用意###
F <- function(par, x, y){
A <- par[1]
B <- par[2]
C <- par[3]
sum((y - A*exp(-(x-B)^2/(2*C^2)))^2)
}

###結果をばしばし入れていく###
mat_rand[i+(j*iter), 1:3] <- optim(par=c(1, 1, 40), F, x=D_rand$x, y=D_rand$y)$par
mat_con[i+(j*iter), 1:3] <- optim(par=c(1, 1, 40), F, x=D_con$x, y=D_con$y)$par
mat_edge1[i+(j*iter), 1:3] <- optim(par=c(1, 1, 40), F, x=D_edge1$x, y=D_edge1$y)$par
mat_edge2[i+(j*iter), 1:3] <- optim(par=c(1, 1, 40), F, x=D_edge2$x, y=D_edge2$y)$par
#}
}
#mat_rand <- mat_rand[mat_rand[, 1]<5,]
#mat_con <- mat_con[mat_con[, 1]<5,]
#mat_edge1 <- mat_edge1[mat_edge1[, 1]<5,]
#mat_edge2 <- mat_edge2[mat_edge2[, 1]<5,]

###結果一覧###
apply(mat_rand[1:iter,], 2, mean);apply(mat_rand[1:iter,], 2, sd)
apply(mat_con[1:iter,], 2, mean);apply(mat_con[1:iter,], 2, sd)
apply(mat_edge1[1:iter,], 2, mean);apply(mat_edge1[1:iter,], 2, sd)
apply(mat_edge2[1:iter,], 2, mean);apply(mat_edge2[1:iter,], 2, sd)

###パラメタごとに整理###
Estimated_a <-
data.frame(
mean = c(apply(mat_rand[1:iter,], 2, mean)[1],
apply(mat_con[1:iter,], 2, mean)[1],
apply(mat_edge1[1:iter,], 2, mean)[1],
apply(mat_edge2[1:iter,], 2, mean)[1]),
sd = c(apply(mat_rand[1:iter,], 2, sd)[1],
apply(mat_con[1:iter,], 2, sd)[1],
apply(mat_edge1[1:iter,], 2, sd)[1],
apply(mat_edge2[1:iter,], 2, sd)[1])
)

Estimated_b <-
data.frame(
mean = c(apply(mat_rand[1:iter,], 2, mean)[2],
apply(mat_con[1:iter,], 2, mean)[2],
apply(mat_edge1[1:iter,], 2, mean)[2],
apply(mat_edge2[1:iter,], 2, mean)[2]),
sd = c(apply(mat_rand[1:iter,], 2, sd)[2],
apply(mat_con[1:iter,], 2, sd)[2],
apply(mat_edge1[1:iter,], 2, sd)[2],
apply(mat_edge2[1:iter,], 2, sd)[2])
)

Estimated_c <-
data.frame(
mean = c(apply(mat_rand[1:iter,], 2, mean)[3],
apply(mat_con[1:iter,], 2, mean)[3],
apply(mat_edge1[1:iter,], 2, mean)[3],
apply(mat_edge2[1:iter,], 2, mean)[3]),
sd = c(apply(mat_rand[1:iter,], 2, sd)[3],
apply(mat_con[1:iter,], 2, sd)[3],
apply(mat_edge1[1:iter,], 2, sd)[3],
apply(mat_edge2[1:iter,], 2, sd)[3])
)

###図示して確認###
#pdf("Sampling_Estimation.pdf")
plot(Estimated_a[, 1], ylim=c(1.5, 2.5), ylab="Value_a", pch=16,
col=c("black", "red", "blue", "tomato")
)
segments(0, 2, 5, 2, lwd=2, col="green")
legend(3.3, 2.5, c("Rand", "Con", "Edge1", "Edge2"), pch=16, col=c("black", "red", "blue", "tomato"))
legend(2.5, 2.5, c("True"), lwd=2, col="green")

for(i in 1 : 4) {
arrows(i, Estimated_a[i, 1],
i, Estimated_a[i, 1] + Estimated_a[i, 2], angle=90)
arrows(i, Estimated_a[i, 1],
i, Estimated_a[i, 1] - Estimated_a[i, 2], angle=90)

}

plot(Estimated_b[, 1], ylim=c(-8, 8), ylab="Value_b", pch=16,
col=c("black", "red", "blue", "tomato")
)
segments(0, 0, 5, 0, lwd=2, col="green")
legend(3.3, 7, c("Rand", "Con", "Edge1", "Edge2"), pch=16, col=c("black", "red", "blue", "tomato"))
legend(2.5, 7, c("True"), lwd=2, col="green")

for(i in 1 : 4) {
arrows(i, Estimated_b[i, 1],
i, Estimated_b[i, 1] + Estimated_b[i, 2], angle=90)
arrows(i, Estimated_b[i, 1],
i, Estimated_b[i, 1] - Estimated_b[i, 2], angle=90)

}

plot(Estimated_c[, 1], ylim=c(35,45), ylab="Value_c", pch=16,
col=c("black", "red", "blue", "tomato")
)
segments(0, 40, 5, 40, lwd=2, col="green")
legend(3.3, 45, c("Rand", "Con", "Edge1", "Edge2"), pch=16, col=c("black", "red", "blue", "tomato"))
legend(2.5, 45, c("True"), lwd=2, col="green")

for(i in 1 : 4) {
arrows(i, Estimated_c[i, 1],
i, Estimated_c[i, 1] + Estimated_c[i, 2], angle=90)
arrows(i, Estimated_c[i, 1],
i, Estimated_c[i, 1] - Estimated_c[i, 2], angle=90)

}
#dev.off()

描けた図
Randがランダムサンプリング。それ以外はバイアスサンプリングのパターン。
Sampling_Estimation.pdf
posted by しばきん at 11:38| 神奈川 ☁| Comment(0) | TrackBack(0) | 統計 | このブログの読者になる | 更新情報をチェックする

2010年10月25日

GLMMの切片

最近Rでglmmを使った解析がお手軽にできるようになりました。でも、僕が解析したときは、得られたパラメータを使って、実際のデータにプロットするとおかしなことになりました。対数正規分布の期待値はその分散にも影響される、ということを知らなかったためです。その事実を知らなくて、随分苦労しました。ちなみに北大の久保先生のところに解説があります。
http://hosho.ees.hokudai.ac.jp/~kubo/memo/seibutukagaku/kubostat.pdf

というわけで自分でも実験してみます。

###################
###glmm切片補正###
###################
set.seed(1)
library(glmmML)
###パラメタの設定###
a <- 1
b <- 3
sigma <- 1.9
N <- 1000

###共変量###
X <- sample(seq(0, 1, length=N), N, rep=T)

###sigmaの分だけ、ばらつかせてやる###
###これで、sigmaを考慮しないモデルでは過分散になる###
Y <- rpois(N, lambda=exp(a + b*X + rnorm(N, 0, sigma)))
id <- seq(1, N)

###glmによる結果###
res <- glm(Y ~ X, family="poisson")
summary(res)

###glmmの結果###
res2 <- glmmML(Y ~ X, cluster=id, family=poisson)
summary(res2)

###sigma^2/2を足してやる。予測値の切片にプラス###
###切片と正規誤差に、傾きで説明できない分が分解されたと考える###
###正規誤差は対数正規分布に従うので、Yを予測するときには###
###正規誤差の寄与分を考慮する必要がある###
plot(X, Y, ylim=c(0, 500))

pred_X <- seq(0, 1, length=N)

###正解のプロット###
points(pred_X,
exp(a + b*pred_X + sigma^2/2), col="green",
type="l", lwd=3)

func_m <- function(X) {
coef(res)[1] + coef(res)[2] * X
}

###GLMによるプロット###
points(pred_X, exp(func_m(X=pred_X)), main="GLM", col="red",
type="l", lwd=3)

###GLMMによるプロット###
###変量効果の分散も期待値に反映###
func_mm <- function(X) {
exp((coef(res2)[1] + coef(res2)[2]*X +
res2$sigma^2/2))
}

###正解値のものとほぼ同一の直線に乗る###
points(pred_X, func_mm(X=pred_X), col="blue",
type="l", lwd=3)

legend(x=0, y=500, c("True", "GLM", "GLMM"),
col=c("green", "red", "blue"), lwd=3)

描けた図glmm.pdf

###################
###glmm切片補正###
###################


やはり分布形に関する基礎知識を固めねば、と思った次第です。
ラベル:GLMM
posted by しばきん at 12:59| 神奈川 ☁| Comment(0) | TrackBack(0) | 統計 | このブログの読者になる | 更新情報をチェックする

2010年10月13日

小学生レベルとは思えない

小学生がまとめたアマガエルの体色変化に関する研究
http://www.higo.ed.jp/edu-c/kagakuten/h19list/pdf/kou11.pdf

すごいすごすぎる。僕が小学生のころにお父さんに手伝ってもらって、光ファイバーで火山噴火の様子を作ったり、川底をさらって水中生物の特徴を調べたのより、ずっとすごい。この研究自力でやったんだとしたら末恐ろしい。
posted by しばきん at 11:05| 神奈川 ☁| Comment(2) | TrackBack(0) | 日記 | このブログの読者になる | 更新情報をチェックする

2010年10月12日

Rで並列計算


#################
###Rで並列計算###
#################

###並列計算用に関数作成###
start_simulation <- function() {

###データの生成###
a <- rnorm(400)
b <- rnorm(400)

###回帰###
c <- lm(a ~ b)

###とりあえず、係数でも記録します###

return(c(a = coef(c)[1],
b = coef(c)[2])
)
}

###繰り返し数###
k <- 10

library(snow)
###乱数の都合でこいつも使います。ないと面白いことになります。###
library(rlecuyer)


#Core 2 Duoなので2つです。クアッドコアとかなら4ですかね。開きすぎると閉じさせられます。
#Firewallは無視です。ブロック解除です。
cl <- makeCluster(2, type = "SOCK")

###乱数種はここでセット###
clusterSetupRNG (cl,seed=123)

###並列計算するfunctionをここで明示###
clusterExport(cl, "start_simulation")

###スタート。ついでに時間も計っちゃいましょう###
###係数はpに、時間はtimesに入るようになっています###
times <- system.time(p <- parSapply(cl, 1:k, function(x){start_simulation()}))

###applyするときは転置してからであります###
apply(t(p), 2, mean)
apply(t(p), 2, sd)

#################
###Rで並列計算###
#################


しかし、あんまりにも複雑で長い時間がかかるものを並列計算させると、終わってるはずの計算がいつまで経っても終わらないこともある(経験談)ので、少し注意が必要みたいです。

posted by しばきん at 17:20| 神奈川 ☁| Comment(0) | TrackBack(0) | R | このブログの読者になる | 更新情報をチェックする
×

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