復元抽出のアルゴリズム
もくてき
粒子フィルタ（パーティクル・フィルタ）を実行する際には、粒子のウェイト(weight)に比例する確率でリサンプリングを実行する必要がある。そのためのアルゴリズムとコードを考えたい。ここでは手元にある各粒子のウェイトはK個の要素からなるベクトルだと仮定して、さらにそれを復元を許してN個リサンプリングするという状況で考える。
こんな面倒な状況を考えなくても、要するにこれは
- いろんな色の球が入ってる壺から、１個適当に球を取り出して、その色をメモって、球を戻す
を複数回繰り返すことと同じで、高校生で習う確率の範囲で理解できる計算なわけだ。
アルゴリズム1(逆変換法)
逆変換法のアイディアを使って以下のようにするのが素朴なアイディアでコードも短い（後述）。
1: weightを確率に直し、その累積確率を計算し、これを{Qk, k=1,2,...K}とする
2: [0, 1]の実数値乱数rを一個生成する
3: rを、k=1,2,..., Kまで順番に累積確率Qkと比較していき、r
アルゴリズム2(ソートを使う方法)
ちょいと頭を使うといかのようなアルゴリズムもできる。
1: weightを確率に直し、その累積確率を計算し、これを{Qk, k=1,2,...K}とする
2: [0, 1]の実数値乱数rをN個生成し、{rn, n=1,2,...N}とする
3: {rn}∪{Qk}として、これをソートし、配列unionとする
4: ソートした配列unionの間に挟まってる"乱数rn"の個数を数える
5: 4で数えた個数ぶんだけ、該当するweightの要素を返却する
速度比較
上述のアルゴリズムのオーダーは、
- アルゴリズム1 : (K-1)*N（乱数との比較がボトルネック。NとKの二重ループが走る。）
- アルゴリズム2 : (N+K)log(N+K)（ソートがボトルネック。ソートにはn*log(n)オーダーの計算アルゴリズムを使う）
なので、データ(K) or リサンプリング(N)数が多い場合には、アルゴリズム2の方が効率的であると考えられるので、それを確かめてみる。
Rで書くと以下のような感じか。乱数は一度に生成するなどの計算時間節約は入れてしまっている。
# アルゴリズム1 resampling1 <- function(weight, size){ prob_cum <- cumsum(weight)/sum(weight) index <- sapply(runif(size), function(x)which(x < prob_cum)[1]) weight[index] } # アルゴリズム2 resampling2 <- function(weight, size){ prob_cum <- cumsum(weight)/sum(weight) size_weight <- length(weight) union <- c(prob_cum, runif(size)) union_indexes <- order(union, decreasing=TRUE) index <- numeric(size) value <- union_indexes[1] counter <- 1 for(union_index in union_indexes[-1]){ if(union_index > size_weight){ index[counter] <- value counter <- counter + 1 }else{ value <- union_index } } weight[index] }
速度の比較結果は
> x <- sample(1:10^4) > system.time(resampling1(x, 10^5)) ユーザ システム 経過 9.35 0.00 9.36 > system.time(resampling2(x, 10^5)) ユーザ システム 経過 0.25 0.00 0.25
というわけで、コードの長さとは裏腹にアルゴリズム２の方が遥かに速い。
・・・で、ここまで書いてから、Rに組み込みのsample関数のprob引数を与えれば、同じ処理ができることに気がついた。なので、アルゴリズム２とsample関数の計算速度を比べてみる。
> system.time(resampling2(x, 10^5)) ユーザ システム 経過 0.25 0.00 0.24 > system.time(sample(x, 10^5, replace=TRUE, prob=1:10^4)) ユーザ システム 経過 0 0 0
・・・Rの実装の方が圧倒的に速い！？！？！？！？！？
Rcppで書きなおす
きぃぃい、悔しい！！！負けるのは悔しいのでRcppでアルゴリズム２を書きなおしてみた。RcppとRcppArmadilloパッケージは
install.packages("Rcpp") install.packages("RcppArmadillo")
として突っ込んでおいておく。
library(Rcpp) sourceCpp(code=' #include <RcppArmadillo.h> using namespace Rcpp; // [[Rcpp::depends(RcppArmadillo)]] // [[Rcpp::export]] NumericVector resamplingCpp(arma::vec weight, int size){ RNGScope scope; const unsigned int size_weight = weight.size(); const arma::vec prob_cum = cumsum(weight)/sum(weight); arma::vec unionset = arma::join_cols(prob_cum, as<arma::vec>(runif(size))); arma::uvec union_indexes = arma::sort_index(unionset, "descend"); arma::uvec index = arma::uvec(size); double value = union_indexes[0]; int counter = 0; const arma::uvec::iterator union_indexes_begin = union_indexes.begin()+1; const arma::uvec::iterator union_indexes_end = union_indexes.end(); for(arma::uvec::iterator union_iterator=union_indexes_begin; union_iterator!=union_indexes_end; ++union_iterator){ if(*union_iterator > (size_weight-1)){ index[counter] = value; counter++; }else{ value = *union_iterator; } } arma::vec result = weight.elem(index); return NumericVector(result.begin(), result.end()); }')
そして、このresamplingCpp関数と、アルゴリズム２（R実装）・sample関数の速度比較してみる。
> system.time(resampling2(x, 10^6)) ユーザ システム 経過 2.90 0.00 2.93 > system.time(sample(x, 10^6, replace=TRUE, prob=1:10^4)) ユーザ システム 経過 0.03 0.02 0.04 > system.time(resamplingCpp(x, 10^6)) ユーザ システム 経過 0.22 0.00 0.22
うおお、なんてこった。元のコード（アルゴリズム２）よりも計算速度が10倍以上速くなったRcppのコードでも負けた…R恐るべし。俺のC++実装の問題な気もするのが、生のCっぽく書かくのはしんどいので、ここでおしまい。sample関数使えばいいや、もう。
答えの確認
ここで俺俺実装した結果が、ちゃんと同じ結果を返すことを確認しておく。
> weight <- sample(1:3) > set.seed(100) > x1 <- resampling1(weight, 10^2) > set.seed(100) > x2 <- resampling2(weight, 10^2) > set.seed(100) > x3 <- resamplingCpp(weight, 10^2) > all(x1[order(x1)]==x2[order(x2)]) [1] TRUE > all(x1[order(x1)]==x3[order(x3)]) [1] TRUE
上のように、乱数のシードを揃えれば確かに全ての結果が一致する。
sample関数については、実際に復元抽出した結果から確認しておく。
> val1 <- table(sample(1:5, 10^5, replace=TRUE, prob=1:5)) > val1/min(val1) 1 2 3 4 5 1.000000 2.026021 3.031165 4.059304 5.012103
確かに、モンテカルロ法の誤差のレベルで、指定したウェイト（1,2,3,4,5)に比例した答えになる。
非同次/確率的な強度λを持つ指数分布に従う確率変数τ(死亡/消滅時刻)の生成
非同次な指数分布
非同次なポアソン分布ってのがある。日本語で定義がちゃんと載ってるページがなかったので、英語版のwikipediaを参照すると
が該当する。要するにこれは、今まで、一定だと思っていたポアソン分布の強度(普通のポアソン分布だと平均だと思っていたパラメータ)を時間や空間に依存させてみたりするってもんだ。あと、この拡張として、この
自身を確率過程だと思っても以下の話は大体同じになる。正確な測度論の記述じゃないけど。
そして、
の「3.5 Process with Jumps」や「9.4 Credit Risk」を見ておくと良い。日本語だとたぶん信頼性工学系の書籍を漁ることになるんだとおもう、多分。
非同次な指数分布の数式展開
さて、(非)同次な指数分布に従う確率変数をと書くことにすると、まず、同次な指数分布の分布関数(時刻
までに死亡/消滅していない確率)は
と書くことができるが、今、非同次なケースを考えようとすると、強度が変数(時間・空間）依存性を持つので、以下のように積分形式で書くことができる。
これはが一定値だと思うと、ちょうど同次のケースに戻るので、これで定義としては良さそうだぞと、逆に、これを持って非同次な指数分布の分布関数とするぞとそういうことです。
ここで、ここから先、逆変換法を使っていくが、逆変換法については
Rで学ぶ逆変換（逆関数）法 from tera monagi
あたりを見てください*1。さて、非同次な指数分布に対して逆変換法を適用するために区間[0, 1]の一様乱数を導入し、
と書く。あとは、これをTについて解いてやればよいので、ちょっと変形して
とする。ここから先は解析的に解くのがキツいので、こいつはこのままで放置するとして、この数式を満たすような最小のが非同次な指数分布に従うというわけだ。もう少しちゃんと書くと
という形になる。更に、ここで、
は逆変換法の視点から解釈するに、とした指数分布になっているので、それを
と書くことにすると、
となる。なので、シミュレーションの手順としては
1: (平均1)とした指数分布から一個乱数
を取得する
2: の計算を数値積分なり解析解で頑張って計算する
3: を満たすような
を求めて、こいつを非同次な指数分布からの乱数だと思う
という手順になる。
Rでやってみる
というわけで、上で書いたアルゴリズムを実際に書いてみた。ここでは、
と設定している。深い意味はなく、解析的に計算したかったのでこうしている。この関数の積分は解析的に計算できるので、全体として確率密度関数は
となるんで、これとサンプリング結果を比較する。結果をPLOTすると、以下のようにほぼ一致しているので、これで良さそうだ。
ソースは以下。コメントにあるように平均1の指数分布と強度の積分が一致する点はめんどくさいのでソルバー(uniroot)関数に投げちゃってる。あと、乱数生成は汎用的にしたかったので数値積分で処理しているので遅いぞと。
#int_0^t \lambda(u)du - \xi #な関数。無理やり根を探す形で書くのでこうした objective <- function(t, lambda, xi){ integrate(lambda, 0, t)$value - xi } #適当な強度の関数形 lambda <- function(t){5*3*t^2+1} #描画範囲 range_plot <- c(0, 2) #サンプリング set.seed(100) x <- sapply( 1:10^4, function(i){ uniroot(objective, range_plot, tol=10^(-3), lambda=lambda, xi=rexp(1))$root } ) #解析解との比較 #上のラムダを積分（手でがんばる） lambda_integrated <- function(t){5*t^3+t} #解析解（確率密度関数） analytical_solution <- function(t){ lambda(t)*exp(-lambda_integrated(t)) } #比較PLOT library(ggplot2) qplot(x, geom="blank") + geom_histogram(aes(y=..density..), fill="aquamarine",color="black") + stat_function(fun=analytical_solution, color="red", size=1) + scale_x_continuous(limits=range_plot)
## 追記
「要するにCox回帰モデルで作る指数分布ってこと？」という突っ込みをいただきましたが、その解釈でもよいと思います。なので、Cox回帰モデルに従うサンプルデータを逆に生成できるぞと、データの水増しに使うことができるぞとそういうことです。
*1:他にいい資料がググったら出てくるかもしれないが、知らん
やってみよう分析！F#で強化学習(Q-learning, ε-greedy行動選択)
- やってみよう分析！Rで強化学習(Q-learning, ε-greedy行動選択) - My Life as a Mock Quant
- やってみよう分析！おまけ 2 - 1： Excel VBAで強化学習(Q-learning, ε-greedy / softmax 行動選択)
のF#実装版。強化学習自体の解説は上の記事読んどいたらいい。めんどいのでとりあえずε-greedyのみ。計算結果は
seq [[(1, 1), 4.99999999999999]; [(1, 2), 20]; [(2, 1), 9.99999999999999]; [(2, 2), 4.99984834210433]]
なんで、他２つとあってる。途中で力尽きて雑になったコードを直したい。。。そして、俺はF#にbreakがないって初めて知った。再帰で書けってことか。
open System open System.Collections.Generic open Microsoft.FSharp.Collections //Adhoc... type State = int type Action = int type QTable = Dictionary<State*Action, float> //General Code let filteredQValuesArray (state: State) (qTable: QTable) = qTable|> Seq.filter (fun (KeyValue(k,v)) -> fst k = state) |> Seq.map (fun (KeyValue(k, v)) -> (snd k, v)) |> Seq.toArray let selectAction (random: Random) (state: State) (qTable: QTable) = let qValuesAtState = filteredQValuesArray state qTable let qValues = qValuesAtState |> Array.map snd let maxQValue = qValues |> Array.max let maxIndex = [|0..(Array.length qValuesAtState - 1)|] |> Array.filter (fun i -> qValues.[i]=maxQValue) let size = Array.length maxIndex let selectedIndex = if size = 1 then maxIndex.[0] else random.Next(size) qValuesAtState.[selectedIndex] |> fst let epsilonGreedy (random: Random) epsilon (state: State) (qTable: QTable) = if random.NextDouble() < epsilon then let qValuesAtState = filteredQValuesArray state qTable let index = random.Next(Array.length qValuesAtState) qValuesAtState.[index] |> fst else selectAction random state qTable //Adhoc... let nextState (state: State) (action: Action) : State = if int action = 1 then if int state = 1 then 2 else 1 else state let calcReward (state: State) (action: Action) net = if int action <> 1 && int state = 1 then net else 0.0 let qLearning() = let random = new Random() let numA = 2 let numS = 2 let alpha = 0.3 let gamma = 0.5 let epsilon = 0.1 let trialMax = 10000 let stepMax = 10 let net = 10.0 let mutable state : State = 1 let qTable = new QTable() for s in [1..numS] do for a in [1..numA] do qTable.[(s, a)] <- 0.0 for i in [1..trialMax] do let j = ref 1 let mutable loopFlag = true state <- 1 while loopFlag && !j < stepMax do let action = epsilonGreedy random epsilon state qTable let stateNext = nextState state action let reward = calcReward state action net let qMax = qTable |> Seq.filter (fun (KeyValue(k,v)) -> fst k = stateNext) |> Seq.map (fun (KeyValue(k,v)) -> v) |> Seq.max qTable.[(state, action)] <- (1.0 - alpha)*qTable.[(state, action)] + alpha*(reward + gamma*qMax) if reward > 0.0 then loopFlag <- false else state <- stateNext incr j qTable [<EntryPoint>] let main argv = qLearning() |> printfn "%A" 0
(正規)分布の確率密度関数を乱数から推定する
すごい基本的なお話な気がするけど、数回は書き直しているのでメモっておく。ここでは
の２つの方法を比較して、それらが一致することを確認する。コードは以下。
#乱数のシードの指定 set.seed(100) #確率密度関数が一定とみなせるような適当な幅 dx <- 0.1 #標準偏差3の乱数生成 x <- rnorm(10^6, sd=3) #5 ± dxに入っている|ないを計算 hit <- sapply(x, function(z)abs(z-5) < dx) #確率を計算 sum(hit)/length(hit) #幅×確率密度関数 (2*dx)*dnorm(5, sd=3)
最後の部分だけ実行結果を載せると
> sum(hit)/length(hit) [1] 0.006558 > (2*dx)*dnorm(5, sd=3) [1] 0.006631809
リカレント・ニューラル・ネットワーク(Recurrent Neural Network, RNN)っていう、時系列に対するDeep Learningやりたい
はじめに
前々からDeep Learning系、特にその時系列への応用に興味があって、それにはリカレント・ニューラル・ネットワーク(Recurrent Neural Network, 以下RNN)ってのを理解する必要があるのは知っていたが、@yamano氏から
@teramonagi
こちらです。まずは手軽に使わせていただきました。
http://t.co/WTWmePJ0Hc
— バイアス撲滅委員会 (@yamano357) 2015, 1月 8というありがたいアドバイスを頂戴した。それに基づいてググりはじめると、Stuttgart(シュトゥットガルト)大学というドイツにある、たぶん日本でいう東工大的なポジションの大学が開発している、ニューラルネット関連の実装をたくさん詰め込んだライブラリ
が元々あって、それをラップしてRから叩けるようにしたものが、
というライブラリということがわかった。
RSNNSパッケージを使ってみる
上に書いたように、このパッケージ(ライブラリ)に入っている実装は、RNNに限らないようだが、興味あるのはRNNなのでそれをやる。Rのパッケージにくっついてくるマニュアルを読み込むとどうやら「jordan」関数ってのがあって、これを使うとよさそうだ。似たようなものとしてElman networkというのもあるらしい。詳細は原著
- Jordan, M. I. (1986), ’Serial Order: A Parallel, Distributed Processing Approach’, Advances in Connectionist Theory Speech 121(ICS-8604), 471-495.
- Elman, J. L. (1990), ’Finding structure in time’, Cognitive Science 14(2), 179–211.
を漁るか、Stuttgart Neural Network Simulatorのマニュアルを読まないとだめだな。どちらも「partially recurrent networks」と記述されていて、「なんでpartially？」ってのが気になってる。まぁ、とりあえず
- 興味持つ⇒動かす⇒ロジックを調べる⇒より興味を持つ⇒より興味を持って動かす⇒・・・
ってことで、パッケージに入っている例のコードを少し改変して動かしてみる。
#パッケージ入れる install.packages("RSNNS") library(RSNNS) library(dplyr) #データ読み出す data(snnsData) #snnsDataデータを入力と出力で分解 inputs <- snnsData$laser_1000.pat[,1] outputs <- snnsData$laser_1000.pat[,2] #更にトレーニング・ patterns <- splitForTrainingAndTest(inputs, outputs, ratio=0.15) #size:隠れ層のユニット数 #learnFuncParams:学習関数に対するパラメーター #linOut:活性化関数をリニア(TRUE)にするかロジスティックにするか(FALSE) modelJordan <- patterns %>% with({ jordan(inputsTrain, targetsTrain, size=c(8), learnFuncParams=c(0.1), maxit=100, inputsTest=inputsTest, targetsTest=targetsTest, linOut=FALSE)})
ここで、デフォルトの学習関数は
によると、
- JE_BP: Standard Backpropagation for partial recurrent networks
となっているが、これが単純なバックプロパゲーションでやってますよーという意味だと思われる。この辺を理解するにもRのパッケージじゃなくて直接Stuttgart Neural Network Simulatorのマニュアルを読まないとだめだな。で、この学習に使われるパラメーターがlearnFuncParamsで設定していると。
また、ここで使っている元々のデータは「カオス状態にあるアンモニア(NH3)レーザーの強度のデータ」だそうな。これは、公式サイトからソースを落としてきて、その中に入ってる「laser.README」というファイルに記述がある。
計算結果(modelJordan変数）を描画する用のplot関数も入っているが、中のコードを見た感じ通常のPLOT関数＋αな印象。また、フィッティング結果が結構上方にバイアスかかってるように見えるんだが、こんなもんなのだろうか？
#適当に結果を描画 plotIterativeError(modelJordan) #1 plotRegressionError(patterns$targetsTrain, modelJordan$fitted.values) #2 plotRegressionError(patterns$targetsTest, modelJordan$fittedTestValues) #3 plot(inputs[1:100], type="l") #4 lines(outputs[1:100], col="red") lines(modelJordan$fitted.values[1:100], col="green")
#1の結果（繰り返し毎に加重残差二乗和が減っていきますよー的なもの）
#2の結果（トレーニングデータ(実データ) v.s. その予測値。黒線は45度線(100%HITの場合これ)、赤線は(予測~実データの単回帰直線))
#3の結果（テストデータ(実データ) v.s. その予測値。黒線は45度線(100%HITの場合これ)、赤線は(予測~実データの単回帰直線))
#4の結果（黒線：入力値、赤線：実データ(出力)、緑線：予測値）
二次元ヒートカーネル(熱核)の時間発展の可視化
二次元ヒートカーネルなんていうと難しく聞こえるけど、要するにラプラス方程式の解である、時間発展するガウシアンの可視化です。物理系だと微分方程式の基本解という意味で、熱核じゃなくてグリーン関数って呼んでるものと等価（なはず）。
具体的には
1. 二次元のドリフト付ブラウン運動を複数生成する
2. 1の結果に対し、各時点において二次元カーネル密度推定を実行
3. 2.の結果をアニメーションとして可視化
という流れになる。まず必要となるライブラリを読み込んでおく。なかったらinstall.packages関数しましょう。
library(dplyr) library(MASS) library(pipeR) library(animation)
次に二次元ドリフト付ブラウン運動を生成させる。これは大体以下のコードを使っている。
#二次元ドリフト付ブラウン運動 #タイムステップ数 size <- 35 #適当な運動を作る randomwalk2D <- function(){cbind(x=cumsum(1+rnorm(size)), y=cumsum(1+rnorm(size)))} #本シミュレーションを走らせて結果をdata.frameに統合 df <- lapply(1:100, function(i)data.frame(randomwalk2D(), group=LETTERS[i], time_step=1:size)) %>>% rbind_all
次に、後の描画範囲指定で使う上下限の範囲を設定しつつ、
でやったように、二次元カーネル回帰させる。
#x, yの上下限計算 xlim <- df %>>% (x) %>>% range ylim <- df %>>% (y) %>>% range #カーネル回帰 kde2ds <- lapply(1:size, function(i){ df %>>% filter(time_step==i) %>>% with({ #バンド幅の計算 band_width <- c(bandwidth.nrd(x),bandwidth.nrd(y)) #二次元カーネル密度推定 kde2d(x, y, band_width, n=50, lims=c(xlim, ylim)) }) })
んで、animationパッケージの中ではImageMagickを使うので、それはそれで別途いれておく。入れておくとして、Rとの連携の際には一応こんなかんじでパスを設定しておくとうまくいった。昔使った時にはいらなかった記憶があるんだけどなぁ…
#animationパッケージ用の設定 ani.options(convert = 'C:/Program Files/ImageMagick-6.9.0-Q16/convert.exe')
以下のコードで二次元ヒートカーネル(熱核)の時間発展の可視化を実行している。z方向の値(確率密度)がすべての時間で共通になるようにうまくスケーリングしている。パラメーターの設定は黒魔術的に適当にやった。
#スケーリングあり ani.record(reset=TRUE) zlim_max <- c(kde2ds[[1]]$z, kde2ds[[size]]$z) %>>% max %>>% `*`(0.1) for(i in 1:size){ kde2ds[[i]] %>>% image(xlim=xlim, ylim=ylim, breaks=seq(from=0, to=zlim_max, length.out=13)) ani.record() } ani.options(interval = 0.1) ani.replay() saveGIF(ani.replay(), movie.name="2dheatkernel_scale.gif")
同様に今度はスケーリングなしのケース。従って、、「(異なる時点(アニメーションの1コマ))における図での明るさが等しい = 確率密度が等しい」とはならない。
#スケーリングなし ani.record(reset = TRUE) for(i in 1:size){ kde2ds[[i]] %>>% image(xlim=xlim, ylim=ylim) ani.record() } ani.options(interval = 0.1) ani.replay() saveGIF(ani.replay(), movie.name="2dheatkernel_noscale.gif")
頑張って3Dにしてみた。色がなんかおかしいぞ・・・？
pal <- colorRampPalette(c("red", "yellow"))(50) animation3d <- function(){ zmax <- max(kde2ds[[1]]$z) for(i in 1:size){ kde2ds[[i]] %>>% with({ col.ind <- cut(z,50) persp(x,y,z, zlim=c(0, zmax), col=pal[col.ind], theta=45)}) } } saveGIF(animation3d(), movie.name="2dheatkernel_3d.gif", interval=0.1)
円内に一様分布する乱数を生成する時は、俺、絶対忘れないよヤコビアンのこと
頭出し
ある半径Rの円内に一様分布する乱数を生成する時には注意しないといけないことがありますよというお話。所謂「一度はやってしまうミス」系でもある。この手の話は円に限ったわけではなく、円の高次元版である球、あるいは超球(次元>3)、あるいは任意の座標変換をかませてそこにヤコビアンが出て来るときでも同じ。
使うライブラリは以下の２つ。なければinstall.packages関数でインストールしておく。
library(ggplot2) library(dplyr)
また、以下のように定数を２つ定義しておく。意味はコメントにある通りだ。
#サンプル数 N <- 10^4 #半径のサイズ R <- 4
本題
さて、問題のある半径Rの円内に一様分布する乱数を生成するにはどうしたらいいのかというと、非常に単純に考えた場合、以下のように(俺は)思考する。
- X方向の成分として[-R, R]の間の一様乱数を生成する
- Y方向の成分として[-R, R]の間の一様乱数を生成する
- X^2 + Y^2 < R^2を満たすようなデータだけを残す
こうすることで、半径R内に一様に分布した乱数が得られそうだぞと、で、実際にやってみる。
#一様乱数から生成 df <- data.frame(x=runif(N, min=-R, max=R), y=runif(N, min=-R, max=R)) %>% filter(sqrt(x^2+y^2)<R) ggplot(df, aes(x=x, y=y)) + geom_point(color="blue")
散布図を見る限り、確かに一様に分布していそうだぞと、これでいいぞと。
ただ、このやり方はちょっとだけ無駄があって、「X^2 + Y^2 < R^2」を満たすデータのみを取得するというフィルタリングの処理を噛ませているために、発生させた乱数の全てを使っているわけではなく、一部（円からはみ出てしまった部分）を捨ててしまっている。そこでもうちょっと効率的なやり方を考えたい。高校数学でお勉強した"極座標"を思い出すとx, y座標と動径方向r、角度Θは
という関係で結びつけられていた。従って、円内に一様に分布する乱数を得るために
- 動径方向rの成分として[0, R]の間の一様乱数を生成する
- 角度Θ方向の成分として[0, 2π]の間の一様乱数を生成する
- XとYを上記の数式に従って変換し、生成する
というようにしてやれば半径R内に一様に分布した乱数が得られそうだぞと、で、実際にやってみる。
#動径・角度方向に分けてサンプリング theta <- runif(N, min=0, max=2*pi) r <- runif(N, min=0, max=R) df <- data.frame(x=r*cos(theta), y=r*sin(theta)) ggplot(df, aes(x=x, y=y)) + geom_point(color="blue")
あ、あれ？この結果の図を見ると明らかに中心部分の方にサンプルが集中していて、外側が手薄になっている、つまり、一様に分布していないぞと、これはよくないぞということが散布図からわかる。
この原因は今度は大学の数学で習う、重積分の変数変換公式における面積(要)素の変換、
のrを完全に見ていない、つまりそのヤコビアンを無視しているからだ*1。従って、動径方向に対しては半径rに比例するように乱数を引いてやらないといけなくて、[0, R]の一様な乱数uを逆関数数を用いて、
と変換して動径方向の乱数rを生成すればよいことがわかる。これは動径方向の累積分布関数P(R)が
とかける事から、このP(R)の逆関数をとることで証明される。
これを実際にやってみると
theta <- runif(N, min=0, max=2*pi) r <- sqrt(2*runif(N, min=0, max=R)) df <- data.frame(x=r*cos(theta), y=r*sin(theta)) ggplot(df, aes(x=x, y=y)) + geom_point(color="blue")
散布図より、一様に分布してるのがわかる。
追記(2015/1/18)
ほんぎゃああああああああああああああ！！！えらいミスしてて、上の図を良く見ると半径３くらいの円にしかなってなくて、俺が欲しかった、半径４の円になっていない。これを修正するためには、上のコードの
r <- sqrt(2*runif(N, min=0, max=R))
となっている箇所を、
r <- sqrt(2*runif(N, min=0, max=0.5*R^2))
とすればOK。これは、uの最大値とRの最大値を調正するためにmax引数をいじる必要があったということだ。これでPLOTすれば心穏やかに正しい結果となる。
*1:そもそも"一様"というのは考えている空間における抽象的な意味での面積素上に一様に分布するということなんで、こいつを無視してはいかん
