Rコードの最適化例：レスリー行列による個体群成長の編集

SIZE(18){COLOR(red){Rコードの最適化例：レスリー行列による個体群成長}}  ([[Rコード最適化のコツと実例集]] に戻る)~

勝手ながら、Q&A (初級者コース) への「しまだ」さんのコードを
R コードの改良(おそらく最適化の趣旨からは外れますが)例に使わせて頂き
ます。オリジナルコードの欠点をあげつらうという趣旨ではありま
せん。あくまで参考になれば、ということですからご了解下さい。
もちろん個人の慣用・趣味レベルの改訂も入りますから、
「ここはこう直すべきだ」という意味ではありませんのでお間違えなく。
(ついでに、いまだに問題を完全には理解していないことをお断りしておきます。)

## しまださんのオリジナルコード
 for (trial in 1:25){
 #--
 x<-rbind(c(0,0.57,0.57,0.57,0.57),c(0.46,0,0,0,0),c(0,0.77,0,0,0),
          c(0,0,0.82,0,0),c(0,0,0,0.91,0.65))
 y<-cbind(c(10,10,10,10,30))
 yy<-cbind(c(1:5))
 totaly<-rbind(c(1:50))
 #--
 for (t in 1:50){
 for (j in 1:5){
 for (i in 1:5){
 yy[j]<-yy[j]+rbinom(1,y[i],x[i,j])
 }
 }
 #--replace:y[i]<-yy[ii]
 for (i in 1:5){
 y[i]<-yy[i]
 }
 #--replace:total_y<-y
 totaly[t]=colSums(y)
 #--reset:yy[i]<-0
 for (i in 1:5){
 yy[i]<-0
 }
 }
 #--graph plot
 par(new=T)
 plot(1:50,totaly,type="l")
 #--reset totaly
 for (i in 1:50){
 totaly[i]<-0
 }
 }

COLOR(magenta){改良が望ましい点のリスト}
-変数 y, yy は行列(縦ベクトル)にしなくても、ベクトルのままで良い(行列とのかけ算をしないので)
-y の更新は y <- yy 一発でOK
-yy の初期化は yy <- numeric(5) でOK(長さ 5 のすべて 0 のベクトルになる)
-複数の系列を同じスケールでプロットするには matplot 関数が便利
-一つ続きの作業は関数にしておくのが基本(編集にも便利)
-インデントはコードを見やすくし、間違いを見付けるにも便利。但しインデントのしかた(そしてブロック化用の波括弧の位置)は
宗教論争を招く恐れあり(私としてはできるだけ多くのコードがエディタの画面に収まるという規準を採用)
-せっかく計算した途中結果はとりあえず返り値にして保存が原則(デバッグ、二次的作業に便利)
-返り値が大量になる際は invisible 返り値にしておけば変数に付値しない限りうっかり画面に出力して慌てることがなく安心
-コメントはけちらず(そして日本語で)
-後で色々条件を変更して実行したければ、変更可能箇所は関数引数にしておくと便利（この例ではレスリー行列、世代数、初期人口、試行数、作図パラメータ等）
-c(1:5) は 1:5 だけでOK (結局同じですが)

COLOR(magenta){個人的趣味レベルの追加リスト}
-小さな行列は、行列の形で入力すると間違いが無い(byrow=TRUE を忘れない) 
-関数の最後のコメントは R のヒストリ機能と併用すると、すぐに関数を実行するのに便利(私はエディタ画面からカット＆ペーストすることが多いので)
-長い(そして多重の)ループには閉じたところにコメントしておくとデバッグに便利
-ある関数でだけ使う定数は関数中(冒頭)に定義しておく(ばらばらにならず、編集に便利)

COLOR(magenta){TODOリスト}
- j (そして i)ループは下請け関数化が簡潔さの観点から望ましい？
- trial ループは *apply 関数が定番のコード？
- 変更可能定数の全面引数化？
- 結果だけでなく、シミュレーション条件をリスト返り値で保存しておく(後で良くわからなくことがある)?

