![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
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
}
}
改良が望ましい点のリスト
個人的趣味レベルの追加リスト
TODOリスト
## 改訂例(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(): The style ref(filename,pagename) is ambiguous and become obsolete. Please try ref(pagename/filename)
いまだに残る疑問点
是非皆さんの改訂版(x)を付け加えて下さい。
改訂例(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') 改訂例(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)
改訂例(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) # 付値しなければグラフ出力だけ
改訂例(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)
growth.1.png