時系列データの季節調整をしてみよう

時系列データには、季節バイアスが入りがちです。年中行事はもちろんのこと、天候の変化も周期的に発生します。ゴールデンウィークのある5月よりも6月の方が行楽客が少ないとしても、行楽客が減少し出したとは言えないでしょう。そもそも月ごとに日数も異なりますし。月次データを見るときは、季節バイアスの影響を考慮する必要があります。これを数字の処理で行なうのが、季節調整です。実際にRで試してみて、どのぐらい調整できるか見てみましょう。単純移動平均とloessアルゴリズムを用います。

1. 精度を確認する手順

実データでは季節バイアスがどれぐらい入っているのか真実は誰も知らないので、季節バイアスが入った擬似データを作成します。二次方程式から作った非季節変動値に、乱数から作った真の季節バイアスを加えて、観測データを作ります。この観測データを処理して計算された季節バイアスを取り出し、真の季節バイアスと就き合わせてみましょう。

2. 擬似データを作成する

1996年1月から60ヶ月の月次データとして、12ヶ月周期の擬似データを作成します。

# 結果が同じになるように乱数シードを固定
set.seed(20150119)

# 12ヶ月5年分で観測数60
n <- 12*5

# 季節バイアスを乱数から生成
s_bias <- runif(12, min=0, max=4)
# 平均0に正規化しておく
s_bias <- s_bias - mean(s_bias)

# 季節バイアスを連結して、5年分の季節変化をつくる
# tsで時系列データ型にしておく
s_chg <- ts(rep(s_bias, 5), start=c(1996, 1), frequency=12)

# 二次方程式から非季節変動値を作る
x <- 1:n
e <- 0 # 誤差はゼロとしておく。入れるときは rnorm(n) などn個のベクトルを指定。
nsv <- ts(0.9 + 0.003*(x-2*n/3)^2 + 0.1*x, start=c(1996, 1), frequency=12)

# 非季節変動値と季節バイアスを足して、原数値を作成する
y <- ts(nsv + s_chg + e, start=c(1996, 1), frequency=12)

3. 単純移動平均

12ヶ月の算術平均をとる単純移動平均は、単純な計算の割りには効果的に季節バイアスを除去します。最初の11ヶ月はデータが無い欠損値になること、非季節変動値に変化があってから効果が見えるのが遅れことが欠点になりますが、目視するには悪く無いです。

# 移動平均
# 最初11ヶ月分の値はでない。
ma <- ts(numeric(n-11), start=c(1996, 12), frequency=12) # グラフを描くために、時系列データ型にしておく
for(i in 12:n){
  ma[i-11] <- mean(y[(i-11):i])
}

4. loessアルゴリズム

Rの組み込み関数stl()で使えるloessアルゴリズムは、局所重み付け回帰関数を用いた高機能な季節調整方法です。ここでは詳しい手順は考えずに、Rに計算させましょう。なお、引数で与える変数が時系列データ型でないと動いてくれないので、注意してください。

stl.y <- stl(y, s.window="periodic")

5. グラフを描いて比較してみる

どれぐらいの精度で季節変動を補正できるのか、グラフを描いて比較してみましょう。この例では、かなりの精度で一致します。移動平均も遅れて動くわけですが、十分機能することが分かります。

plot(nsv + e, ylim=c(0, 10), main="季節調整方法の精度", lwd=2, xlab="", ylab="")
lines(y, lty=2, lwd=1)
lines(y - stl.y$time.series[,1], lty=3, col="red", lwd=2)
lines(ma, lty=4, col="blue", lwd=2)

legend("topleft", c("真の非季節変動値", "原数値", "loessアルゴリズム", "12ヶ月移動平均"), col=c("black", "black", "red", "blue"), lty=c(1, 2, 3, 4), lwd=c(2, 1, 2, 2), bg='white', box.col = 'black', bty="n")

ただし、これは非季節変動値が滑らかな二次関数だからで、ゼロにしている誤差が大きくなると当てはまりが悪くなります。季節変動の標準偏差と同じ大きさの標準偏差を持つ誤差項を入れると、トレンドの変化を掴むのに問題があるほどでは無い*1ですが、真の値と計算値の乖離は目で見て分かる程度になります。

6. X-12-ARIMA

最近はX-13-ARIMAが出てきたようですが、官公庁などでは米国商務省が開発した移動平均型季節調整法であるX-12-ARIMAを使っています。これはアプリケーションが配布されていて、Rからも扱うためのx12パッケージもあります。
Windows機の場合、C:/WinX12にX-12-ARIMAを置いたら、Rでx12パッケージをインストールしたあと、以下のように使います。

