Hatena::ブログ(Diary)

驚異のアニヲタ社会復帰への道

Prima Project

2017-07-11

MikuHatsune2017-07-11

ガンマ分布:癌にかかるまでの時間

ガンマ分布は、期間 ¥mu ごとに1回くらい起こるランダムな事象N 回起こるまでの時間の分布としてモデル化されることがあるらしい。

ガンマ分布自体は、shape とscale (もしくは逆数のrate) のふたつのパラメータで決定される。

 

悪性腫瘍(いっぱんにがん)は、遺伝子に変異が蓄積してがんとして発病する、と言われる。ここで、ひとりのヒトに変異が1年あたり¥mu 回入るとき、これがN 回蓄積すると発癌する、とする。ここで、変異の修復は考えない。

1年あたりに1回以下の珍しい頻度なので、変異が入る事自体はポアソン分布から得られる。

ベタにシミュレーションしつつ、ある年で生じた変異を翌年以降も引き継ぎつつ、すべての変異発生数を推移行列で格納しておく。

ここで、ガンマ分布でさくっとやってしまうと、shape が発癌に必要なヒット回数、scale(rate) が1年あたりのヒット回数¥mu になる。

f:id:MikuHatsune:20170711203834p:image

# simulation
n.ind <- 10000 # 人数
n.hit <- 10    # 発癌に必要なヒット回数
age.box <- 1:100
m2 <- 0.2      # 1年あたりのヒット回数。ポアソン分布の平均パラメータ
first.cancer <- rep(-9,n.ind)
for(i in 1:n.ind){
  events <- rpois(length(age.box),m2)
  cumsum.events <- cumsum(events)
  tmp <- which(cumsum.events >= n.hit)
  if(length(tmp)>0){
    first.cancer[i] <- tmp[1] - 1
  }
}

# poisson model
d <- dpois(0:(n.hit-1), m2) # less than n.hit probability
d <- c(d, 1-sum(d))         # add more than n.hit probability

# transition matrix
# All transition patterns are considered.
idx <- mapply(function(z) seq(z, n.hit+z, by=1), 0:n.hit)
idx <- replace(idx, idx>n.hit, n.hit)

dage <- 0:100
res <- matrix(0, 0:n.hit, nr=length(dage), nc=n.hit+1)
res[1,] <- d
for(i in 1:(length(dage)-1)){
  trans <- outer(res[i,], d)
  #trans[lower.tri(trans)] <- 0
  res[i+1,] <- tapply(c(trans), c(idx), sum)
}

par(mar=c(5, 5, 2.2, 0.5), cex.axis=1.5, cex.lab=1.5, font=2)
hist(first.cancer,breaks=seq(from=-10,to=100,by=10), freq=FALSE, col=2, density=10, xlab="first cancer age")
lines(head(dage, -1), diff(res[, n.hit+1]), col=4, lty=3, lwd=3)
dg <- dgamma(dage, n.hit, rate=m2)
lines(dage, dg, col=3, lwd=3)
legend("topleft", legend=c("simulation", "pois model", "gamma model"), col=c(2,4,3), pch=15)

 

さて、ここで実際の癌の発生状況を確認してみる。国立がん研究センターから集計表ダウンロードして適当にネ申エクセルを使う部分だけtxt にする。

f:id:MikuHatsune:20170711204133p:image

f:id:MikuHatsune:20170711204210p:image

d <- read.delim("cancer.txt", check.names=FALSE, stringsAsFactors=FALSE)
ye <- unique(d$診断年)
idx <- grep("^C\\d+", d)

m <- mapply(function(z) d[(idx[z]-2):(idx[z+1]-3)], head(seq(idx), -1))

ageidx <- grep("歳", colnames(d))
agename <- colnames(d)[ageidx]
age <- sapply(sapply(str_extract_all(agename, "\\d+"), as.numeric), mean)
dat <- mapply(function(z) split(z, z$性別), split(d, d$コード))
dat[[14]] <- list(dat[[14]][[1]], dat[[15]][[1]], dat[[16]][[1]])
dat[[13]] <- list(dat[[13]][[1]], dat[[17]][[1]], dat[[18]][[1]])
dat[15:16] <- dat[17:18] <- NULL

