2013-06-01

Rのグラフにスプライン曲線を描きこむ

R言語

TTFなどのフォントもそうですし、CADなどグラフィックス分野の方が使う事が多いと思いますが、間欠的なデータから連続した線を描画したいときは実用上は多々あります。ラグランジュ補間、ニュートン補間などバリエーションは幾つもありますが、Rではスプライン補間をする事が容易です。
「R でグラフ中にプロットされてない値を予測してみる」を見ていたら簡単そうなので、試してみました。 $y=\frac{\sin(x)}{x}$ と言う簡単そうで積分するとトリッキーで大変なルベーグ積分不可の関数の極値をつないでみます。
まずは、 $y=\frac{\sin(x)}{x}$ を描いてみます。

次に、極値を求めます。

最後に、極大値と極小値をスプライン曲線で結びます。

意味が良く分からないものの、それらしいグラフになりました。ソースコードは以下になります。

# 関数を定義
f <- function(x){
  sin(x)/x
}

# fの一階微分
df <- function(x){
  (cos(x)*x-sin(x))/x^2
}

# fの二階微分
ddf <- function(x){
  (-sin(x)*x*2*x - (cos(x)*x-sin(x))*2*x)/x^4
}

# xの変分
dx <- 0.01
# xを生成（始点を僅かに正にすることでyが∞を回避）
x <- seq(1e-6, 30, dx)
# yを計算
y <- f(x)

getExtremum <- function(x){
  dy <- df(x)
  n <- length(dy)
  # 微分係数dy[i]*dy[i+1]がゼロ以下の場合は符号が逆転しており極値近傍と見なせる
  x[dy[-n]*dy[-1] <= 0]
}
extremum <- getExtremum(x)

# ddf(x[i])が負ならば、極大値
ulim <- c(extremum[ddf(extremum)<0])
# ddf(x[i])が正ならば、極小値
llim <- c(extremum[ddf(extremum)>0])

# スプライン関数を導出
sp_u <- smooth.spline(ulim, f(ulim))
sp_l <- smooth.spline(llim, f(llim))

# xとyのそれぞれ上下限
ylim <- c(-0.5, 1)
xlim <- c(min(x), max(x))

# 余白位置などを調整
par(mgp=c(1, 1, 0), mai=c(1,1,1,1), mar=c(4,3,3,1))

# 元データをプロット
curve(f(x), xlim=xlim, ylim=ylim, type="l", xlab="x", ylab=expression(sin(x)/x), axes=FALSE, main="スプライン関数のお試し")

# スプラインを描画／predict()関数で補間点を計算
par(new=T); plot(predict(sp_u, x), xlim=xlim, ylim=ylim, type="l", lty="dotted", xlab="", ylab="", axes=FALSE, col="red2")
par(new=T); plot(predict(sp_l, x), xlim=xlim, ylim=ylim, type="l", lty="dotted", xlab="", ylab="", axes=FALSE, col="red2")

# スプライン関数の元もプロットしておく
par(new=T); plot(ulim, f(ulim), xlim=xlim, ylim=ylim, type="p", lty="dotted", xlab="", ylab="", axes=FALSE, bg="gray")
par(new=T); plot(llim, f(llim), xlim=xlim, ylim=ylim, type="p", lty="dotted", xlab="", ylab="", axes=FALSE, bg="gray")

# 軸を描画
axis(1, at=seq(0, 30, 10), pos=ylim[1])
axis(2, at=c(seq(ylim[1], ylim[2], 0.25)), pos=0)

極大値と極小値を計算するほうが手間がかかっていますね。

2013-05-22

Rで偏微分方程式

R言語

速度的に実用性が低い感じもするのですが、Rで偏微分方程式の数値解を計算することもできます。微分方程式一般に言えることですが、せっせと計算して閉じた解を求めるのは苦労するので、数値解は実用的には重要です。
deSolveパッケージを使えば三元連立式ぐらいまでは容易で、使い方はサンプル・コードを見ればよいのですが、恐らく最も見慣れている偏微分方程式である波動方程式のサンプルを書いてみました。
deSolveの旧バージョンであるReacTranパッケージのリファレンスに例としては上がっているわけですが、パッケージの構成が変わってそのままでは動かないので、もしかしたら参考になるかも知れません。

1. 計算の概要

