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

!元ネタ

みなさん、こんばんは。「アタシだけのヒト」みつかりましたか?[[天丼です|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

DNAを増幅するサーマルサイクラーを自作してみたよ

DNAをPCR法で増幅するために必要なサーマルサイクラーを自作してみました。自作と言っても、いわゆる、PCの自作と同じでパーツを組み立てていく感じです。購入から組み立ての様子を簡単に紹介します。 モチベーション ラボには様々なレクリエーションがあります。例えば、単にどこかに遊びに行ったり、スポーツ大会したり、ひたすら合宿形式でプログレスのプレゼンをするミーティングするなどがあります。それもよいのですが、せっかくなので、普段の研究時間ではトライできないが、研究に関わる hack を行う、というイベントを企画してみました。夏休みの自由研究や社会科見学的なノリです。   うちのラボでは、PCRを使ったウェットの実験技術の開発をしてきました。しかし、サーマルサイクラーのハードウェアの仕組みを体験的に理解している訳ではありません。そこで、サーマルサイクラーを作ってみました。   欧米で始まっている、自宅のガレージやキッチンでバイオロジーを行うムーブメント、バイオパンク、DIYbio を体験しておきたいというのもありますし、Arduino などオープンハードウェア、Maker のムーブメントを体験するのも目的の一つです。ハードウェア開発が思っているほどハードルが下っていることを体験できて、かつ、将来、ウェットの開発だけでなく、装置開発などもできたら、ラッキー、ぐらいの気持ちでやってみました。   購入 今回作ったのは、組み立て式で、かつ、仕様などや設計図が公開されているOpenPCRというサーマルサイクラーです。ハードウェアの仕様・設計図、制御ソフトウェアなどの情報がすべて公開されており、部品からも自作することが可能です。今回は、「設計図から部品や回路のパーツを作り、それらを組み立てる直前のもの」を購入しました。   ChaiBio https://www.chaibio.com/   OpenPCR https://www.chaibio.com/products/openpcr   なぜか http://openpcr.org/  で購入できなかったので、eBay にある ChaiBio で買いました。   OpenPCR - eBay http://www.ebay.com/itm/111096418574   本体価格は

R でいまどきなパッケージ開発 (devtools, testthat, roxygen2)

追記 (2012/04/21): 以下のコードは S4 classes で書いていますが、R5 reference classes で書き直してみました。こちらもどうぞ。 http://blog.hackingisbelieving.org/2012/04/r5-reference-class-r-devtools-testthat.html R のパッケージ開発の情報があまりないので、自分はこんな感じでやってます、というのを書いてみます。パッケージ開発支援の devtools と単体テスト支援の testthat, そしてドキュメント生成支援の roxygen を使うのがいまどきっぽいです。 そもそもパッケージを作製しているひとをあまりみたことがないので、もっとこうすべき、というのがあれば教えてほしいです。 今回はデモケースとして S4 OOP で、Idol クラスを定義し、とある身体的特徴の統計量を計算するパッケージを作ります。R のプロンプトは > で、シェルのプロンプトは $ で示しています。 0. 準備 必要になるパッケージをインストールします。 $ sudo R > install.packages(devtools) > install.packages(testthat) > q() devtools の設定をします。~/.Rpackages に設定を記述します。 $ emacs ~/.Rpackages list(   default = function(x) {     file.path("~/Project/dev/R/", x, x)   },   "idol" = "~/Projects/dev/R/idol/idol" ) 以下の行は今回パッケージを作製する作業ディレクトリになります。   "idol" = "~/Projects/dev/R/idol/idol" 1. ともあれ実装を始める 作業ディレクトリに移動します。 $mkdir -p ~/Project/dev/R/idol $ cd ~