yl <- range(0, d[,ageidx]/rowSums(d[,ageidx]))
# 確率密度
cols <- jet.colors(nrow(tmp))
i5 <-ye %% 5 == 0
par(mfrow=c(length(dat), 3), mar=c(2.5, 3, 2.2, 0.5), cex.axis=1.5, font=2, las=1)
for(i in seq(dat)){
  for(j in c(1,3,2)){
  tmp <- dat[[i]][[j]]
  matplot(age, t(tmp[,ageidx]/rowSums(tmp[,ageidx])), type="l", lty=1, col=cols, ylab="Density", ylim=yl, lwd=3)
  title(sprintf("%s: %s", tmp[1,2], tmp[1,4]), cex.main=2)
  legend("topleft", legend=ye[i5], col=cols[i5], pch=15, ncol=2, cex=1.2)
  }
}

# 発生数
par(mfrow=c(length(dat), 3), mar=c(2.5, 3, 2.2, 0.5), cex.axis=1.5, font=2)
for(i in seq(dat)){
#par(mfrow=c(1, 3))
  yl <- range(0, do.call(rbind, mapply(function(z) z[,ageidx], dat[[i]], SIMPLIFY=FALSE)))
  for(j in c(1,3,2)){
  tmp <- dat[[i]][[j]]
  matplot(age, t(tmp[,ageidx]), type="l", lty=1, col=cols, ylab="Density", ylim=yl, lwd=3)
  title(sprintf("%s: %s", tmp[1,2], tmp[1,4]), cex.main=2)
  legend("topleft", legend=ye[i5], col=cols[i5], pch=15, ncol=2, cex=1.2)
  }
}

 

さてここで、癌発生年齢の経験分布から、ガンマ分布のパラメータを適当に当てはめて推定してみる。ガンマ分布のパラメータは発癌までの回数と、1年あたりのヒット回数なので、これを比較すると癌になるまでに要する変異数的なものがわかる。

実際には、変異はDNA損傷修復機構で修復されるし、実際に発癌に寄与する変異(driver mutation) は固有の遺伝子クリティカルな変異数個であることが多い(大腸がんや脳腫瘍など)。

若年での発症が極端には多くない大腸がんだけ取り出してみる。というのも、例えば白血病や脳・中枢神経系は10歳以下くらいの小児に後発するタイプの悪性腫瘍があるため、単純なガンマ分布だけでのフィッテイングだと当てはまりがよくない。また、若年発症するものは多くが遺伝が強く関係するタイプ(乳癌子宮体癌のBRCA1/2 など)(ならば大腸がんだってFAP があるのでは、と思ったけど統計上は若年層はあまりいないようである)。皮膚がんや肝臓がんあたりは遺伝というより紫外線肝炎ウイルスの寄与が大きいのでこれらはよりよくガンマ分布に当てはまりそうだと思ったけどそうでもなかった。

f:id:MikuHatsune:20170711205016p:image

推定の初期値によるのかもしれないが、ヒット回数は12回から25回くらいが多いようである。ヒット速度は0.2 から0.3 くらい。

f:id:MikuHatsune:20170711205337p:image

par(mfrow=c(1, 3), mar=c(2.5, 3, 2.2, 0.5), cex.axis=1.5, font=2, las=1)
#for(i in seq(dat)){
  for(j in c(1,3,2)){
  tmp <- dat[[i]][[j]]
  y <- t(tmp[,ageidx]/rowSums(tmp[,ageidx]))
  s1 <- mapply(function(z){
    ytmp <- y[,z]
    resid <- function(par){
      yhat <- dgamma(age, par[1], par[2])
      sum((ytmp-yhat)^2)
    }
    optim(init, resid, method="BFGS")$par}, seq(ncol(y)))
  dg <- mapply(function(z) dgamma(age, s1[1, z], s1[2, z]), seq(ye))
  dg <- sweep(dg, 2, colSums(dg), "/")
  yl <- range(0, do.call(rbind, mapply(function(z) z[,ageidx]/rowSums(z[,ageidx]), dat[[i]], SIMPLIFY=FALSE)))
  tmp <- dat[[i]][[j]]
  hoge <- t(tmp[,ageidx]/rowSums(tmp[,ageidx]))
  matplot(age, dg, type="l", lty=1, col=cols, ylab="Density", ylim=yl, lwd=3)
  title(sprintf("%s: %s", tmp[1,2], tmp[1,4]), cex.main=2)
  legend("topleft", legend=ye[i5], col=cols[i5], pch=15, ncol=2, cex=1.2)
  }
}

