プログラム書法-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 = "保存量")

プログラム書法-2

コンピュータ上の実数は近似値である。ただし，2 のべき乗数の和で表される数は正確な値を保持している。整数は全て2のべき乗数の和であるから正確な値を保持している。

> options(digits=20)
> x <- 0.001 # 0.001 は 10 進数では正確な値であるが，2 進数では循環小数となり，正確な値を保持していない
> s <- 0
> for (i in 1:1000) s <- s+x
> s
[1] 1.0000000000000006661
> s == 1 # よって，0.001 を 1000 回足しても 1 にはならない
[1] FALSE

> (x <- 2^(-10)) # 2 の マイナス 10 乗は 10 進数でも，2 進数でも正確な値を保持している
[1] 0.0009765625
> s <- 0
> for (i in 1:1024) s <- s+x
> s
[1] 1
> s == 1 # 1024 回加えると，正確に 1 になる
[1] TRUE

結論：実数変数がある値と等しいかどうかを判定するようなプログラムを書いてはいけない

R 用の indent その3

たぶん，最終的に用いているのは皆同じプログラムなんだろうけど，Sweave を使う例も挙げておこう。

<<>>=
test<-function(x,y,z) {
n<-length(x) # 日本語のコメント
result<-numeric(n)
for (i in 1:n) {
result[i]<-(x[i]+y[i])*z[i]
}
return(result)
}
@

を含む *.Rnw を \SweaveOpts{keep.source=FALSE} で Sweave すると，以下のようなチャンクを含む *.tex ファイル（もちろん最終的には，*.pdf ファイルなど）ができる。
">" や "+"  は options(prompt = " ", continue = " ") で，空白にすることができる（空にはできないようだ）。
ただしこの場合も，コメント行はすっかりとなくなる。

> test <- function(x, y, z) {
+     n <- length(x)
+     result <- numeric(n)
+     for (i in 1:n) {
+         result[i] <- (x[i] + y[i]) * z[i]
+     }
+     return(result)
+ }

 

R 用の indent その2

formatR パッケージの tidy.souce 関数を使う。この方法では，コメントはちゃんと残る。

元のプログラム（encoding: utf-8）

test<-function(x,y,z) {
n<-length(x) # 日本語のコメント
result<-numeric(n)
for (i in 1:n) {
result[i]<-(x[i]+y[i])*z[i]
}
return(result)
}

> library(formatR)
> tidy.source("testsource.R")
test <- function(x, y, z) {
    n <- length(x)  # 日本語のコメント
    result <- numeric(n)
    for (i in 1:n) {
        result[i] <- (x[i] + y[i]) * z[i]

    }
    return(result)
}

整形すべきプログラムをクリップボードにコピーして，tidy.source() とするのが最も簡単。

tidy.dir は，指定したディレクトリ中の全ての R ファイルを tidy.source を使って整形する

tidy.gui は，tidy.source を呼ぶ GUI システムを提供する

R 用の indent その1

source("プログラムファイル.R", keep.source=FALSE)
dump("関数名", file="新プログラムファイル.R")

難点　元のプログラムについている注釈が全部なくなる

元のプログラム　testsource.R

test<-function(x,y,z){
n<-length(x)
result<-numeric(n)
for(i in 1:n){result[i]<-(x[i]+y[i])*z[i]}
return(result)}

source("testsource.R", keep.source=FALSE)
dump("test", file="new.testsource.R")

インデントされたプログラム new.testsource.R

test <-
function (x, y, z)
{
    n <- length(x)
    result <- numeric(n)
    for (i in 1:n) {
        result[i] <- (x[i] + y[i]) * z[i]
    }
    return(result)
}

同じようなことであるが，以下のようにしてもよい

> source("testsource.R", keep.source=TRUE)
> print(test)
function(x,y,z){
n<-length(x)
result<-numeric(n)
for(i in 1:n){result[i]<-(x[i]+y[i])*z[i]}
return(result)}
> print(test, useSource=FALSE)
function (x, y, z)
{
    n <- length(x)
    result <- numeric(n)
    for (i in 1:n) {
        result[i] <- (x[i] + y[i]) * z[i]
    }
    return(result)
}

プログラム書法-1

定数の引用

その言語処理システムで定義されているなら，それを使う。

R では，円周率は pi で引用できる。

もし定義されていない場合には，pi <- atan(1)*4 で正確に定義できる。

pi <- 3.14159  などと，いい加減に定義しないこと。

ネイピア定数は定義されていないけど e <- exp(1) とすればよい。

プログラミング書法

Rプログラミングのレクチャーなどするお方は，標準的なプログラミング書法に従わないと，受講生は皆，悪しきお手本に従って，悪しきプログラミングをする ようになるという実例があるわけだから，早々に改める方がよいと思うけどなあ。なぜ，悪しきプログラミング書法に固執するのか訳が分からない。

単に，頑迷固陋なのだろうか。そんなお年でもないと思うけど。

参考

Google's R Style Guide