library("x12")
x12path("C:/WinX12/x12a/x12a.exe")
x12.y <- x12(y)
# 季節調整値を表示
slot(x12.y, "d11")

オプションや同時計算される系列が多いので扱いが難しいのですが、季節バイアスの変化にも対応できるし精度は高いです。何十年分の統計を抱えている官公庁以外、使い道が限られそうですが。

*1:実際、回帰分析を行なっても、真の値でy=0.003*x^2 - 0.14*x、loessアルゴリズムによる推定値でy=0.002998*x^2 - 0.1402*xとなり、ほとんど差が無い。

「年/月」の形式を月末日でDate型に変更

月次データで「年/月」と言う形式（e.g. 2013/11、2014/2）は良く見かけると思いますが、Rのread.table関数でデータフレームに読み込むと文字列型*1になってしまいます。日付型でないと、subset関数などで絞り込むときに不便ですね*2。しかし、as.Date関数はこの形式を日付型に変えてくれません。そこで、ちょっとした細工をしてみましょう。
日経平均.txtと言うファイルを読み込みます。

df1 <- read.table("日経平均.txt", header=TRUE, sep="\t")

中身は、こんな感じです。

年月	始値	高値	安値	終値
1994/01	17369.74	20229.12	17369.74	20229.12
1994/02	20416.34	20416.34	18931.39	19997.20
1994/03	20216.62	20677.77	19111.92	19111.92

as.Dateをすると、こんな悲劇に。日付がないのがいけない模様です。

as.Date(df1$年月,"%Y/%m")
[1] NA NA NA NA NA NA NA NA NA

gsubで置換して日付を足せば変換できます。

df1$年月日 <- as.Date(gsub("([0-9]+)/([0-9]+)", "\\1/\\2/1", df1$年月))

これをsubset(df1, 年月日>="2013-1-1" & 年月日<="2013-12-31")のように絞り込めます。

日付を月末にする

少し応用して、2013-12-31のように月末値を代入しましょう。月末日は毎月かわりますし、うるう年の処理もいることに注意してください。
まずは月始日*3で、Date型ではなく、POSIXlt型の変数を作ります。

tmpDate <- as.POSIXlt(gsub("([0-9]+)/([0-9]+)", "\\1/\\2/1", df1$年月))

tmpDate$yearで年から1900を引いたものが、tmpDate$monで月から1を引いたものが、tmpDate$mdayで日が取得できます。これから翌月の月始日の日付を作成します。12月の翌月は、年が一つ増えて、1月に戻ることに注意してください。

nextMonth <- ISOdate(tmpDate$year+1900+(tmpDate$mon+1==12)*1, ((tmpDate$mon+1) %% 12)+1, 1)

他の環境に習熟している人は、mktimeではなくISOdateかと思うかも知れません。
翌月の月始日の24時間前（=3600*24秒）が当月の月末です。

df1$年月日 <- as.Date(nextMonth - 3600*24)

別解

見直して、コードの見通しを改善しました。

tmpDate <- as.POSIXlt(gsub("([0-9]+)/([0-9]+)", "\\1/\\2/1", df1$年月))
tmpDate$mon <- tmpDate$mon + 1 # 1ヶ月ずらす（桁上がり処理はされる）
tmpDate$mday <- 0 # 翌月の0日目は、当月の月末
df1$年月日 <- as.Date(tmpDate)

*1:正確には因子になります。

*2:文字列として2013/2と2013/11を比較すると、2013/11のほうが小さくなります。2013/02としてあれば良いのですが。

*3:必ず1日です。

Excelでもプログラミングができる事を忘れた経済学徒へ

Excel/VBAで大量のワークシートを一括してCSVで保存できます。
とりあえずコードだけ。Excelブック内にマクロを置く/書く方法は「Excel 2010 VBA の基礎知識」（Excel for Mac 2011）などを参照してください。

Sub J_YAMASAKI()
For Each mySheet In Worksheets
mySheet.Activate
' 保存先ディレクトリは C:\ を ~/J_YAMASAKI のように適当に変えてお使いください
ActiveWorkbook.SaveAs Filename:="C:\" & mySheet.Name & ".csv", FileFormat:=xlCSV, _
CreateBackup:=False
Next
End Sub

パネルデータで全員が一斉に説明変数の値を変化させたときのタイムダミーの係数