par(mar=c(5, 5, 2.2, 0.5), cex.axis=1.5, cex.lab=1.5, font=2)
plot(t(s1), pch=15, col=cols, xlab="発癌までのヒット回数", ylab="1年あたりのヒット回数")
legend("bottomright", legend=ye[i5], col=cols[i5], pch=15, ncol=2, cex=1.2)
title(sprintf("%s: %s", tmp[1,2], tmp[1,4]), cex.main=2)

2017-07-10

MikuHatsune2017-07-10

賭ケグルイの投票じゃんけんで蛇喰さんはどうやって芽亜里の奴隷であるクラスメイトの割合を推定したか

賭ケグルイを見ている。

黒髪ロングのプリーツスカート黒タイツのはやみんなので視聴意欲がやばい。

 

ここで、賭ケグルイの1話で、投票じゃんけんという変則じゃんけんで勝負する。

クラスN=30 人がグー G、チョキ C、パー P のいずれかの手を投票する。

蛇喰さんと芽亜里はその投票箱から重複なしで3枚カードを引く。

各々の持ち手はその3枚となり、あとはじゃんけんのルールに従って勝敗が決まる。出した手は捨てられる。

 

ということである。ここで、両者120万円(120枚のチップ)から始めて、出した手の順に(他はアニメ演出上、蛇喰さんが把握できた手について)

1回戦:2枚賭け 蛇喰勝(G??) 芽亜里(C??) アニメ絵では残りの2枚はG であることが分かっているが蛇喰さんには見えない

2回戦:50枚賭け 蛇喰負(GCC) 芽亜里(GG?)

3回戦:2枚賭け 蛇喰勝(P??) 芽亜里(G??) アニメ絵では残りの2枚はC と(たぶん) G であることが分かっているが蛇喰さんには見えない。

4回戦:50枚賭け 蛇喰負(PG?) 芽亜里(PP?)

その後、じゃんけんぽん、が3回は行われて、蛇喰さんがすべてのチップを奪われるが、即金で

N回戦:1000万賭け 蛇喰勝(C??) 芽亜里(PPG)

というふうになって最終的に蛇喰さんが勝利する。

 

ここでこの投票じゃんけんイカサマのネタバレだが、芽亜里はクラスメイト30人中21人の投票を操り、主人公の鈴井くんの提示した手(N回戦ではチョキC) 以外の手(つまりグーG とパーP)のどちらかに投票する。こうすると、クラスのp の割合が投票を操作されているとき、GCP の投票確率は

(G,C,P)=(¥frac{2+p}{6}, ¥frac{1-p}{3},¥frac{2+p}{6})となる。p=0 ですべて¥frac{1}{3}p=1 のときグーとパーのみ出るようになって各々¥frac{1}{2} である。

このとき、C が出る確率はかなり低いから、P を出しておけば少なくとも負けることはないだろう、というのが芽亜里の作戦である、たぶん。

 

ここで蛇喰さんは、2枚賭けのときは芽亜里が適当に手を出しているのに50枚賭けという大勝負のときには策を講じた手を出していること、自分の真後ろにたっていた主人公の手を手鏡で観察していたことから、N回戦のときに芽亜里にイカサマの存在を伝え、なおかつそれが、「10人から20人が芽亜里の指示通りに投票する」ということを見抜いた。

実際にこれを見抜けるかやってみる。

 

いま、GCP がある確率で出るとき(等確率 ¥frac{1}{3} とする)、自分の手札の揃え方は10通りある。この10通りが出る確率は多項分布であり、以下の通りである。

