n<-3
ds<-runif(n)
us<-sample(1:10,n)
Niter<-100000
dt<-0.01
m<-matrix(0,Niter,n)
p<-rep(0,Niter)
m[1,]<-runif(n)
p[1]<-0
for(i in 1:n){
p[1]<-p[1]+us[i]*m[1,i]
}
#pが1からスタートするように補正
m[1,]<-m[1,]/p[1]
p[1]<-0
for(i in 1:n){
p[1]<-p[1]+us[i]*m[1,i]
}
for(i in 2:Niter){
xyz<-1
for(j in 1:n){
xyz<-xyz*m[i-1,j]
}
for(j in 1:n){
pre<-j-1
if(pre<1)pre<-n
post<-j+1
if(post>n)post<-1
m[i,j]<-m[i-1,j]+dt*xyz*(-ds[j]/m[i-1,pre] + us[pre]/us[j]*ds[pre]/m[i-1,post])
}
for(j in 1:n){
p[i]<-p[i]+us[j]*m[i,j]
}
}
par(mfcol=c(2,2))
plot(m[,1],m[,2],type="l",main="1 vs. 2")
plot(m[,2],m[,3],type="l",main="2 vs. 3")
plot(m[,3],m[,1],type="l",main="3 vs. 1")
plot(p,ylim=c(0,2),main="保存量")
これを書き換える。
1. ベクトル演算を活用する(ループを使わない)
2. R にはスカラーはないが,添え字が不要なときにはベクトル・行列は使わない
3. for ループの中で不変であるものは,外に出す
dt*xyz の掛け算は for (j の外で計算する
p[i] の計算は for (i の外で計算する
注:pre, post を計算する箇所はこのままの方が速度は速い(フォーマットは変える)
以下のようになる。速度は旧版の 2 割り増し。
n <- 3
ds <- runif(n)
us <- sample(10, n) # sample(1:10, n) と同じだけど
Niter <- 100000L # あまり関係ないけど,整数は L を付ける
dt <- 0.01
m <- matrix(0, Niter, n)
m1 <- runif(n)
p1 <- sum(us*m1) # ベクトル演算を活用(ループを使わない)
# p が 1 からスタートするように補正
m[1, ] <- m1 / p1
for (i in 2:Niter) {
dt.xyz <- dt * prod(m[i - 1, ]) # ループの外で計算, ベクトル演算を活用する(ループを使わない)
for (j in 1:n) {
pre <- j-1
if (pre < 1) {
pre <- n
}
post <- j+1
if (post > n) {
post <- 1
}
m[i, j] <- m[i - 1, j] + dt.xyz * (-ds[j]/m[i - 1, pre] + us[pre]/us[j] * ds[pre]/m[i - 1, post])
}
}
p <- m %*% us # まとめて計算,行列掛け算を使う
par(mfcol = c(2, 2))
plot(m[, 1], m[, 2], type = "l", main = "1 vs. 2")
plot(m[, 2], m[, 3], type = "l", main = "2 vs. 3")
plot(m[, 3], m[, 1], type = "l", main = "3 vs. 1")
plot(p, ylim = c(0, 2), main = "保存量")
※コメント投稿者のブログIDはブログ作成者のみに通知されます