## 改訂例(1)
 growth.1 <- function (n = 1) {                    # 引数 n (既定値 1)は初期人口への倍数
   x <- matrix( c(0,    0.57, 0.57, 0.57, 0.57,    # レスリー行列(各世代の繁殖率と生存率)
                  0.46, 0,    0,    0,    0,       
                  0,    0.77, 0,    0,    0,
                  0,    0,    0.82, 0,    0,
                  0,    0,    0,    0.91, 0.65 ), nrow=5, byrow=TRUE)
   y0 <- c(10,10,10,10,30)*n                    # 各世代の初期人口
   yy0 <- numeric(5)                            # 作業用変数 yy の初期値定数
   total.y <- matrix(0, nr = 25, nc = 50)       # 25回の結果を記録する変数(0 で初期化)   
   for (trial in 1:25){                         # 25 回のシミュレーション
     y <- y0                                    # 何度も使う定数は変数にいれておき参照
     for (t in 1:50){                           # 50 期間経過させる
       yy <- yy0                                # 作業変数(0 で初期化) 
       for (j in 1:5) {
         for (i in 1:5) yy[j] <- yy[j] + rbinom(1,y[i],x[i,j]) 
       }
       total.y[trial, t] <- sum(yy)             # trial 回目の期間 t の結果を記録
       y <- yy                                  # y を yy で更新
     } # t ループの終り
   } # trial ループの終り
 
   matplot(t(total.y), type="l")  # 複数の系列を同時にプロット(行列 total.y を転置する必要)
   invisible(total.y)             # 返り値(変数に付値しない限りコンソールには出力しない)
 } 
 # total.y <- growth.1(1)

growth.1(1), growth.1(10), growth.1(100) の結果のグラフ
#ref(growth.1.png, left)

COLOR(magenta){いまだに残る疑問点}
-このままでは繁殖と死亡を区別していないような?(死亡率ではなく生存率だからこのままで良いのかな)
-初期人口が少なすぎるような？少集団の挙動に興味があるのか、それともいずれ
大集団で実行する計画なのか？ 後者の場合(そして比率を変えて何度も実行する気)なら、コードの高速化の余地。
-大集団でシミュレーション回数を増やした結果を時刻毎に平均すれば結局レスリー行列によるかけ算結果を近似する（大数の法則）？

COLOR(magenta){是非皆さんの改訂版(x)を付け加えて下さい。}
//-失礼しました。ご指摘の幾つかの件見直し中に気づき、変更作業中にコメントを頂くことになりました。 --  &new{2004-06-21 (月) 16:29:30};
-なお、コードクリニックを開業する気はありませんのであしからず。私が(そして他の人にも)面白いと思われるものだけを取り上げていく方針です。 --  &new{2004-06-21 (月) 16:33:31};
-仕事中に覗いたら、なんと！ビックリしました。と、同時に感謝します。やっぱし違いますね。勉強になります。 -- [[しまだ]] &new{2004-06-21 (月) 17:03:19};
-改訂版(1) の二行 total.y[trial, t] <- sum(yy); y <- yy は一行 total.y[trial, t] <- sum(y <- yy) で良い。R では付値演算子 <- はそれ自身(付値された値を返す)関数であるという事実を使う。初心者にはわかりにくいが、なれるととても便利。 --  &new{2004-06-22 (火) 10:09:24};

COLOR(magenta){改訂例(2)}

## 改訂例(2)　　改訂例(1)の二倍ほどの時間がかかるので注意
 rbinom2 <- function(n, p) { # 偽のベクトル化
 	sapply(1:length(p), function(i) rbinom(1, n[i], p[i]))
 }
 growth.2 <- function (n = 1) {                    # 引数 n (既定値 1)は初期人口への倍数
   total.y <- matrix(0, nr = 25, nc = 50)          # 25回の結果を記録する変数(0 で初期化)
   x <- matrix( c(0,    0.57, 0.57, 0.57, 0.57,    # レスリー行列(各世代の繁殖率と生存率)
                  0.46, 0,    0,    0,    0,       
                  0,    0.77, 0,    0,    0,
                  0,    0,    0.82, 0,    0,
                  0,    0,    0,    0.91, 0.65 ), nrow=5, byrow=TRUE)
   rnd <- matrix(0, nr=5, nc=5)
   y0 <- c(10,10,10,10,30)*n                    # 各世代の初期人口
   for (trial in 1:25){                         # 25 回のシミュレーション
     y <- y0                                    # 何度も使う定数は変数にいれておき参照
     for (t in 1:50){                           # 50 期間経過させる
       rnd <- matrix(rbinom2(rep(y, 5), x), nr=5, nc=5)
       y <- colSums(rnd)                        # y の更新
       total.y[trial, t] <- sum(y)              # trial 回目の時刻 t の結果を記録
     } # t ループの終り
   } # trial ループの終り
 
   matplot(t(total.y), type="l")  # 複数の系列を同時にプロット(行列 total.y を転置する必要)
   invisible(total.y)             # 返り値(変数に付値しない限りコンソールには出力しない)
 } 
 # total.y <- growth.2(1)