まずは波動方程式を眺めてみましょう。波動関数は $u(t,x)$ で、 $x$ を位置、 $t$ を時間とします。 $c$ は波動の速度。
$\frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2}$
偏微分方程式のままだと解けないので、常微分方程式に変形します。
$\frac{\partial u_1}{\partial t} = u_2$
$\frac{\partial u_2}{\partial t} = c^2 \frac{\partial^2 u_1}{\partial x^2}$
初期パラメーターを以下のようにします。力学的振動のときと異なって、ベクトルで大量の値がセットされることに注意してください。
$u_1(0,x)=exp{-0.05x^2}$
$u_2(0,x)=u(t,\infty)=u(t,-\infty)=0$
$c=1$
この式からそのままコーディングできる人は、かなり数値演算に慣れていると思います。

2. 実際の演算コード

以下が実際のコードになります。ポイントは恐らく、二階微分を行うところになると思います。数値微分のちょっとした知識が要ります。deSolveパッケージの癖で、複数の変数のベクトルが、一つのベクトルにまとめられているのも、注意しないと混乱するかも知れません。

# deSolveパッケージを読み出す
library("deSolve")
# 位置xの刻み幅は0.2
dx <- 0.2
# xの最小値と最大値
xlim <- c(-70, 70)
# -40から40まで位置xのベクトルを得る
x <- seq(xlim[1], xlim[2], dx)
# xの数を保存しておく
N <- length(x)
# 波の高さの初期値
uini <- exp(-0.05 * x^2)
# 波の高さの変化分の初期値
vini <- rep(0, N)
# パラメーターは同じベクターにしておく
yini <- c(uini, vini)
# 時間tのベクトルを得る
times <- seq (from = 0, to = 50, by = 1)
# 常微分方程式の形式で、波動方程式を表現する
wave <- function(t, y, parms, dx, N) {
# ベクターを二つに分ける
u1 <- y[1:N]
u2 <- y[(N+1):(2*N)]
# 波の高さの変化を計算
du1 <- u2
# 波の高さの変化分の変化の計算
# 二階微分を行う。両端の変化率は境界条件によりゼロ。
du2 <- c(0, (diff(u1[-1],1) - diff(u1[-N],1))/dx^2, 0)
# 変化を戻り値で返す
return(list(c(du1, du2)))
}
# ode.1D関数は一次元の微分方程式を扱う。
# nspecは計算する変数の数で、ここではu1とu2で二つ。
# 計算手法はデフォルトだと収束しなくなるので、ode45
out <- ode.1D(yini, times, wave, parms = 0, nspec = 2, dx=dx, N=N, method="ode45", atol=1e-14, rtol=1e-14)

ちょっと時間がかかると思いますが、気長に待てば計算できるでしょう。

3. 演算結果をプロットしてみる

出てきた計算結果は以下のようにプロットしたり、中身を見ることができます。

# 波の変化をプロットしてみる
matplot.1D(out, which=1, subset = time %in% seq(0, 50, by = 10), type = "l", col = c("black", rep("darkgrey", 5)), lwd = 2, grid = x, xlim = xlim, main="u(x, t)")
legend("topright", lty = 1:6, lwd = 2, col = c("black", rep("darkgrey", 5)), legend = paste("t = ",seq(0, 50, by = 10)))

# 等高線プロットをしてみる
image(out, xlab = "t", ylab = "x", main = "PDE", add.contour = FALSE)

左が $u(t, x)$ 、右が $du$ になります。色が波動の高さで、赤色が高く青色が低くなります。表題などを綺麗にする方法は調査中。

# t=30のときの大きさを波の形状を見てみる。ylim重要。
plot(out[30,][1:N], ylim=c(0, 1), xlab="x", ylab="u")
# 中央地点x=2*N/3の波の高さの変化を見る
plot(out[,2*N/3][1:length(times)], ylim=c(0, 1), xlab="t", ylab="u")

コンピューターを使っている気になれますね。

2013-05-07

Rで常微分方程式

R言語

自発的には揺れないボロビルの力学的振動は、単純化すると以下のように表されます。
$\frac{d^2y}{d^2t} + k\frac{dy}{dt} + \omega^2 y = 0$
$y$ が揺れ幅（変位）、 $\omega$ が固有振動数で、 $k$ が抵抗です。 $t$ は時間を表します。これをプロットしてみましょう。

1. 閉じた解を求める

定係数二階同次方程式なので、数学的に解くことができます。細かい計算は省略します。
$y = C_1 e^{-\frac{k}{2} t} \cos{ \Bigg ( t \sqrt{ \omega^2 - \frac{k^2}{4}} \Bigg )} + C_2 e^{-\frac{k}{2}t} \sin{ \Bigg ( t \sqrt{ \omega^2 - \frac{k^2}{4}} \Bigg ) }$
ほぼこのまま、curveで描画できます。

