Rコードの最適化例：レスリー行列による個体群成長
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
検索
|
最終更新
|
ヘルプ
]
開始行:
SIZE(18){COLOR(red){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....
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=TRU...
-関数の最後のコメントは R のヒストリ機能と併用すると、す...
-長い(そして多重の)ループには閉じたところにコメントしてお...
-ある関数でだけ使う定数は関数中(冒頭)に定義しておく(ばら...
COLOR(magenta){TODOリスト}
- j (そして i)ループは下請け関数化が簡潔さの観点から望ま...
- trial ループは *apply 関数が定番のコード？
- 変更可能定数の全面引数化？
- 結果だけでなく、シミュレーション条件をリスト返り値で保...
## 改訂例(1)
growth.1 <- function (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,...
y0 <- c(10,10,10,10,30)*n # 各世代...
yy0 <- numeric(5) # 作業用...
total.y <- matrix(0, nr = 25, nc = 50) # 25回の...
for (trial in 1:25){ # 25 回の...
y <- y0 # 何度も...
for (t in 1:50){ # 50 期間...
yy <- yy0 # 作業変...
for (j in 1:5) {
for (i in 1:5) yy[j] <- yy[j] + rbinom(1,y[i],x[...
}
total.y[trial, t] <- sum(yy) # trial ...
y <- yy # y を yy...
} # t ループの終り
} # trial ループの終り
matplot(t(total.y), type="l") # 複数の系列を同時にプ...
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)を付け加えて下さい。}
//-失礼しました。ご指摘の幾つかの件見直し中に気づき、変更...
-なお、コードクリニックを開業する気はありませんのであしか...
-仕事中に覗いたら、なんと！ビックリしました。と、同時に感...
-改訂版(1) の二行 total.y[trial, t] <- sum(yy); y <- yy ...
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) { # 引数...
total.y <- matrix(0, nr = 25, nc = 50) # 25回...
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,...
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 ループの終り
} # trial ループの終り
matplot(t(total.y), type="l") # 複数の系列を同時にプ...
invisible(total.y) # 返り値(変数に付値しな...
}
# total.y <- growth.2(1)
-改訂例(2)でも，t ループの中の3行は，total.y[trial, t] <-...
## 改訂例(2') 改訂例(1)の二倍ほどの時間がかかるので注意
## レスリー行列(x), 初期人口(y0), 試行回数(TRIAL), 期間...
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回の結...
# 改訂例...
for (trial in 1:TRIAL) { # TRIAL 回のシ...
y <- y0 # 何度も使う定数は変...
for (t in 1:T) { # T 期間経過させる
total.y[t,trial] <- sum(y<-colSums(matrix(rbinom...
}
}
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, b...
y <- c(10,10,10,10,30) # 各世代...
total.y <- growth.2(x, y)
~
~
COLOR(magenta){改訂例(3)} 引数大増量、返り値しっかり版 (...
作業用定数ベクトル yy0 をなくした(0*y とすれば良し、yy <-...
# 引数 LeslieMatrix レスリー行列
# init.pop 世代別初期人口ベクトル
# period 経過期間 (既定値 50)
# sim.no シミュレーション回数 (既定値 25)
growth.3 <- function (LeslieMatrix, init.pop, period = 5...
ngen <- length(init.pop) # 世代数
total.y <- array(0, c(sim.no, period)) # sim....
for (trial in 1:sim.no){ # sim....
y <- init.pop # 作業...
for (t in 1:period){ # peri...
yy <- 0*y # 作業...
for (j in 1:ngen)
for (i in 1:ngen) yy[j] <- yy[j] + rbinom(1,y[i]...
total.y[trial, t] <- sum(y <- yy) # tria...
} # t ループ終り
} # trial ループ終り
matplot(t(total.y), type="l") # 複数の系列を同時にプ...
## 返り値(変数に付値しない限りコンソールには出力しない)
invisible( list(total.y=total.y, LM=LeslieMatrix,
ip=init.pop, period=period, sim.no=sim...
}
# res <- growth.3(LeslieMatrix, init.pop, period = 50, s...
> 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) # ...
> str(res) # ...
List of 5
$ total.y: num [1:25, 1:50] 81 100 89 82 89 87 86 76 89...
$ 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...
[20] 94 107 109 111 110 101 106 105 114 118 131 135 137...
[39] 161 174 168 168 158 144 160 166 169 171 163 149
> growth.3(LM, IP, period = 50, sim.no = 25) # 付値しな...
-期間経過初期で一旦総人口が落ち込むように見えるのはのはな...
-初期値に依存しているからです。その後は安定して増加してい...
-レスリー行列の一行目0, 0.57, 0.57, 0.57, 0.57は繁殖率で...
COLOR(magenta){改訂例(3-1)}
確率 0 の場合が多いのでパスするようにした (実行時間約半分)
# 引数 LeslieMatrix レスリー行列
# init.pop 世代別初期人口ベクトル
# period 経過期間 (既定値 50)
# sim.no シミュレーション回数 (既定値 25)
growth.31 <- function (LeslieMatrix, init.pop, period = ...
ngen <- length(init.pop) # 世代数
total.y <- array(0, c(period, sim.no)) # sim....
for (trial in 1:sim.no){ # sim....
y <- init.pop # 作業...
for (t in 1:period){ # peri...
yy <- 0*y # 作業...
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) # tria...
} # t ループ終り
} # trial ループ終り
matplot(total.y, type="l") # 複数の系列を同時にプロット
## 返り値(変数に付値しない限りコンソールには出力しない)
invisible( list(total.y=total.y, LM=LeslieMatrix,
ip=init.pop, period=period, sim.no=sim...
}
# res <- growth.31(LeslieMatrix, init.pop, period = 50, ...
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
終了行:
SIZE(18){COLOR(red){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....
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=TRU...
-関数の最後のコメントは R のヒストリ機能と併用すると、す...
-長い(そして多重の)ループには閉じたところにコメントしてお...
-ある関数でだけ使う定数は関数中(冒頭)に定義しておく(ばら...
COLOR(magenta){TODOリスト}
- j (そして i)ループは下請け関数化が簡潔さの観点から望ま...
- trial ループは *apply 関数が定番のコード？
- 変更可能定数の全面引数化？
- 結果だけでなく、シミュレーション条件をリスト返り値で保...
## 改訂例(1)
growth.1 <- function (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,...
y0 <- c(10,10,10,10,30)*n # 各世代...
yy0 <- numeric(5) # 作業用...
total.y <- matrix(0, nr = 25, nc = 50) # 25回の...
for (trial in 1:25){ # 25 回の...
y <- y0 # 何度も...
for (t in 1:50){ # 50 期間...
yy <- yy0 # 作業変...
for (j in 1:5) {
for (i in 1:5) yy[j] <- yy[j] + rbinom(1,y[i],x[...
}
total.y[trial, t] <- sum(yy) # trial ...
y <- yy # y を yy...
} # t ループの終り
} # trial ループの終り
matplot(t(total.y), type="l") # 複数の系列を同時にプ...
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)を付け加えて下さい。}
//-失礼しました。ご指摘の幾つかの件見直し中に気づき、変更...
-なお、コードクリニックを開業する気はありませんのであしか...
-仕事中に覗いたら、なんと！ビックリしました。と、同時に感...
-改訂版(1) の二行 total.y[trial, t] <- sum(yy); y <- yy ...
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) { # 引数...
total.y <- matrix(0, nr = 25, nc = 50) # 25回...
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,...
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 ループの終り
} # trial ループの終り
matplot(t(total.y), type="l") # 複数の系列を同時にプ...
invisible(total.y) # 返り値(変数に付値しな...
}
# total.y <- growth.2(1)
-改訂例(2)でも，t ループの中の3行は，total.y[trial, t] <-...
## 改訂例(2') 改訂例(1)の二倍ほどの時間がかかるので注意
## レスリー行列(x), 初期人口(y0), 試行回数(TRIAL), 期間...
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回の結...
# 改訂例...
for (trial in 1:TRIAL) { # TRIAL 回のシ...
y <- y0 # 何度も使う定数は変...
for (t in 1:T) { # T 期間経過させる
total.y[t,trial] <- sum(y<-colSums(matrix(rbinom...
}
}
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, b...
y <- c(10,10,10,10,30) # 各世代...
total.y <- growth.2(x, y)
~
~
COLOR(magenta){改訂例(3)} 引数大増量、返り値しっかり版 (...
作業用定数ベクトル yy0 をなくした(0*y とすれば良し、yy <-...
# 引数 LeslieMatrix レスリー行列
# init.pop 世代別初期人口ベクトル
# period 経過期間 (既定値 50)
# sim.no シミュレーション回数 (既定値 25)
growth.3 <- function (LeslieMatrix, init.pop, period = 5...
ngen <- length(init.pop) # 世代数
total.y <- array(0, c(sim.no, period)) # sim....
for (trial in 1:sim.no){ # sim....
y <- init.pop # 作業...
for (t in 1:period){ # peri...
yy <- 0*y # 作業...
for (j in 1:ngen)
for (i in 1:ngen) yy[j] <- yy[j] + rbinom(1,y[i]...
total.y[trial, t] <- sum(y <- yy) # tria...
} # t ループ終り
} # trial ループ終り
matplot(t(total.y), type="l") # 複数の系列を同時にプ...
## 返り値(変数に付値しない限りコンソールには出力しない)
invisible( list(total.y=total.y, LM=LeslieMatrix,
ip=init.pop, period=period, sim.no=sim...
}
# res <- growth.3(LeslieMatrix, init.pop, period = 50, s...
> 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) # ...
> str(res) # ...
List of 5
$ total.y: num [1:25, 1:50] 81 100 89 82 89 87 86 76 89...
$ 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...
[20] 94 107 109 111 110 101 106 105 114 118 131 135 137...
[39] 161 174 168 168 158 144 160 166 169 171 163 149
> growth.3(LM, IP, period = 50, sim.no = 25) # 付値しな...
-期間経過初期で一旦総人口が落ち込むように見えるのはのはな...
-初期値に依存しているからです。その後は安定して増加してい...
-レスリー行列の一行目0, 0.57, 0.57, 0.57, 0.57は繁殖率で...
COLOR(magenta){改訂例(3-1)}
確率 0 の場合が多いのでパスするようにした (実行時間約半分)
# 引数 LeslieMatrix レスリー行列
# init.pop 世代別初期人口ベクトル
# period 経過期間 (既定値 50)
# sim.no シミュレーション回数 (既定値 25)
growth.31 <- function (LeslieMatrix, init.pop, period = ...
ngen <- length(init.pop) # 世代数
total.y <- array(0, c(period, sim.no)) # sim....
for (trial in 1:sim.no){ # sim....
y <- init.pop # 作業...
for (t in 1:period){ # peri...
yy <- 0*y # 作業...
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) # tria...
} # t ループ終り
} # trial ループ終り
matplot(total.y, type="l") # 複数の系列を同時にプロット
## 返り値(変数に付値しない限りコンソールには出力しない)
invisible( list(total.y=total.y, LM=LeslieMatrix,
ip=init.pop, period=period, sim.no=sim...
}
# res <- growth.31(LeslieMatrix, init.pop, period = 50, ...
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
ページ名:
![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