立命館大学の柴田悠氏への私信なのですが、パネルデータ分析の固定効果モデルに年次などタイムダミーが入っているときに、ある瞬間に個々が同じような行動をとって、説明変数の値を一斉に変化させたときに、それがタイムダミーの係数として観測されない事を、シミュレーションで示したいと思います*1。

1. モデル

計量モデルは以下になります。

Yが被説明変数、Xが説明変数、aが固定効果、βとγは推定する係数、Iが個々の数、iが個々を表す添字、Tが期間の長さ、tが期を表す添字とします。εは誤差項です。

2. 係数の真の値

Iは10、Tは5として、真の値はβが0.015、が{0.00 0.05 0.04 0.06 0.02}、が{0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0}を設定します。

3. 観測値の生成

観測値になるXは一様分布から乱数で0から1を生成し、t=3のときだけ+1するようにしましょう。ある瞬間に個々が同じような行動をとった例とします。誤差項εは標準偏差0.1の正規分布から生成します。Xとεが決まれば、真の値から観測値Yも定まります。
Rで生成するためのコードは以下になります。

set.seed(201408)
I <- 10 # 個体の数
ai <- (1:i)/10 # 固定効果
T <- 5 # 期間の長さ
gamma <- c(0, (5^(1:(T-1)) %% 7)/100) # タイムダミーの効果
epsilin <- rnorm(I*T, sd=0.1)
df <- data.frame(
  i = numeric(I*T), t = numeric(I*T), Y = numeric(I*T), X = numeric(I*T)
)
for(i in 1:I){
  for(t in 1:T){
    l <- (i-1)*T + t
    df$i[l] <- i
    df$t[l] <- t
    df$X[l] <- runif(1) + ifelse(t==3, 1, 0)
    df$Y[l] <- ai[m] + 0.15*df$X[l] + gamma[t] + epsilin[l]
  }
}

生成した観測値の一部を見ると、こうなります。t=3のときに、明らかにXが大きくなっていることを確認してください。

i	t	Y	X
1	1	1.0242424	0.90365560
1	2	1.2447202	0.74466442
1	3	1.2130497	1.03975109
1	4	1.3497088	0.29219997
1	5	1.0248503	0.02165759
2	1	1.3007814	0.63130932
2	2	1.3448006	0.86396715
2	3	1.2383986	1.17330402
2	4	1.2283682	0.66754810
2	5	1.1475519	0.75067122

4. 推定結果

ちゃっちゃと固定効果モデルで推定してみましょう。plmパッケージを使います。

library(plm)
plmdata <- plm.data(df, index=c("i", "t"))
td <- as.factor(df$t)
rfe <- plm(Y ~ X + td, model="within", data=plmdata)
summary(rfe)

推定結果は以下になります。

Oneway (individual) effect Within Model

Call:
plm(formula = Y ~ X + td, data = plmdata, model = "within")

Balanced Panel: n=10, T=5, N=50

Residuals :
Min. 1st Qu. Median 3rd Qu. Max.

• 0.16200 -0.05120 -0.00315 0.04550 0.16500

Coefficients :
Estimate Std. Error t-value Pr(>|t|)
X 0.188583 0.049956 3.7750 0.0005953 ***
td2 0.067277 0.041042 1.6392 0.1101244
td3 -0.015109 0.059223 -0.2551 0.8001284
td4 0.113756 0.040291 2.8233 0.0077890 **
td5 0.043061 0.040308 1.0683 0.2926920

• • -

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Total Sum of Squares: 0.54784
Residual Sum of Squares: 0.284
R-Squared : 0.4816
Adj. R-Squared : 0.33712
F-statistic: 6.50321 on 5 and 35 DF, p-value: 0.00022576

Xが0.188583と推定されている一方で、t=3のタイムダミーの係数は-0.015109となっています。観測数が十分ではないのか精度に欠けますが、タイムダミーの値は過少に推定されていますから、タイムダミーの係数として観測されない事が分かります。
また、t=3のインパクトを大きくとっても、推定結果が大きく変わることはありません。観測数を十分増やすと、推定値は真の値に近づきます。この例ではI=300ぐらいにすると、Xの係数βが0.1541469になりました*2。

*1:数理的には説明変数とタイムダミーの共分散がゼロである限りは、十分なサンプルサイズで推定結果に影響を及ぼさないことは、ほぼ自明です。

*2:係数の割りには分散が大きすぎたかも知れません。

