プログラム書法-3

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 = "保存量")