x <- expand.grid(0:3, 0:3, 0:3)
x <- x[rowSums(x) == 3,]
colnames(x) <- c("G", "C", "P")
prob <- c(1, 1, 1)/3
dm <- apply(x, 1, dmultinom, prob=prob)
cbind(x, dm)
   G C P         dm
4  3 0 0 0.03703704
7  2 1 0 0.11111111
10 1 2 0 0.11111111
13 0 3 0 0.03703704
19 2 0 1 0.11111111
22 1 1 1 0.22222222
25 0 2 1 0.11111111
34 1 0 2 0.11111111
37 0 1 2 0.11111111
49 0 0 3 0.03703704

ここで、蛇喰さんがイカサマを見抜けたのは、「(G,C,P)の手札が揃う勝負が極端に少ない」からという仮説にしておく。というのも、投票が偏らなければ投票箱は各々¥frac{1}{3} になるはずだから、とする。たぶん。

ここで、上の各手の出方のうち、(G,C,P)以外に絞れば、余事象的に「手札がすべて同じ(3枚)」と「2枚と1枚(と0枚)」は計算できる。

そしてこれは(G,C,P)か、それ以外にすれば、二項分布で計算できる。

二項分布から95%信頼区間を計算してもよいが、蛇喰さんは「10人から、万全を期すなら20人」といっているので、適当にrstan を使って事後分布で評価してみる。

 

例えば、1回の勝負で(G,C,P) が来たときの、投票箱の手の分布はこのような感じになる。

ここで、(G,C,P) の出る確率は足して1 なので、G+C+P=1 の3次元平面になるが、平面なので2次元平面に写像した図にしている。三角形頂点のGCP の点に近いほど、その手の出る確率が1 に近くなる。

f:id:MikuHatsune:20170710120034p:image

1回の勝負で2枚と1枚、例えば(G,G,C) のときはこうなる。

f:id:MikuHatsune:20170710120028p:image

1回の勝負で2枚と1枚、例えば(G,G,G) のときはこうなる。

f:id:MikuHatsune:20170710120025p:image

library(rstan)
code <- "
data{
  int trial;
  int gcp[3, trial];
}
parameters{
  simplex[3] a;
}
model{
  for(i in 1:trial)
    gcp[, i] ~ multinomial(a);
}

"
m <- stan_model(model_code=code)
trial <- 3
standata <- list(trial=trial, gcp=replicate(trial, c(1, 1, 1)))
fit <- sampling(m, standata, iter=30000)


M <- cbind(c(0,1,0.5), c(0,0,sqrt(3)/2)) # 2次元に写像する行列
triangle <- diag(1, 3) %*% M
xy <- ex$a %*% M
kd <- kde2d(xy[,1], xy[,2], c(bandwidth.nrd(xy[,1]), bandwidth.nrd(xy[,2])), n=500)

cols <- matlab::jet.colors(100)
i <- which(kd$z == max(kd$z), arr.ind=TRUE)
out <- c(kd$x[ i[1] ], kd$y[ i[2] ], 1) %*% solve(cbind(M, 1))
par(mar=c(3, 3, 2, 2), cex.axis=1.5)
image(kd$z, col=cols)
polygon(rbind(triangle[c(1:3,1),],c(-5,0),c(-5,5),c(5,5),c(5,-5),c(-5,-5),c(-5,0),triangle[1,]), col="white", border=NA)
polygon(triangle)
points(triangle, pch=16, cex=3, xpd=TRUE, col="hotpink")
text(triangle, c("G", "C", "P"), col="yellow", font=2, xpd=TRUE)
points(out %*% M, pch="★", col="hotpink", cex=1.5, xpd=TRUE)
legend("topright", legend=mapply(sprintf, "%s: %.4f", c("G", "C", "P"), out), bty="n", cex=2)

 

ここで大事なのが、確率的な問題なので、たとえG が出る確率が低くても(青の部分)、引く手が(G,G,G) になることはありえる、ということである。ましてや、等確率¥frac{1}{3} はそこそこありうる確率点である。

しかし、この分布も、勝負の回数が増えることで、等確率¥frac{1}{3} はおかしい、という状況が出てくる。

 