C1 <- 5; C2 <- 5; k <- 0.1; omega <- 1
curve(C1*exp(-k/2*t)*cos(t*sqrt(omega^2-k^2/4)) + C2*exp(-k/2*t)*sin(t*sqrt(omega^2-k^2/4)), 0, 100, n=1000, xname="t", ylab="y")

もっと複雑な式になると、辛い思いをします。

2. 微分方程式のまま処理する

$y=y_1$ 、 $dy_1/dt=y_2$ と置いて、連立微分方程式に展開してみましょう。
$\frac{dy_1}{dt} = y_2$
$\frac{dy_2}{dt} = - ky_2 -\omega^2 y_1$
この状態で、deSolveパッケージに渡すことができます。

library("deSolve")
model <- function(t, Y, params){
# params["omega"]としなくていいようにwithを使う
  with(as.list(params),{
    dy1 = Y[2]
    dy2 = -omega^2*Y[1] - k*Y[2]
    list(c(dy1, dy2))
  })
}
# 変位Y1とその微分の初期値
Y <- c(y1=8, y2=0)
# t=0からt=100まで0.1刻みで計算する
times <- seq(0, 100, 0.1)
# 常微分方程式を計算する
out1 <- ode(Y, times, model, parms = c(omega=1, k=0.1, F0=0.35, beta=1))
# プロットする上下の限界
ylim <- c(-10, 10)
# プロットする。軸は書かない
plot(y1 ~ time, data=out1, type="l", ylim=ylim, axes=FALSE, ylab="揺れ幅（＝変位）", xlab="時間", lwd=1, main="固有振動数と揺れ幅")
# 軸を描画
axis(1, labels=seq(0, 100, 25), at=seq(0, 100, 25), pos=ylim[1])
axis(2, labels=seq(-10, 10, 5), at=seq(-10, 10, 5), pos=0)

プロット結果は以下のようになります。

3. deSolveパッケージに関して

初期値問題、微分代数方程式、偏微分方程式も解ける優れモノです。詳しい説明は以下を参照してください。

Soetaert, Petzoldt, and Setzer (2010) "Solving Differential Equations in R: Package deSolve," Journal of Statistical Software, Vol.33(9)

2013-04-13

Rでパスワードを作る

R言語

最近はウェブサービスが乱立している事もあり、パスワードが漏洩する事件も良く聞くようになりました。この時に複数のサイトでパスワードが共通にしている人は、被害が他のサイトにも拡大していく可能性があります。
それを避けるために複数サイトで異なるパスワードを使うように推奨されていますが、人間がパスワードを生成すると何かの単語になりがちです。Appleのように厳しいパスワードを求めてくるサイトもあるので、自動で生成してみましょう。

#
# 'A', 'B', 'C' ... と言うパスワードに用いる文字のベクターを作る関数
# begin: 開始文字
# len: 長さ
# exception: 除外する文字
#
make_set <- function(begin, len, exception=c()){
  v <- strsplit(rawToChar(as.raw(0:(len-1) + as.integer(charToRaw(begin)))), "")[[1]]
  v[!v %in% exception]
}

#
# A〜Z、a〜z、0〜9のベクターを作る
# ただし紛らわしいIとl、Oと0は除外
#
chr_set_A <- make_set("A", 26, c("I", "O"))
chr_set_a <- make_set("a", 26, c("l"))
chr_set_0 <- make_set("0", 10, c("0"))
# 三つを合成する
chr_set <- c(chr_set_A, chr_set_a, chr_set_0)
# 乱数を初期化
# 注意：Seedの生成方法が分かるとパスワードを特定されやすくなるので、ある程度の長さの任意の数字（secret_code）を加えるなどした方が安全。
secret_code <- 192837465
set.seed(as.integer((1000*as.numeric(format(Sys.time(),"%y%j%H%M%OS3")) + secret_code) %% .Machine$integer.max))
#
# Apple IDを意識して、三種類の文字がないと脆弱パスワードと見なす関数
# clst: パスワードに使う文字のベクター
#
is_strong <- function(clst){
  if(!any(clst %in% chr_set_A)){
    return(FALSE)
  }
  if(!any(clst %in% chr_set_a)){
    return(FALSE)
  }
  if(!any(clst %in% chr_set_0)){
    return(FALSE)
  }
  # 三連続同一文字の禁止
  len <- length(clst)
  if(3 > len){
    return(FALSE)
  }
  for(i in 1:(len-2)){
    if(clst[i]==clst[i+1] & clst[i]==clst[i+2]){
      return(FALSE)
    }
  }
  TRUE
}

