スキップしてメイン コンテンツに移動

!元ネタ

みなさん、こんばんは。「アタシだけのヒト」みつかりましたか?[[天丼です|http://ja.wikipedia.org/wiki/%E5%A4%A9%E4%B8%BC#.E3.81.8A.E7.AC.91.E3.81.84.E3.81.AE.E6.A5.AD.E7.95.8C.E7.94.A8.E8.AA.9E.E3.81.A8.E3.81.97.E3.81.A6.E3.81.AE.E5.A4.A9.E4.B8.BC]]。完全に「定番の挨拶」として流行らせようとしてますね。

いや、なかなかみつからないものですね、僕なんかはもうそろそろ探索範囲を広げてみようと考えていますよ。2次元の世界あたりはどうでしょうか? 3次元より狭いので探索空間はぐっと狭くなります、これは良いアイディアです。次元を落すのは高次元のデータ解析の基本ですからね。

さてそもそも「アタシだけのヒト」とはどんなヒト? きっと自分に対して誠実なひとなのではないか? と私は考えました。誠実が服を脱いで歩いていると巷で評判の私の言うことなので、どうかと思いますが、今日はそういうことにして、この問題についてベイズ的に考えてみます。

誠実なヒトは、約束を守るヒトです。アナタからみてアノヒトは
*約束を守る性格
*約束を破る性格
のどちからを取ります。また、アナタは
*約束を守られる
*約束を破られる
のどちらかを体験します。

このとき、アノヒトが「約束を守るひとである、と確信していく 」= 「アタシだけのヒトであると確信していく」様子をベイズの定理とベイズ更新 (bayes updating) を使ってシミュレーションしています。その筋では reputation とか呼ばれている話です。

まず誠実なヒトを定義しましょう。
*誠実: 約束を10回のうち8回まで守る
*不実: 約束を10回のうち3回しか守らない
誠実なヒトでも失敗もありますので8/10ということにします。マザーテレサも「愛したいなら許すことを知れ」的なことを言っていますしね。逆に不実のひととは3/10とします、嫌なやつです。

アタシだけのヒトも定義します。
*確信が0.9以上になったら「アタシだけのヒト」
*確信が0.1以下になったら「アタシだけのヒトじゃない」
このように決断することにします。肝心の確信 (belief) の度合いですが、これは
*P(性格)が事前確率
*P(約束|性格)が尤度関数
としたときの事後確率P(性格|約束)が確信度になります。t回目の約束のときにはt-1の事後確率を事前確率として使うことで確率を更新していきます。これがベイズ更新です。tDiaryで数式書く設定してないので、コードみてください(あとで追記するかも)

さて、アノヒトと100回ほど約束をします。そのたびに「アタシだけのヒト」であるという「確信」がどう移り変わるでしょうか?

ここでは、仮にアノヒトは7/10の割合で約束を守る人だとしましょう。でもあなたはこれを知りえないので、実際に約束をしてみて、経験から確信度を変化させてみるわけですね。

もし最初の約束をしたときに、「アタシだけのヒト確信度」が五分五分だと思っているとします。以下のグラフの黒い線が確信の変化です。平均して20回ぐらい約束を繰り返してみて、「アタシだけのヒト」と確信しています。この線は、何度も人生をやりなおしみて、その平均を見ています。1000回生れ変わった設定です。灰色の線は1回ごとのシナリオを示しています。セーブしておいてそこから何度も始めてみた、という感じです。bad end になったシナリオがいくつかありますね。


{{image 0, 'シナリオ1', nil, [480,480]}}



確信度が2/10 から始めてみます。32回目にしてやっと確信していますね。また dead end なシナリオも増えています。第一印象重要!!! 第一印象最悪でも最後はくっつく、みたいなおいしい話はアニメでは良く聞くのに、なぜか僕がそのようなめに会わないのもわかります。でも本当に良いひとだったらねばる価値はありそうですね。また確信が形成されるまでの間に、灰色の線が暴れています。この期間のみかけ確信度が高いところで決断してしまうと、判断すると見誤ることが可能性が出てきますね。


{{image 1, 'シナリオ2', nil, [480,480]}}



4.5/10回の嘘つきに、最初の確信度0.8で一目惚れです。4.5ってところがアメとムチな感じでやらしいです。

平均して50回の約束ぐらいでゆっくりと気付くのですが、個別のシナリオ、つまり灰色の線をみると暴れまくりです。速攻バレてます。まあ、ずるずると付き合ってしまう、みたいなシナリオもないことないようですがレアのようです。これもお話の世界だけなのでしょうかね。ちなみに5.5回の嘘つきさんだと100回約束の間には「アタシだけのヒトじゃない」という決断には至りませんでした、やらしいですね。



{{image 2, 'シナリオ3', nil, [480,480]}}



全体を通して面白いのは一度確信に至ると決断が変わらないことが多いようです。この現象は、ウォールステッターさんという国際政治学者が cry-wolf syndrome と呼び、真珠湾攻撃、キューバ危機のような政治的誤算が生じる理由として議論したそうです。

まあ、わーわー言いましたが、そもそも相手が自分のことを「アタシだけのヒト」だと確信してくれるかどうかは、また別の話ですけどね! 仮に誠実と思われても「良い人」どまりってこともあるし。ままならね〜

以下、Rのコードを載せておきます (つづく...かもしれない)。こなれてないコードですが...
<<<
rm(list=ls())

wolfboy <- function (p, simNum, w0, true_prob) {

w <- rep(0, simNum)
z <- rep(0, simNum)
w[1] <- w0

for ( i in 2:simNum ) {
z[i] <- ifelse( runif(1) > true_prob, 1, 0 )
#cat(i, z, "\n")
if (z[i] == 0) {
w[i] <- (w[i-1]*p[1,1]) / ( w[i-1]*p[1,1] + (1-w[i-1])*p[2,1] )
} else {
w[i] <- (w[i-1]*p[1,2]) / ( w[i-1]*p[1,2] + (1-w[i-1])*p[2,2] )
}
}
return(w)
}

simNum <- 100

# def honest man
p <- matrix(c(8/10, 2/10, 3/10, 7/10), nrow=2, byrow=T)
threshold.heaven <- 0.9
threshold.hell <- 0.1

# your mind
w0 <- 0.5
#w0 <- 0.2
#w0 <- 0.8

# man
true_prob <- 0.7
#true_prob <- 0.45

# sim.
ws <- replicate( 1000, wolfboy(p, simNum, w0, true_prob) )
w.mean <- rowMeans(ws)
pass.index <- max( which(w.mean <= threshold.heaven) )
cat(pass.index, "\n")

# plot
png("wolf.png")
matplot(1:100, ws[,1:30], type="l", col="gray", ylim=c(0,1),
xlab="num. of date", ylab="Prob. of honest")
lines(w.mean, type="l", col="black")
abline(h=threshold.heaven, col="red")
abline(h=threshold.hell, col="blue")
abline(v=pass.index, col="red")
dev.off()
>>>

以下の2冊に元ネタの「羊飼いとオオカミ」(オオカミがきたぞーって嘘をつく少年のアレ)の例が載っています。これを参考にしました。1のほうがやわらかい日本語で書かれていて読みやすいと思います。

#[[入門ベイズ統計—意思決定の理論と発展|http://www.amazon.co.jp/gp/product/4489020368?ie=UTF8&tag=itoshitv-22&linkCode=as2&camp=247&creative=7399&creativeASIN=4489020368]]
#[[意思決定の基礎 (シリーズ意思決定の科学)|http://www.amazon.co.jp/gp/product/4254295111?ie=UTF8&tag=itoshitv-22&linkCode=as2&camp=247&creative=7399&creativeASIN=4254295111]]

参考文献: {{amazon "4063343839"}}
きっと、ちぃにはベイズ意思決定回路が実装されるに違いない。

このブログの人気の投稿

シーケンスアダプタ配列除去ツールまとめ

FASTQ/A file からシーケンスアダプター配列やプライマー配列を除くためのプログラムをまとめてみる。

まず、配列の除去には大別して2つの方向性がある。ひとつは、アダプター配列を含む「リード」を除いてしまう方法。もうひとつは除きたい配列をリードからトリムする方法である。後者のほうが有効リードが増えるメリットが、綺麗に除ききれない場合は、ゲノムへのマップ率が下がる。
気をつける点としては、アダプター/プライマーの reverse complement を検索するかどうか。paired end の際には大事になる。クオリティでトリムできるものや、Paired-end を考慮するものなどもある。アダプター/プライマー配列の文字列を引数として直接入力するものと、multi fasta 形式で指定できるももある。

From Evernote: シーケンスアダプタ配列除去ツールまとめTagDust
http://genome.gsc.riken.jp/osc/english/software/src/nexalign-1.3.5.tgz http://bioinformatics.oxfordjournals.org/content/25/21/2839.full
インストール: curl -O http://genome.gsc.riken.jp/osc/english/software/src/tagdust.tgztar zxvf tagdust.tgz cd tagdust/ make sudo make install rehash
使いかた: tagdust adapter.fasta input.fastq -fdr 0.05 -o output.clean.fastq -a output.artifactual.fastq
解説: 入出力形式は fastq/a が使える。リード全体を除く。速い。アダプター配列を fasta 形式で入力できるのが地味に便利で、これに対応しているものがなかなかない。Muth–Manber algorithm (Approximate multiple string search) を利用。FDRを指定できる。GPL3

Quartz-Seqで1細胞/微量RNA-Seqを始めたい方へ

はじめに 新しい高精度な1細胞RNA-Seq, Quartz-Seq論文を出してから、各方面から多く相談を受けています。
Sasagawa Y and Nikaido I, et. al. Quartz-Seq: a highly reproducible and sensitive single-cell RNA-Seq reveals non-genetic gene expression heterogeneity. Genome Biology. 14. 2013 
そこで、新しく1細胞RNA-Seqを始める方へ、僕達が理想だと考えている技術導入の手順を紹介したいと思います。また我々の方法は1細胞(6-14 pg Total RNA)だけでなく pg-ng オーダーの少量RNAからシーケンスが可能です。そのような方も以下の手順が参考になると思います。 0. 1細胞/微量RNA-Seqが本当に必要なのか検討する 1細胞/微量RNA-Seqでは、現時点でQuartz-Seqが世界最高の性能を持っている訳ですが、十分なサンプルを用意し、通常のRNA-Seqしたほうが、より精度の高いデータが得られます。なので、基本的には、サンプルをたくさん集める方法をしっかり検討すべきです。まずは、戦略面と技術面で1細胞/微量RNA-Seqが本当に必要かを検討する基準について書きます。 0.1. 戦略面での検討 あなたが抱えているプロジェクトが、1細胞/微量RNA-Seqでなければアプローチできないかどうかを問い直すことが重要です。
基本的には以下の2つの状況で、1細胞/微量RNA-Seqが役に立ちます。
a. 細胞状態が連続的に変化し、さまざまな細胞状態が、細胞集団に含まれている場合 (振動現象、ゆらぎなど) b. 細胞状態を特定するマーカーがほどんどわかっていない場合
最初から細胞状態が2状態しかないことが明らかで、しかも細胞状態を代表する遺伝子が分かっている、という状況では、FACSなどで cell sorting し、目的の細胞を採取することを考えるべきです。そして、微量RNA-Seqや通常のRNA-Seqで、しっかりと biological replication を取る方が良いでしょう。微量になると、テクニカルなノイズが増えるために、生物学的な差を知るためには、n を…