統計モデルに観測値と観測値の割り算値を入れたときのバイアス

「データ解析のための統計モデリング入門」の著者の久保拓弥氏へのお返事なのですが、別所の「統計モデルに観測値と観測値の割り算値を入れても問題ない」に対して以下のようなコメントをもらったので、具体的な例を作ってみました*1。

…まあ，みなさん割算値だい好きだし，これからも割算値いりの統計モデルが使われるのでしょう．とりあえず「この場合は明らかにかなりまずいでしょ」とわかる例を何か作ってみようかな…

誤差項がコーシー分布のようなものになったときに、どのぐらいバイアスがかかるのかに関して、直観的な話です。

1. 真のモデルに割り算が無いケースを考える

以下のようなx、y、zの関係を現す真のモデルを考えます。

「観測値は非確率変数で、未知の何かが誤差項を作っている」と仮定すると終了するので、xの観測値x'はxと正規分布に従う観測誤差、yの観測値y'はと正規分布に従う観測誤差の合算とします。後のシミュレーションでパラメーターを減らすために、xとyの観測誤差はi.i.dとします。
変形してみましょう。

単純に推定すると1/yの影響を忘れてかなり大きな誤差が出ます。x/yがモデルに入るときは、1/yの影響が無い事が理論的に演繹されている気もしますが、今回は特定化バイアスではなく誤差の分布の影響を比較するために残します。
推定モデルを整理しておきます。εとηは誤差項です。
割り算なし：
割り算あり：
一様分布からxとyを作り出しzを計算して、xとyに誤差項を加えてx'とy'を作ってから、x'、y'、zでデータセットを作成し、そこからの推定量を計算していきます。真の値は2です。
さて、割り算なしと割り算ありの違いを見てみましょう。

観測値の誤差が小さい場合は差は目立ちませんが、誤差が大きくなると割り算ありモデルの推定結果が不安定になります。「この場合は明らかにかなりまずいでしょ」とわかる例に思えるかも知れませんが、実は説明力を失って、真の値から乖離している領域であることに注意してください。なお測定誤差が大きくなると、切片項以外の推定される係数はゼロに近づきます。

なお割り算ありモデルの残差にシャピロ-ウィルク検定をかけると、ほぼ正規分布にならない事が分かります。一方で、割り算なしモデルの残差はほとんど正規分布していました。無理に割り算を入れるのは、よした方が良いようです。

2. 真のモデルに割り算があるケースを考える

以下のようなx、y、zの関係を現す真のモデルを考えます。

xとyの観測値の誤差は正規分布だと考えるので、x/yの観測誤差の分布はコーシー分布になります。
推定モデルを整理しておきます。
割り算なし：
割り算あり：

これも「この場合は明らかにかなりまずいでしょ」とわかる例に思えるかも知れませんが、ある程度は信頼のおける推定値である青点線の左側では、割り算のありなしで大差は見られません。
なお、こちらの推定の残差にシャピロ-ウィルク検定をかけると、どちらもほとんど正規分布ではないことになります。割り算ありモデルは割り算で、割り算なしモデルは掛け算で分布が歪んでしまうようです。

2.1 本当にヤバいぐらい壊れるケース

後から気づいたので追記します。setup関数内のxとyの初期値を以下にして、第2節のグラフを描くと以下のように割り算で派手に壊れます。

x <- rnorm(n, sd=5) # xの真の値（0にかかる正規分布）
y <- rnorm(n, sd=5) # yの真の値（0にかかる正規分布）

何が起きたのかと思っていたのですが、偶発的に異常値がx/yにできるので、それに推定結果が引っ張られてしまったようです。以下がループ中で推定値が飛んだとき（setupの引数がsd=0.03）のヒストグラム。

(；ﾟдﾟ) -200以下と、+100以上！？
xとyの真の値の分布がrnorm(n, sd=5) + 30だったりすると問題が出ないので、常に問題になるわけではないですが、x/yの分布図は見ておいたほうが無難みたいです。

3. モデルの特定化の方が大事に思える

誤差の大きさや特定化の精度の程度によりますが、割り算による影響が大きいとは限りません。どこかの演習問題でx/yの誤差の分布は正規分布で近似できるとありましたが、ある程度はそういう性質があるようです。
少なくとも観測値が非確率的と見なせるほど精度が高ければ気にしなくて良いですし、質的選択モデルのような計量モデルの都合で割り算を排除できないときに、観測値と観測値の割り算をためらう理由も無い気がします。心配性の人は、残差の分布を確認しましょう。
割り算を排除したとしても、それによってモデルの特定化に失敗していたら、誤差は大きくなります。変数をかけたり*2、対数をとったりしても誤差の分布は歪むわけです。現実的に完璧な推定はあり得ませんから、力の入れどころは考えていく必要があると思います。