さてここで、買収されている(奴隷)の割合をp とすると、(G,C,P) が出る確率は¥frac{1-p}{3}、それ以外(2枚同じか、3枚とも同じ) 手が出る確率は¥frac{2+p}{3} となる。勝負が進んで蛇喰さんが手の内訳を確認できた、すなわち、蛇喰さん自身の手(これは勝負の回数分)と、芽亜里の手(初手で勝負が決まれば不明だが、2回目以降じゃんけんが進めばわかる)を足した回数、ということにする。つまり、アニメの4回戦が終わったときには、手を6回確認できていることとする。

芽亜里の奴隷の割合p が一様な事前分布を仮定したときの、N回手を確認したあとの事後分布は以下のようになる。

f:id:MikuHatsune:20170710220934p:image

f:id:MikuHatsune:20170710223536p:image

回が増えると分布は極端にp=1、つまりみんな奴隷に傾く。

 

¥alpha% 信用区間プロットした。12回くらい手のうちがわかると、10/30 人が95% CI、20/30 人が50% CI、つまり中央値になるので、ここらへんのことを言っているのだろうか。

f:id:MikuHatsune:20170710120754p:image

 

とはいっても、クラスメイト全員が芽亜里の奴隷であることはあるのはあるだろうが、ちょっと現実的ではない、かもしれない。このとき、奴隷割合p の事前分布をちょっといじってみる。

いま、一様分布を事前分布として仮定すると、p=0p=1 も等確率で生じてしまう。30人のクラスのうち、さすがに芽亜里が奴隷にできるのは半分もいかないだろう、という仮定をおけば、違う事前分布を設定できる。ここで、p は確率なので、適当にベータ分布からサンプリングしてみよう。ベータ分布のパラメータはふたつで決まる。ここで、ふたつのパラメータをどちらも1 にすれば、一様分布になる。

今回は適当に3と10を選んだ。おそらくp=0.2 あたりが最頻値になる。

f:id:MikuHatsune:20170710224623p:image

 

これで同じようにすると、50回手のうちを確認してもたかだかp=0.6 程度までしか事後確率は変動しない。

f:id:MikuHatsune:20170710225320p:image

 

¥alpha% CIプロットしてみると、10/30 人の奴隷、20/30 人の奴隷を推定値として採用するのは難しいような気がする。というのも、50枚賭けを2回負けた時点で、持ちチップが20枚くらいしかないので、1枚ずつ賭けて負けても50回手のうちを確認することができないからである。

とすれば、p=0.2 より多い割合で奴隷がいる、と事前分布を仮定するのが無難っぽい。というのも、転校してきて最初の勝負を自信満々にふっかけられたら何か裏があると思っておいたほうがいいっぽい。

f:id:MikuHatsune:20170710225821p:image

 

統計ガチ勢の人がいたらぜひとも数式的にごにょごにょ解いてください。

code <- "
data{
  int n;
}
parameters{
  real<lower=0, upper=1> p[n];
}
model{
  p ~ beta(3, 10);
  for(i in 1:n)
    i ~ binomial(i, (2+p[i])/3);
}
"

m <- stan_model(model_code=code)
standata <- list(n=50)
fit <- sampling(m, standata, iter=6000, warmup=1000)
ex <- extract(fit)
med <- colMeans(ex$a)
cint <- c(99.9, 99, 95, 90, 80, 75, 50)/100 # 信用区間
ci <- apply(ex$p, 2, quantile, (1-cint)/2)

# 累積確率分布
cm <- apply(ex$p, 2, sort)
par(mar=c(5, 5, 2, 2), cex.axis=1.5, cex.lab=1.5)
matplot(seq(0, 1, length=nrow(cm)), cm, type="l", col=matlab::jet.colors(ncol(cm)), lwd=3, lty=1, xlab="index", ylab="cumulative probability")
abline(0, 1, lty=3)