#
# 強度のあるパスワードを作る関数
# args: [パスワードの長さ#1],[パスワードの長さ#2]...
#
mkpasswd <- function(...){
  digit <- length(chr_set) # パスワードの各桁の文字種の数
  # Rらしく引数の数は不定
  len <- c(...)
  num <- length(len)
  # 戻り値を初期化しておく
  r <- character(num)
  for(i in 1:num){
    # 3以下のパラメーターは計算できないので、除外
    if(3>len[i]){
      r[i] <- NaN
      next
    }
    # 強度を満たすまで生成しなおす
    clst <- c()
    while(!is_strong(clst)){
      # +0.5がないと、最初の一文字の発生確率が半分になってしまう
      clst <- chr_set[round(runif(len[i], 0, digit) + 0.5)]
    }
    # collapse引数が無いとベクターは合体しない
    r[i] <- paste(clst, collapse="")
  }
  r
}

# 8文字のパスワードを作ってみる
sprintf("生成されたパスワード: %s", mkpasswd(8))

本当の問題は、どうやって大量のパスワードを管理するかなのですが、強度のあるパスワードを簡単に作る事ができます。

2013-01-15

Rでlmコマンドの結果からStd. Errorを取得する

R言語

たまに回帰分析の係数の区間推定を求めたいときがあります。
lmコマンドの推定結果から標準誤差（Std. Error）を取得するときは、以下のようにすると簡単です。

r <- lm(expression) # expressionは任意の推定式です
coef(summary(r))[, 2]

例えば3番目の係数の95%信頼区間を求めたいときは、以下のようにしてください。

a <- 0.05/2 # 両側検定
df <- summary(r)$df[2] # 自由度
b <- coef(r)[3] # 3番目の係数
se <- coef(summary(r))[, 2][3] # 3番目の標準誤差
sprintf("%.3f(95%%信頼区間%.3f〜%.3f)", b, b-se*qt(a, df), b+se*qt(a, df))

[3]を無くせば、全係数の区間推定量が出ます。

2012-03-18

Account Auto-DiscoveryをBloggerテンプレートに埋め込む

小ネタ

ブログ・サービスのBloggerはテンプレートが使えるのですが、ダブル・コーテーション（"）で書いたタグが、シングル・コーテーション（'）に勝手に変換される癖があります。仕様上は特に問題が無いのですが、はてなブックマークのAccount Auto-Discoveryを埋め込むときに苦労したのでメモを書きます。

Account Auto-Discoveryは、はてな以外のサービスにあるコンテンツに関して、はてなのブックマークなどのサービスの制御を許すプロトコルです。ブックマークされた後にタイトル等を変更したくなるときがある場合は、対応をしておくと便利です。

基本的にBloggerのテンプレートの〜の間に、以下のようなXMLをエスケープ内に埋め込めば良くなっています。

しかし以前から何度か試しているのですが、どうも上手く動きません。別のサイトの静的HTMLに同様のコードを埋めたときは上手く行くのにです。

おかしいなと思いつつ半年ぐらい経過したある日、別のサイトの静的HTMLをダブルコーテーションからシングルコーテーションに変えてみたら、うまく動かなくなりました。XHTMLの仕様上は、属性のブルコーテーションとシングルコーテーションに差異はないはずなのにです。

2月29日にはてなに要望を出したところ、3月15日に対応がされ、Account Auto-Discoveryが使えるようになりました。Bloggerユーザーで新テンプレートを使っていて、以前にAccount Auto-Discoveryの設定で悩んだ人は、今、リトライを行うと上手く行くかも知れません。

2012-02-18

Rで解くFizzBuzz問題

R言語

前のエントリーでFizzBuzz問題がもっと速くなるのでは無いかと指摘したら、“エレガントではないけどと速いコード”を書いて頂けたので、他のコードとあわせて紹介したいと思います。遅い順に合計7種類。

1. @kos59125 氏 Type I

@kos59125 氏が説明用に示したオーソドックスなコードです（Rで解くFizzBuzz問題）。C言語等で実装すると、それなりの速度が出るパターンですが、Rでは遅いです。

fb1<-function(n){
  fb<-character(n)
  for(i in 1:n){
    if(i%%3==0&&i%%5==0)
      fb[i]<-"FizzBuzz"
    else
      if(i%%3==0)
        fb[i]<-"Fizz"
    else
      if(i%%5==0)
        fb[i]<-"Buzz"
    else
      fb[i]<-i
  }
  fb
}