A. ソースコード

上のグラフ二枚のプロットに使ったコードは以下になります。

#
# n: 観測数, sd: 標準偏差, hg: 不均一分散(TRUE/FALSE)
# model: 真のモデルの数式表記
# return: x, y, zの観測値
#
setup <- function(n=100, sd=1, hg=FALSE, model=parse(text="1+2*x-3*y")){
  x <- runif(n, 1, 3) # xの真の値
  y <- runif(n, 1, 2) # yの真の値
  u <- rnorm(n, sd=sd) # yの観測誤差
  oy <- y + u # yの観測値
  e <- rnorm(n, sd=sd) # xの観測誤差
  # 不均一分散のときは上書き
  if(hg){
    for(i in 1:n){
      e[i] <- rnorm(1, sd=max(1e-6, sd*y[i])) # 標準偏差は0より僅かに正にする
    }
  }
  ox <- x + e # xの観測値
  # 被説明変数zは真の値に従う
  z <- eval(model)
  data.frame(x = ox, y = oy, z=z)
}


#
# シミュレーションで比較してプロット
# n: 観測数, min_sd: 標準偏差の最小値, max_sd: 標準偏差の最大値
# hg: TRUEならば不均一分散を入れる（未使用）
# model: 真のモデルのテキスト表記
# plegend: 凡例の位置
#
compare <- function(n=100, min_sd=1e-06, max_sd=1, hg=FALSE, model="1+2*x-3*y", plegend="topright"){
  len <- 100
  compared <- data.frame(sd=numeric(len), non_divided=numeric(len), divided=numeric(len), nd_sw_test=numeric(len), d_sw_test=numeric(len), t_test=integer(len))
  compared$sd <- seq(max(1e-06, min_sd), max_sd, length.out=len)
  for(i in 1:len){
    df <- setup(n=n, sd=compared$sd[i], hg=hg, model=parse(text=model))
    # ベタベタな分岐
    if(model=="1+2*x-3*y"){
      # 普通に推定
      r <- lm(z ~ x + y, data=df)
      # 割り算推定
      z_p_oy <- df$z/df$y
      x_p_oy <- df$x/df$y
      c_p_oy <- 1/df$y
      r_p_oy <- lm(z_p_oy ~ x_p_oy + c_p_oy)
    }else{
      # 割り算あり
      yz = df$z*df$y
      r <- lm(yz ~ y + x + 0, data=df)
      # 割り算なし
      x_p_oy <- df$x/df$y
      r_p_oy <- lm(df$z ~ x_p_oy)
    }
    # 切片項
    compared$non_divided[i] <- coef(r)[2]
    compared$divided[i] <- coef(r_p_oy)[2]
    # 正規性の検定結果
    compared$nd_sw_test[i] <- shapiro.test(residuals(r))$p.value
    compared$d_sw_test[i] <- shapiro.test(residuals(r_p_oy))$p.value
    # 95%信頼区間に2.0が含まれるか
    compared$t_test[i] <- ifelse(abs(2.0 - coef(r)[2])<coef(summary(r))[, 2][2]*qt(1-0.05/2, r$df.residual), 1, 0)
  }
#  ylim <- c(min(compared$non_divided, compared$divided), max(compared$non_divided, compared$divided))
  ylim <- c(0, 4)
  plot(non_divided ~ sd, data=compared, xlim=c(min_sd, max_sd), ylim=ylim, type="l", main=sprintf("観測データ同士の割り算の影響\n観測数:%d%s", n, ifelse(hg,"，不均一分散あり","")), xlab=sprintf("xとyの観測誤差の標準偏差%s", ifelse(hg, "(y=1のとき)", "")), ylab=parse(text = "beta[1]"))
  lines(compared$sd, compared$divided, lty=2)
  abline(h=2, col="red", lty=3)
  # ある程度の信頼性の推定値を出した最大の標準偏差
  rsd <- max(subset(compared, t_test==1)$sd)
  abline(v=rsd, col="blue", lty=4)
  text(rsd, 0, "推定値の95%信頼区間に2が入る限界線", col="blue", adj=c(0, 0))
  legend(plegend, lty=c(1, 2, 3), pch=-1, legend=c("割り算なし", "割り算あり", "真の値"), col=c("black","black","red"), bty="n")
  return(compared)
}