# 信用区間幅のプロット
meari <- 10:20
cols <- matlab::jet.colors(length(cint))
par(mar=c(5, 5, 2, 2), cex.axis=1.5, cex.lab=1.5)
matplot(t(ci), type="l", lwd=3, col=cols, lty=1, ylim=c(0, 1), xlab="カード内訳を確認した回数", ylab="クラス内の寝返り割合")
abline(h=meari/N, lty=3)
abline(h=n0/N)
legend("topleft", legend=sprintf("%.1f%s", cint*100, "%"), title="CI", lty=1, lwd=3, col=cols, bg="white")

# 分布のベタ書き
dxy <- mapply(function(z){
  dx <- density(unlist(z))
  i0 <- 0 < dx$x & dx$x < 1
  return(list(x=dx$x[i0], y=dx$y[i0]))
}, apply(ex$p, 2, list), SIMPLIFY=FALSE)

par(mar=c(5, 5, 2, 2), cex.axis=1.5, cex.lab=1.5)
yl <- range(0, dxy)
cols <- jet.colors(length(dxy))
plot(0, type="n", ylim=yl, xlim=c(0, 1), xlab="probability", ylab="density")
for(i in seq(dxy)){
  lines(dxy[[i]]$x, dxy[[i]]$y, col=cols[i], lwd=2)
}

2017-07-03

MikuHatsune2017-07-03

糖質制限ダイエットを始めたらたった1日で体重が2kg減った話をしたらもっと詳しい体重推移データをもらった

ダイエットのデータを昔使って遊んだけど、元データを持っている人が更に長期のデータを公開してくれた

さっそく遊んでみる。

 

2013年は測っていない時期が多すぎて、推定がクソ。

以前遊んだときは、ダイエットを意識して計画的に体重を測定していた2016年のものだったらしい。2015年とかは体重が増えている。

2016年ダイエット意識高い系()以降は油断して体重がまた増えている。

 

ダイエットさぼったとかそういう異常検知系のことができませんかね? というディスカッションになったけど頑張ればできるかも?

f:id:MikuHatsune:20170703210345p:image

2017-06-27

MikuHatsune2017-06-27

声優統計コーパスを使ってみる

声優統計コーパスというものがある。

日本声優統計学会

プロの女性声優 3 名が 3 パターンの感情表現で読み上げた音声 2 時間分 を「声優統計コーパス」として無料公開します - 糞ネット弁慶

 

音素バランス文という、音声言語研究では非常になんかいい例文があって、それをプロの声優に読み上げてもらうことで、テキストマイニング、音声研究に役立てようというデータベース

基本的に利用、解析、ダウンロードは無料で、「同人誌論文などで利用される場合」となぜか同人誌のほうが論文に先んじて書かれる始末。

声優は女性声優3人が、100の音素バランス文を普通に、喜んで、怒って、の3つの感情パターンで読み上げているため、音声の感情の研究にも使える。

BGM のない、アフレコ音声なので、音声合成などにも使えそうである。

ここで、音素バランス文の構築自体は、別の話なので上のリンクからたぶん関連記事に辿れるとおもう。

 

例えば21番の音素バランス文は、

 

痛みは点滴より鎮痛薬を静脈投与することで鎮痛を行う

イタミワテンテキヨリチンツーヤクヲジョーミャクトーヨスルコトデチンツーヲオコナウ

i,t/a,m/i,w/a,t/e,N,t/e,k/i,y/o,r/i,ch/i,N,ts/u:,y/a,k/u,o,j/o:,my/a,k/u,t/o:,y/o,s/u,r/u,k/o,t/o,d/e,ch/i,N,ts/u:,o,o,k/o,n/a,u

 

という文を、3人の女性声優が3感情パターンの9通りで読んでくれるので、これだけで一生遊んで暮らせるくらいのレベルである。

 

ひとまず、これらをダウンロードしてcorpus というディレクトリにいれたとして、3人の声優の声と、感情と、それらの9パターンの組み合わせがなんとなくもやっと分離できるかやってみる。やっていることはここの焼き直し。

 

21番目の音素バランス文の、藤東知夏の怒りをスペクトログラムにした。実はコレ、土谷麻貴の怒りと一緒のように聴こえるのだが??

f:id:MikuHatsune:20170627154201p:image

library(seewave)
library(tuneR)
library(phonTools)
library(stringr)