-改訂例(2)でも，t ループの中の3行は，total.y[trial, t] <- sum(y <- colSums(matrix(rbinom2(rep(y, 5), x), nr=5, nc=5))) 　とも書けますが（上の方の  rnd の定義も不要になる），上のような方がわかりやすいでしょうね。 --  &new{2004-06-22 (火) 10:28:39};

## 改訂例(2')　　改訂例(1)の二倍ほどの時間がかかるので注意
 ## レスリー行列(x), 初期人口(y0), 試行回数(TRIAL), 期間数(T) 
 growth.2 <- function (x, y0, TRIAL=25, T=50) {  
 	rbinom2 <- function(n, p) {              # 偽のベクトル化
 		sapply(1:length(p), function(i) rbinom(1, n[i], p[i]))
 	}
 	DIM <- nrow(x)
 	total.y <- matrix(0, nr=T, nc=TRIAL)     # TRIAL回の結果を記録する変数(0 で初期化) 
                                                 # 改訂例(1)とは行数・列数が逆
 	for (trial in 1:TRIAL) {                 # TRIAL 回のシミュレーション
 		y <- y0                          # 何度も使う定数は変数にいれておき参照
 		for (t in 1:T) {                 # T 期間経過させる
 		      total.y[t,trial] <- sum(y<-colSums(matrix(rbinom2(rep(y, DIM),x),nr=DIM))) 
 		}
 	}
 	matplot(total.y, type="l")               # 複数の系列を同時にプロット
 	invisible(total.y)                       # 返り値(変数に付値しない限りコンソールには出力しない)
 } 
 x <- matrix( c(0,    0.57, 0.57, 0.57, 0.57,    # レスリー行列(各世代の繁殖率と生存率)
                0.46, 0,    0,    0,    0,       
                0,    0.77, 0,    0,    0,
                0,    0,    0.82, 0,    0,
                0,    0,    0,    0.91, 0.65 ), nrow=5, byrow=TRUE)
 y <- c(10,10,10,10,30)                          # 各世代の初期人口
 total.y <- growth.2(x, y)
~
~
COLOR(magenta){改訂例(3)} 引数大増量、返り値しっかり版 (さらに引数の整合性チェック(存在、サイズ等)を含める余地)~
作業用定数ベクトル yy0 をなくした(0*y とすれば良し、yy <- numeric(5) より少し早くなる)

# 引数 LeslieMatrix レスリー行列
 #      init.pop     世代別初期人口ベクトル
 #      period       経過期間 (既定値 50)
 #      sim.no       シミュレーション回数 (既定値 25)
 growth.3 <- function (LeslieMatrix, init.pop, period = 50, sim.no = 25) {
   ngen <- length(init.pop)                        # 世代数
   total.y <- array(0, c(sim.no, period))          # sim.no 回の結果を記録する行列(0 で初期化)   
   for (trial in 1:sim.no){                        # sim.no 回のシミュレーション
     y <- init.pop                                 # 作業変数(init.pop で初期化) 
     for (t in 1:period){                          # period 期間経過させる
       yy <- 0*y                                   # 作業変数(0 で初期化) 
       for (j in 1:ngen) 
         for (i in 1:ngen) yy[j] <- yy[j] + rbinom(1,y[i],LeslieMatrix[i,j]) 
       total.y[trial, t] <- sum(y <- yy)           # trial 回目の期間 t の結果を記録、 y を yy で更新
     } # t ループ終り
   } # trial ループ終り
 
   matplot(t(total.y), type="l")  # 複数の系列を同時にプロット(行列 total.y を転置する必要)
   ## 返り値(変数に付値しない限りコンソールには出力しない)
   invisible( list(total.y=total.y, LM=LeslieMatrix, 
                   ip=init.pop, period=period, sim.no=sim.no))
 } 
 # res <- growth.3(LeslieMatrix, init.pop, period = 50, sim.no = 25)