# 第1節のグラフ
set.seed(201408)
rd <- compare(max_sd=0.5, n=1000)


# 第2節のグラフ
set.seed(201408)
rd <- compare(max_sd=0.5, n=1000, model="1+2*(x/y)")

B. 不均一分散

こうすると不均一分散があるときに、割り算があるときのほうが、マシな推定量になる事が分かります。まぁ、クマが体長10cmから50mまであったりしないでしょうから、社会科学系のデータでないと気にする必要はないかも知れませんが。

set.seed(201408)
rd <- compare(max_sd=0.5, n=1000, hg=TRUE)

*1:ソースコードは公開するので、ミスがあったときは御勘弁を

*2:第2種の変形されたベッセル関数になるそうです。残念ながら、モーメント母関数などを持つかは確認していません。

Rの拡張でライフゲームを作ってみる

わざわざC言語で書くほどのものではなく、単に行列の受渡しを確認しただけですが、Rの中で使えるライフゲームを作ってみたので公開します。
以下のソースコードを LifeGame.c と名づけて保存します。

#include <R.h>
#include <Rinternals.h>

SEXP LifeGame(SEXP m)
{
  SEXP  ans, dim;
  int  nor, noc, i, j, nol, x, y, r, c;
/* 第一引数が行列か確認 */
  if(!isMatrix(m)){
    error("A matrix is required for the first argument.");
  }
/* 行列の行数、列数を得る */
  dim = getAttrib(m, R_DimSymbol);
  nor = INTEGER(dim)[0]; /* 行数 */
  noc = INTEGER(dim)[1]; /* 列数 */
/* 戻り値の行列を初期化（注意：REALSXPではなくINTSXPを指定） */
  PROTECT(ans = allocMatrix(INTSXP, nor, noc));
/* 端は消去 */
  for(i=0; i<nor; i++){
    INTEGER(ans)[i] = 0;
    INTEGER(ans)[i + nor*(noc-1)] = 0;
  }
  for(j=0; j<noc; j++){
    INTEGER(ans)[nor*j] = 0;
    INTEGER(ans)[nor-1 + nor*j] = 0;
  }
/* 中央部分だけを処理 */
  for(i=1; i<nor-1; i++){
    for(j=1; j<noc-1; j++){
/* 周辺の8マスの状態を確認 */
      nol = 0;
      for(y=-1; y<=1; y++){
        for(x=-1; x<=1; x++){
          if(!y && !x)
            continue;
          if(0<INTEGER(m)[i+y + nor*(j+x)])
            nol++;
        }
      }
/* 周囲のマスの状態で生死を決定 */
      if(0<INTEGER(m)[i + nor*j])
        INTEGER(ans)[i + nor*j] = nol<=1 ? 0 : (nol>=4 ? 0 : 1);
/* 周囲のマスの状態で発生を決定 */
      else
        INTEGER(ans)[i + nor*j] = nol==3 ? 1 : 0;
    }
  }
  UNPROTECT(1);
  return(ans);
}

コンパイルしてDLL化します。

R CMD SHLIB LifeGame.c

Rを起動して、例えばグライダーを飛ばして遊びます。行列内の1と0で識別しているので、分かりづらいですね。

# ライブラリを呼び出す
dyn.load(paste("LifeGame", .Platform$dynlib.ext, sep = ""))
# Rの関数でラッピングして使う
lgame <- function(m){
  .Call("LifeGame", m)
}

# グライダー
mg <-matrix(as.integer(c(
1,1,1,
1,0,0,
0,1,0)), 3, 3, byrow=TRUE)
# 移動物体なので広めの行列を作る
m <- matrix(as.integer(rep(0,12*12)),12,12)
# 広い行列にグライダーを写す
for(r in 1:nrow(mg)){
  for(c in 1:ncol(mg)){
    m[r+nrow(m)-nrow(mg)-1,c+ncol(m)-ncol(mg)-1] <- mg[r,c]
  }
}
# 移動させてみる
for(i in 1:15){
  print(m)
  m <- lgame(m)
}

インストーラー作成ツールNSISをいじってみる

MS-Windowsのアプリケーションのためにインストーラーを作るときは、Microsoft Windows InstallerかInstallShieldが定番になっていると思いますが、Qtなどオープンソースで開発しているときはライセンス料の問題などで使いづらいかも知れません*1。しかし、そういう時は、フリーソフトのNSIS (Nullsoft Scriptable Install System)があるので、これを使えばいいようです。スクリプト・ファイルを使わないといけないので取っ付きづらい感じもあるのですが、試しに使ってみました。