# ID, 声優名、感情を取り出しておく
l <- do.call(cbind, mapply(strsplit, list.files(), "_"))
f <- list.files(recursive=TRUE, full.names=TRUE)
l <- list(cv=unique(l[1,]), em=unique(l[2,]), id=unique(str_extract(f, "\\d+")))

wav <- mapply(readWave, f)
idx <- which(grepl("021", f))
w <- wav[[ idx[1] ]]
spectro(w, flim=c(0, 5)) # スペクトログラム

 

収録時間をプロットしてみた。土谷麻貴の収録時間はほかの二人に比べて短い気がする。といっても、後ろの空白時間がけっこうあって、これは空白時間を削除するプログラムの仕様なのかもしれない。

f:id:MikuHatsune:20170627154452p:image

sec <- mapply(function(z) length(z@left) / z@samp.rate, wav)
s0 <- mapply(function(z) sec[grep(paste(l[[1]][z], l[[2]][z], sep="_"), f)], seq(l), SIMPLIFY=FALSE)
par(mar=c(5, 5, 2, 2), las=1, cex.lab=1.5, cex.axis=1.5)
matplot(do.call(cbind, s0), pch=16, xlab="Voice ID", ylab="Recoding second")
legend("topleft", legend=l[[1]], pch=16, col=seq(l[[1]]))

 

さてここで、100文3声優3感情の900音声について一括してなんかやってみる。音声データから各音声の特徴を表しているであろうフォルマントをごっそり上から数個取ってきて、適当に3次元に次元を落としていい感じに分かれるのかどうかやってみる。

ここでフォルマントとはという話は一切しない。

 

適当にbandwidth をいじるとすべての音声について少なくとも6つフォルマントがとれたので、6次元データをtsne により3次元にして3D プロットした。なんかいい感じにクラスターはできそうである。

f:id:MikuHatsune:20170627154806p:image

 

声優ごとにわける。藤東知夏土谷麻貴上村彩子である。たいしていい感じにわかれていないようである。

f:id:MikuHatsune:20170627154803p:image

 

感情ごとにわける。がangry、happyがnormalである。たいしていい感じにわかれていないようである。

f:id:MikuHatsune:20170627154800p:image

 

声優と感情ごとにわける。藤東知夏angry、マゼンダ藤東知夏happyピンク藤東知夏normal、

土谷麻貴angry、オレンジ土谷麻貴happyマルーン土谷麻貴normal、

暗緑上村彩子angry、上村彩子happy黄緑上村彩子normal である。たいしていい感じにわかれていないようである。

f:id:MikuHatsune:20170627154757p:image

 

今回は旬が命でたいした解析はできなかったが、音声、感情、テキスト、となんでも使えそうなので使っていこう。

library(foreach)
library(doSNOW)
library(tsne)

# 少し時間がかかるので並列処理
cl <- makeCluster(8, type="SOCK")
registerDoSNOW(cl)

frm <- foreach(i = seq(wav), .combine=c, .packages="phonTools") %dopar% {
  findformants(wav[[i]]@left, fs=wav[[i]]@samp.rate, verify=FALSE, maxbw=1000)
}
stopCluster(cl)

res <- frm[seq(frm) %% 2 ==1]

# tsne で適当に3次元に落とす
h <- t(sapply(res, head, min(sapply(res, length))))
sne <- tsne(h, k=3)

cols <- c("red", "yellow", "darkseagreen")
# 声優ごと
idx <- apply(mapply(function(z) grepl(z, f), l[[1]]), 1, which)
plot3d(sne, col=cols[idx], type="s", radius=2)

# 感情ごと
idx <- apply(mapply(function(z) grepl(z, f), l[[2]]), 1, which)
plot3d(sne, col=cols[idx], type="s", radius=2)

# 声優と感情の9パターン
cols0 <- c("red", "magenta", "hotpink", "yellow", "orange", "maroon", "darkseagreen", "green", "lightgreen")
idx9 <- c(outer(l[[1]], l[[2]], paste, sep="_"))

idx <- apply(mapply(function(z) grepl(z, f), idx9), 1, which)
plot3d(sne, col=cols0[idx], type="s", radius=2)


 

