餡子付゛録゛

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

「年/月」の形式を月末日で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_t = a_i + \beta X_t + \sum^T_{i=1} \gamma_i TD_i + \epsilon
Yが被説明変数、Xが説明変数、aが固定効果、βとγは推定する係数、Iが個々の数、iが個々を表す添字、Tが期間の長さ、tが期を表す添字とします。εは誤差項です。

2. 係数の真の値

Iは10、Tは5として、真の値はβが0.015、\gamma_tが{0.00 0.05 0.04 0.06 0.02}、a_iが{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の関係を現す真のモデルを考えます。
z = 1 + 2 \cdot x + 3 \cdot y
「観測値は非確率変数で、未知の何かが誤差項を作っている」と仮定すると終了するので、xの観測値x'はxと正規分布に従う観測誤差、yの観測値y'はと正規分布に従う観測誤差の合算とします。後のシミュレーションでパラメーターを減らすために、xとyの観測誤差はi.i.dとします。
変形してみましょう。
\frac{z}{y} = \frac{1}{y} + 2 \cdot \frac{x}{y} + 3
単純に推定すると1/yの影響を忘れてかなり大きな誤差が出ます。x/yがモデルに入るときは、1/yの影響が無い事が理論的に演繹されている気もしますが、今回は特定化バイアスではなく誤差の分布の影響を比較するために残します。
推定モデルを整理しておきます。εとηは誤差項です。
割り算なし:z = \beta_0 + \beta_1 x' + \beta_2 y' + \epsilon
割り算あり:\frac{z}{y'} = \beta_2 + \beta_1 \frac{x'}{y'} + \beta_0 \frac{1}{y'} + \eta
一様分布からxとyを作り出しzを計算して、xとyに誤差項を加えてx'とy'を作ってから、x'、y'、zでデータセットを作成し、そこから\beta_1の推定量を計算していきます。真の値は2です。
さて、割り算なしと割り算ありの違いを見てみましょう。

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

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

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

以下のようなx、y、zの関係を現す真のモデルを考えます。
z = 1 + 2 \cdot \frac{x}{y}
xとyの観測値の誤差は正規分布だと考えるので、x/yの観測誤差の分布はコーシー分布になります。
推定モデルを整理しておきます。
割り算なし:z \cdot y' = \beta_0 y' + \beta_1 x' + \epsilon
割り算あり:z = \beta_0 + \beta_1 \frac{x'}{y'} + \eta

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

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を使いたくない時もあります。

ポアソン回帰で推定しているモノはλの式

某所の(1)ポアソン回帰モデルの説明が、(2)対数変換OLSと同じになっている気がします。違うものだと思うのですが、シミュレーションをして(1)と(2)の推定をして確認してみました。

1. モデル

ポアソン分布はパラメーター\lambdaで決定されるわけですが、\lambdaを説明変数で説明するモデルになります。n個のパラメータがあり、xを説明変数、\betaを係数として、以下のような式ですね。
\lambda = e^{\beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_n x_n}
P(y=k) = \frac{\lambda^k e^{-\lambda}}{k!}
被説明変数yの値が0以上の整数kのときの確率を、xは間接的に決定するわけですね。教科書的には最尤法を用いて求めることになるみたいですが、実用的にはリンク関数を用いて一般化線形回帰モデルで推定できるようです。

2. データ作成

まずはポアソン回帰モデル用のデータxとyを作成します。

set.seed(20130919)
x <- round(runif(100, max=3))
lambda <- exp(1.1 + 1*x)
y <- numeric(length(x))
for(i in 1:length(x)){
  y[i] <- rpois(1, lambda[i])
}

3. Rで推定してみる

glm関数を使うだけですけど、推定できます。

r_pois <- glm(y ~ x, family=poisson)
summary(r_pois)

4. 対数変換OLSと比較してみる

比較のために、対数変換OLSをかけてみましょう*1

r_ols <- lm(log(y) ~ x)
summary(r_ols)

当たり前ですが、残差平方和も異なります。

sum((y - fitted.values(r_pois))^2)
sum((y - exp(fitted.values(r_ols)))^2)

5. over-dispersionのテスト

tの切片項が0を棄却されたら、つまり有意に0以外ならばover-dispersionでポアソン分布でないことになります。

t <- fitted.values(r_pois)
z <- ((y - t)^2 - y)/(t*sqrt(2))
summary(lm(z ~ 0 + t))

人工的なデータなのですが、P値が0.323でした。

*1:実は対数化できるようにyはゼロ以上の数字になっています。