1. 準備する

まずはNSISをこちらからダウンロードしてきます。数年前にダウンロードしてあったバージョン2.46と、最新版の3.0βで試しましたが、大した事をしなければ差異は無いようです。
次にWindows Applicationを用意します。ただのテキスト・ファイルでも問題ないのですが、昔書いた巡回セールスマン問題解法プログラムを引っ張り出してきます。

2. NSISをインストールする

ダウンロードしてきたnsis-2.46-setup.exeを実行します。インストーラーのメッセージは英語ですが、ライセンス条項に I agree をして、デフォルト構成でインストールして問題ないと思います。コマンド・プロンプトで makensis コマンドを使う事になりますが、C:\Program Files (x86)\NSIS などのインストール先にパスを通してくれないので注意しましょう。

3. スクリプトファイルを記述する

巡回セールスマン問題解法プログラムは、本体とドキュメントとQt関連のDLLのみで構成されたシンプルなものです。今回は、本体とドキュメントを%PROGRAMFILES%\TSP以下に置き、Qt関連のDLLをWindowsのシステム・ディレクトリに置くようにインストーラを書いてみました。

; ライブラリ操作を行うマクロを読み込む
!include "Library.nsh"

Name "TSP Installer" ; インストーラーの名前
OutFile "TSPInst.exe" ; インストーラーのファイル名

; インストール先ディレクトリ
InstallDir "$PROGRAMFILES\TSP"

; 必要な実行権限（user/admin）
RequestExecutionLevel admin

; インストール時に表示するものを指定
Page license ; ライセンスを表示（他のページより先に書く必要あり）
Page directory ; インストール先ディレクトリを表示
Page instfiles ; インストールしたファイル名を表示

; アンインストール時に表示するものを指定
UninstPage uninstConfirm ; アンインストールの確認
UninstPage instfiles ; アンインストールしたファイル名を表示

; 日英で国際化する。
LoadLanguageFile "${NSISDIR}\Contrib\Language files\English.nlf"
LoadLanguageFile "${NSISDIR}\Contrib\Language files\Japanese.nlf"

; ライセンス文書は日英別（rtfも使える）
LicenseLangString license ${LANG_ENGLISH} license-english.txt
LicenseLangString license ${LANG_JAPANESE} license-japanese.txt
LicenseData $(license)

; 日英でディレクトリ名を変える
LangString ShortcutDirectory ${LANG_ENGLISH} "Travelling Salesman Problem"
LangString ShortcutDirectory ${LANG_JAPANESE} "巡回セールスマン問題解法プログラム"

; インストール時の挙動
Section "Install" ; この名前は重要ではない

; インストール先ディレクトリを指定
SetOutPath $INSTDIR

; 変数定義
!define ReleaseDirectory 巡回セールスマン問題解法プログラム
!define RegistryKey "Software\Microsoft\Windows\CurrentVersion\Uninstall\TSP"

; インストールするファイルを指定
File ..\${ReleaseDirectory}\TSP.exe
File ..\${ReleaseDirectory}\readme.txt
File ..\${ReleaseDirectory}\TSP_ja_JP.qm

; 上書きインストールの場合は、DLLは展開しない
Var /GLOBAL ALREADY_INSTALLED
IfFileExists "$INSTDIR\TSP.exe" 0 new_installation
StrCpy $ALREADY_INSTALLED 1
new_installation:
; DLLをシステム・ディレクトリに展開
!insertmacro InstallLib DLL $ALREADY_INSTALLED REBOOT_NOTPROTECTED ..\${ReleaseDirectory}\QtCore4.dll $SYSDIR\QtCore4.dll $SYSDIR
!insertmacro InstallLib DLL $ALREADY_INSTALLED REBOOT_NOTPROTECTED ..\${ReleaseDirectory}\QtGui4.dll $SYSDIR\QtGui4.dll $SYSDIR
!insertmacro InstallLib DLL $ALREADY_INSTALLED REBOOT_NOTPROTECTED ..\${ReleaseDirectory}\libgcc_s_dw2-1.dll $SYSDIR\libgcc_s_dw2-1.dll $SYSDIR
!insertmacro InstallLib DLL $ALREADY_INSTALLED REBOOT_NOTPROTECTED ..\${ReleaseDirectory}\mingwm10.dll $SYSDIR\mingwm10.dll $SYSDIR