2017-06-19

MikuHatsune2017-06-19

合格率60% の認定試験でほとんどの国民は一般人になれるか

虚構新聞のネタなので、一般人を試験により認定するという話は存在しません。

 

こんな記事を見かけた。

一般人認定試験、来年度実施を検討 「共謀罪」成立受け

共謀罪について、一般人であるかどうかの認定試験を行って、一般人か否かが決まる。ここで、「合格率は60%程度。繰り返し受検すればほとんどの人が一般人になれる」と言っている。

 

ここで、合格率が60% の試験を何回受けたらほとんどの人が合格して一般人になれるか考える。

10万人の受検者がいて、1回の試験の合格率が60% のとき、このシミュレーションでは13回目で全員が合格した。

n <- 100000
p <- 0.6
x <- 0:15
res <- c(n, rep(0, length(x)-1))
for(i in 1:(length(res)-1)){
  m <- rbinom(res[i], 1, p)
  tab <- table(factor(m, 0:1))
  res[i+1] <- tab[1] # 合格していない人
}

par(mar=c(4.5, 5, 2, 2), cex.lab=1.5, las=1)
plot(head(x+1, -1), tail(res/n, -1), xlab="初めて合格したときに受検していた回数", ylab="合格できていない受験者の割合")
 [1] 100000  40084  16161   6473   2627   1057    424    166     56     25      5      1      1      0
 [15]      0      0

f:id:MikuHatsune:20170619200320p:image

 

各施行がベルヌーイ施行のとき、1回の試験での0/1 binary アウトカムは二項分布から得られた。ここで、二項分布の拡張で負の二項分布というものがある。これは、r 回の成功を得るまでに必要な施行回数、とか、r 回の成功を刷る前に失敗した試行回数、とか言われる。

ここで、1回目にいきなり成功する確率と、1回目は失敗して2回目に成功する確率の和を考える。成功確率が0.6 のとき、r+(1-r)r で与えられる。

欲しい成功回数が1回のときは、負の二項分布は幾何分布 *geom と同じである。

p <- 0.6
p + (1-p)*p
[1] 0.84

ここで、今回の状況では、1回の成功(合格)までに何回受検するかだから、dnbinom では


dnbinom(x, 1, p)
 [1] 6.000000e-01 2.400000e-01 9.600000e-02 3.840000e-02 1.536000e-02 6.144000e-03 2.457600e-03
 [8] 9.830400e-04 3.932160e-04 1.572864e-04 6.291456e-05 2.516582e-05 1.006633e-05 4.026532e-06
[15] 1.610613e-06 6.442451e-07

pnbinom では

pnbinom(x, 1, p)
 [1] 0.6000000 0.8400000 0.9360000 0.9744000 0.9897600 0.9959040 0.9983616 0.9993446 0.9997379 0.9998951
[11] 0.9999581 0.9999832 0.9999933 0.9999973 0.9999989 0.9999996

となる。ここで、x は0 から始めているので、1回の成功までに失敗する回数が返ってくる。

 

理論的に求めると、10回受検しても10^{-4}、つまり1万人に1人は合格できないし、10万人に1人にするなら13回必要になる。

h <- pnbinom(x, 1, p)
e10 <- 10^(-(2:5))

par(mar=c(4.5, 5, 2, 2), cex.lab=1.5, las=1)
plot(x+1, log10(1-h), xlab="初めて合格したときに受検していた回数", ylab="合格できていない受験者の割合 [log10]")
abline(h=log10(e10), lty=3)

f:id:MikuHatsune:20170619200317p:image

 

シミュレーションと重ねてみると、受検者が減った10回目くらいまでの試験では理論値とよく一致する。

par(mar=c(4.5, 5, 2, 2), cex.lab=1.5, las=1)
plot(x+1, log10(1-h), xlab="初めて合格したときに受検していた回数", ylab="合格できていない受験者の割合 [log10]")
points(head(x+1, -1), tail(log10(res/n), -1), pch=15)
abline(h=log10(e10), lty=3)

f:id:MikuHatsune:20170619200311p:image