2. @kos59125 氏 Type II

@kos59125 氏が説明用に示したapply関数を使ってみたパターンのようです（Rで解くFizzBuzz問題）。(1)より高速ですが、もっと高速化が可能です。

fb2<-function(n){
  sapply(1:n,function(i){
    if(i%%3==0&&i%%5==0)
      "FizzBuzz"
    else if(i%%3==0)
      "Fizz"
    else if(i%%5==0)
      "Buzz"
    else
      i
  })
}

3. @kos59125 氏 Type III

高速版として紹介されていたコードです（Rで解くFizzBuzz問題）。ややトリッキーですが、rep()を使って剰余演算なしにしている所が改良に思えます。しかし、まだ高速化できます。

fb3<-function(n){
  isFizz<-rep(c(F,F,T),length=n)
  isBuzz<-rep(c(F,F,F,F,T),length=n)
  isNumber<-!isFizz&!isBuzz
  fizz<-ifelse(isFizz,"Fizz","")
  buzz<-ifelse(isBuzz,"Buzz","")
  number<-ifelse(isNumber,1:n,"")
  paste(fizz,buzz,number,sep="")
}

4. 裏 RjpWiki 短縮コード

ワンライナーになっていた、裏 RjpWikiの短縮版コードです（R らしく FizzBuzz）。ベクトル演算ですが、C言語のスカラー演算みたいにあっさりしていますね。

fb4 <- function(n){
i <- 1:n
ifelse(i %% 15 == 0, "FizzBuzz", ifelse(i %% 3 == 0, "Fizz", ifelse(i %% 5 == 0, "Buzz", i)))
}

5. uncorrelatedパターン

私が速度的に改良してみたコードです。(4)からifelse()を減らしてみたパターンが「RでFizzBuzz問題を高速化」だったのですが、まだifelse()が残っていたのでそれを削ります。

fb5 <- function(n){
  i <- 1:n
  l <- c(0, "Fizz", "Buzz", "FizzBuzz")
  ii <- (0==i%%5)*2 + (0==i%%3)*1 + 1
  r <- l[ii]; ns <- which(1==ii); r[ns] <- ns; r
}

6. 裏 RjpWiki 高速コード

速度的ならより高速にできると裏 RjpWikiの人が指摘してくれたコードです（R らしいかもしれないがエレガントではない FizzBuzz）。条件演算やループ処理が無いので高速です。またr[1:(l%/%3)*3]は、r[0==i%%3]とするより高速になるようです。

fb6 <- function(n){
  r <- 1:n
  l <- length(r)
  r[1:(l%/%3)*3] <- "Fizz"
  r[1:(l%/%5)*5] <- "Buzz"
  r[1:(l%/%15)*15] <- "FizzBuzz"
  r
}

なお1:(l%/%3)*3はseq(3,l,3)と書き直すと遅くなります。Rは速度的に実行してみるまで分からない所が多いですね。

7. 裏 RjpWiki 変態チック版

マニアックに改良がされていたので、追記します（もっと変態チックな FizzBuzz）。

fb7 <- function(limit){
  ans3 <- rep(c("d", "d", "Fizz", "d", "Buzz", "Fizz",
"d", "d", "Fizz", "Buzz", "d", "Fizz",
"d", "d", "FizzBuzz"), ceiling(limit/15))[1:limit]
  temp <- ans3=="d"
  ans3[temp]<- which(temp)
  return(ans3)
}

3と5と15の最小公倍数が15である事を利用しているのと、3と5の倍数で無い位置にdを配置し、which関数を使って後で置き換えているのが特色です。which関数の使いどころを初めて見た気がします。

8. ベンチマーク

100万回、ループしてみて速度を比較してみました。(6)がやはり高速です。添字のために大量の配列を作っているわけですが、ループ演算とifelse()のオーバーヘッドを考慮すると、微々たるコストしかないようです。

N <- 1000000
gc();gc();system.time({ fb1(N) })
gc();gc();system.time({ fb2(N) })
gc();gc();system.time({ fb3(N) })
gc();gc();system.time({ fb4(N) })
gc();gc();system.time({ fb5(N) })
gc();gc();system.time({ fb6(N) })
gc();gc();system.time({ fb7(N) })

fb1	fb2	fb3	fb4	fb5	fb6	fb7
11.08秒	10.36秒	4.39秒	5.52秒	0.87秒	0.79秒	0.58秒

餡子付゛録゛

ソフトウェア開発ツールの便利な使い方を紹介。