; アンインストーラーを書き出す
WriteUninstaller $INSTDIR\uninstaller.exe

; スタートメニューのショートカットを作る
CreateDirectory "$SMPROGRAMS\$(ShortcutDirectory)"
CreateShortCut "$SMPROGRAMS\$(ShortcutDirectory)\TSP.lnk" $INSTDIR\TSP.exe
CreateShortCut "$SMPROGRAMS\$(ShortcutDirectory)\uninstaller.lnk" $INSTDIR\uninstaller.exe

; レジストリに追加して「プログラムと機能」or「プログラム」に登録する
WriteRegStr HKLM ${RegistryKey} "DisplayName" "TSP:Travelling Salesman Problem"
WriteRegStr HKLM ${RegistryKey} "UninstallString" '"$INSTDIR\uninstaller.exe"'
WriteRegStr HKLM ${RegistryKey} "Publisher" 'uncorrelated@yahoo.co.jp'
WriteRegStr HKLM ${RegistryKey} "DisplayVersion" '1.0'

SectionEnd ; end the section

; アンインストール時の挙動
Section "Uninstall"

; ファイルを消す
Delete $INSTDIR\TSP.exe
Delete $INSTDIR\readme.txt
Delete $INSTDIR\TSP_ja_JP.qm

; システム・ディレクトリからDLLを消す
!insertmacro UnInstallLib DLL SHARED REBOOT_NOTPROTECTED $SYSDIR\QtCore4.dll
!insertmacro UnInstallLib DLL SHARED REBOOT_NOTPROTECTED $SYSDIR\QtGui4.dll
!insertmacro UnInstallLib DLL SHARED REBOOT_NOTPROTECTED $SYSDIR\libgcc_s_dw2-1.dll
!insertmacro UnInstallLib DLL SHARED REBOOT_NOTPROTECTED $SYSDIR\mingwm10.dll

; アンインストーラーも消去
Delete $INSTDIR\uninstaller.exe

; ディレクトリ消去。なお/rをつけると再帰的に、/REBOOTOKをつけると再起動時に消す。
RMDir $INSTDIR

; スタートメニューのショートカットを消す
Delete "$SMPROGRAMS\$(ShortcutDirectory)\TSP.lnk"
Delete "$SMPROGRAMS\$(ShortcutDirectory)\uninstaller.lnk"
RMDir "$SMPROGRAMS\$(ShortcutDirectory)"

; レジストリを消しておく
DeleteRegKey HKLM "Software\Microsoft\Windows\CurrentVersion\Uninstall\TSP"

SectionEnd

TSPInstaller.nsiと言うファイル名で保存しました。

4. 日英のライセンス・ファイルを記述する

国際化しているので、license-english.txtとlicense-japanese.txtを用意します。動かすだけなら、空ファイルでも構いません。

5. コンパイルする

以下のようにコンパイルします。

makensis TSPInstaller.nsi

これでTSPInst.exeが生成されます。

6. TSPInst.exeを動かす

動かすと以下のような今風ではない殺風景なインストーラーが起動します。

7. インストール状態を確認する

%PROGRAMFILES%\TSP以下にファイルがちゃんと展開されています。

DLLはC:\Windows\SysWOW64以下に配置されました。

これでスタート・メニューから起動とアンインストールができます。

また、「プログラムと機能」からアンインストールが可能になって、Windowsアプリケーションらしく見えます。

8. まとめ

小さなアプリケーションだと、いちいちインストールするよりもZIPファイルを展開したら使えるようにして欲しいと言われることもあるのでインストーラーを用意する必要性は無い事も多いのですが、たまにインストーラーを要望される事もあります。Windows環境でのNSISはWinampの他、OpenOfficeなどでも実績があり、国際化に対応していたり、外部実行ファイルを呼び出せたりと何かと高機能なので、そういう時は有力な選択肢の一つになると思います。
付属ドキュメントが英語で、ドキュメント中の例が部分的なところが扱いづらいと感じるかも知れませんが、付属のサンプル・ソースコードやオンライン・ドキュメントを検索して完全版のソースコードを見ていると分かることも多いので、オープンソースのツール等でWindowsアプリケーションを作っていてインストーラーを書く必要に迫られている人は、ぜひトライしてみてください。

*1:オープンソース／フリーソフトのWiXを使う選択肢もありますが、宗教的理由などで.NETを使いたくない時もあります。