> LM <- rbind(c( 0,    0.57, 0.57, 0.57, 0.57 ),      # レスリー行列(各世代の繁殖率と生存率)
               c( 0.46, 0,    0,    0,    0    ),
               c( 0,    0.77, 0,    0,    0    ),
               c( 0,    0,    0.82, 0,    0    ),
               c( 0,    0,    0,    0.91, 0.65 ) )
 > IP <- c(10,10,10,10,30)                             # 世代別初期人口
 > res <- growth.3(LM, IP, period = 50, sim.no = 25)   # 実行例 (結果を変数 res に保存)
 > str(res)                                            # 返り値リストの中身確認 
 List of 5
  $ total.y: num [1:25, 1:50] 81 100 89 82 89 87 86 76 89 94 ...
  $ LM     : num [1:5, 1:5] 0 0.46 0 0 0 0.57 0 0.77 0 0 ...
  $ ip     : num [1:5] 10 10 10 10 30
  $ period : num 50
  $ sim.no : num 25
 > res[[1]][1,]        # 例えば第一回シミュレーション結果は
  [1]  81  80  71  62  66  74  79  73  78  88  94  81  78  86  79  88  95  79  88
 [20]  94 107 109 111 110 101 106 105 114 118 131 135 137 146 139 139 144 148 155
 [39] 161 174 168 168 158 144 160 166 169 171 163 149
 > growth.3(LM, IP, period = 50, sim.no = 25) # 付値しなければグラフ出力だけ
 
-期間経過初期で一旦総人口が落ち込むように見えるのはのはなぜ？ --  &new{2004-06-22 (火) 23:23:55};
-初期値に依存しているからです。その後は安定して増加していきますので、本来なら安定した増加部分‐‐安定齢分布‐‐を初期値に使えば、綺麗に増加（もしくは減少）していきます。 -- [[しまだ]] &new{2004-06-23 (水) 23:59:19};
-レスリー行列の一行目0, 0.57, 0.57, 0.57, 0.57は繁殖率です。２行目以降は生存率です -- [[しまだ]] &new{2004-06-24 (木) 00:00:41};

COLOR(magenta){改訂例(3-1)}  
確率 0 の場合が多いのでパスするようにした (実行時間約半分)

# 引数 LeslieMatrix レスリー行列
 #      init.pop     世代別初期人口ベクトル
 #      period       経過期間 (既定値 50)
 #      sim.no       シミュレーション回数 (既定値 25)
 growth.31 <- function (LeslieMatrix, init.pop, period = 50, sim.no = 25) {
   ngen <- length(init.pop)                        # 世代数
   total.y <- array(0, c(period, sim.no))          # sim.no 回の結果を記録する行列(0 で初期化)   
   for (trial in 1:sim.no){                        # sim.no 回のシミュレーション
     y <- init.pop                                 # 作業変数(init.pop で初期化) 
     for (t in 1:period){                          # period 期間経過させる
       yy <- 0*y                                   # 作業変数(0 で初期化) 
       for (j in 1:ngen) 
         for (i in 1:ngen) 
           if ((p <- LeslieMatrix[i,j]) != 0)       # 確率零ならパス     
             yy[j] <- yy[j] + rbinom(1,y[i],p) 
       total.y[t, trial] <- sum(y <- yy)           # trial 回目の期間 t の結果を記録、 y を yy で更新
     } # t ループ終り
   } # trial ループ終り
 
   matplot(total.y, type="l")  # 複数の系列を同時にプロット
   ## 返り値(変数に付値しない限りコンソールには出力しない)
   invisible( list(total.y=total.y, LM=LeslieMatrix, 
                   ip=init.pop, period=period, sim.no=sim.no))
 } 
 # res <- growth.31(LeslieMatrix, init.pop, period = 50, sim.no = 25) 
 
 LM <- rbind(c( 0,    0.57, 0.57, 0.57, 0.57 ),      # レスリー行列(各世代の繁殖率と生存率)
             c( 0.46, 0,    0,    0,    0    ),
             c( 0,    0.77, 0,    0,    0    ),
             c( 0,    0,    0.82, 0,    0    ),
             c( 0,    0,    0,    0.91, 0.65 ) )
 IP <- c(10,10,10,10,30)                             # 世代別初期人口
 ## 実行例
 res <- growth.31(LM, IP)

#comment

タイムスタンプを変更しない