実行環境 Windows 7 + PostgreSQL 9.5.0にPL/v8をインストール PostGIS 2.2のインストールは9.5 RC1の時と同様（以上、1月23日） PNGファイルをインポートしてPostGISラスタのサンプルデータを作る インポートしたラスタを再びPNGファイルに出力して確認（以上、1月24日） PostGISラスタで画像処理、続けてPNGファイルに出力するテンプレート PostGISの関数のソースを見て、PL/v8への置き換えを考える（以上、1月25日） コールバック関数ST_Range4maのPL/v8版を作る PL/pgSQLの単純化バージョンも作る SQLでもコールバック関数を作る テストの全体（以上、1月26日） コールバック関数4つ × 周辺範囲5つ × 複数回テスト、を一つのクエリで テスト結果のTSVをPostgreSQLにインポート PL/v8の方が、元からあるPL/pgSQLの数倍速い（以上、1月28日）
Contents


テスト結果のPostgreSQLテーブルをRで読み込む

数日空いてしまいましたが、前回PostgreSQLのテーブルにしたテスト結果を今日はRに読み込みます。いつものPL/Rは、Windows + Postgres 9.5用のコンパイル済みファイルがまだ公式サイトにないので使えず、代わりにConEmuでRtermを起動。Rのバージョンは ↓ のとおりです。
as.matrix(R.Version())
               [,1]
platform       "x86_64-w64-mingw32"
arch           "x86_64"
os             "mingw32"
system         "x86_64, mingw32"
status         ""
major          "3"
minor          "2.2"
year           "2015"
month          "08"
day            "14"
svn rev        "69053"
language       "R"
version.string "R version 3.2.2 (2015-08-14)"
nickname       "Fire Safety"


↓ パッケージRPostgreSQLをロードして、変数conにコネクションを格納。結果の画面は省略します。
library(RPostgreSQL)
con = dbConnect(PostgreSQL(), host='localhost', port=5432,
    dbname='test_postgis22', user='postgres', password='****')
# to close
# dbDisconnect(con)

dbGetQueryでコネクションとSQLを渡せば、PostgreSQLへのクエリ結果が自動的にデータフレームになります。SQLはコメントや改行もそのままで良く、ヒアドキュメント風に書けるのですごく便利。文字列のクォートとバックスラッシュのエスケープには気を使いますが、今回は無関係。

↓SQLでピクセル範囲と関数別の所要時間平均値・エラーバー用の幅を出し、Rに読み込んだところ。エラーバー用の幅は、平均値と95%信頼区間の上下限の間隔です。Rで信頼区間を出してもいいけど、最初からデータに入ってる方が便利なのでSQLで簡易に（t値を直接入れて）計算。
(dat = dbGetQuery(con, "
    select pixel_range, func_pre, avg(sec),
        stddev(sec) / sqrt(count(*))
        * 2.776445 -- qt(0.975, 5-1) on R
        as err
    from result_range4ma
    group by 1, 2
    order by 1 desc, 2 desc
"))

   pixel_range func_pre       avg         err
1            4       st 92.361922 11.49552451
2            4      sql 28.644770  0.16963925
3            4   simple 71.117605  0.16328329
4            4     plv8 14.807546  0.07551800
5            3       st 57.345901  0.46703252
6            3      sql 20.567276  0.24028493
7            3   simple 46.195001  0.24517174
8            3     plv8 10.111938  0.09228771
9            2       st 32.750697  0.10027578
10           2      sql 14.402145  0.18311965
11           2   simple 25.933485  0.26790897
12           2     plv8  6.430331  0.07915714
13           1       st 16.810590  0.55494263
14           1      sql 10.171218  0.07117002
15           1   simple 12.832783  0.06684748
16           1     plv8  3.859446  0.02598734
17           0       st  7.497573  0.09384754
18           0      sql  7.531693  0.12204661
19           0   simple  5.419449  0.03518697
20           0     plv8  2.218324  0.02873016


今回のテストと外れ値、95%信頼区間について

今回、各ケース（関数×ピクセル範囲）につき5回しか試行してませんが、1ケースを除いて所要時間はごく近く、95%信頼区間も狭いです。上のデータのerr列（平均値と95%信頼区間の上下限の間隔）のとおり。ただ気になるのは、前回の最後に書いた「関数=ビルトインのst、ピクセル範囲=4」のケース。1回だけ目立って所要時間が多い時があり、その影響で95%信頼区間も桁違いに広いです（上のデータの1行目）。

その1回を外れ値扱いできるか検討したところ、確かに一般的な基準「四分位範囲の上下1.5倍」の外なのですが、これを用いると他のケースでも多くの巻き添え（外れ値）が出てしまう（四分位範囲が極端に狭いので）という困った結果に。より極端な基準として「四分位範囲の上下3倍」（参考: How to Calculate Outliers）を用いても、他ケースで巻き添えが起きるのは同じ。釈然としないけど「その1回」だけ除く根拠は見つかりませんでした。

もう一つ困ったのは、このように「ある1回だけ上に飛び出た」測定値群に対しt分布で信頼区間を出すと、上だけでなく下にも同じ幅で広がること。どう考えても実感と違う…。そこで、よりロバストな方法としてMAD（Median absolute deviation、中央値からの偏差の絶対値）による信頼区間を試すと、もっともらしい結果になりました。

以上をきちんと書くと長くなるので別の機会とし、とりあえず今回は「外れ値」扱いなし、t分布を仮定した95%信頼区間でやってみた結果を書きます。


全体を一つの棒グラフにまとめる（エラーバー付）

以下、Rだけの作業です。前々項で読み込んだデータから棒グラフ用のarray↓ を作成。第1次元（行）はPL/v8などの関数の略称、第2次元（列）はピクセル範囲、第3次元は平均値／エラーバー幅。普通は行・列のmatrixで済むところ、今回は一つ一つの棒（平均値）にエラーバーを付けるのでarray化しました。
(ary = array(c(dat$avg, dat$err), c(4,5,2),
    dimnames = list(
        unique(dat$func_pre),
        unique(dat$pixel_range),
        c('avg', 'err')
)))

, , avg
              4        3         2         1        0
st     92.36192 57.34590 32.750697 16.810590 7.497573
sql    28.64477 20.56728 14.402145 10.171218 7.531693
simple 71.11761 46.19500 25.933485 12.832783 5.419449
plv8   14.80755 10.11194  6.430331  3.859446 2.218324

, , err
                4          3          2          1          0
st     11.4955245 0.46703252 0.10027578 0.55494263 0.09384754
sql     0.1696393 0.24028493 0.18311965 0.07117002 0.12204661
simple  0.1632833 0.24517174 0.26790897 0.06684748 0.03518697
plv8    0.0755180 0.09228771 0.07915714 0.02598734 0.02873016


もしエラーバー不要なら、ピクセル範囲でグループ化した棒グラフは1行で書けます。↓ 
barplot(ary[,,1], beside=T, horiz=T)


これにエラーバーを足し、ラベルとか目盛りとか色々整えようとすると途端に面倒なのがRならでは^^;でも一回作れば別の機会に役立つので、今回は ↓ こんな風にしました。出力結果が、冒頭の画像です。
avg = ary[,,'avg']
cols = adjustcolor(terrain.colors(4), alpha=0.5)
xlim = c(0, max(apply(ary, 1:2, sum)))
yloc = barplot(avg, beside=T, plot=F)
ylim = range(yloc) + c(-1, 1) * 1.5

# set plot area
par(lend=1, mfrow=c(1,1), plt=c(12, 98, 15, 98)/100)
plot(axes=F, type='n', x=0, xlim=xlim, ylim=ylim, xaxs='i', yaxs='i',
    xlab='', ylab='')
axs = axis(1)
usr = par('usr')

# plot
barplot(add=T, axes=F, beside=T, col=cols, height=avg, horiz=T,
    las=1, names.arg=paste(dimnames(ary)[[2]], 'pixel\nrange'),
    xlim=xlim, ylim=ylim)
for (r in seq(nrow(ary))) {
    for (c in seq(ncol(ary))) {
        tmp = ary[r,c,]
        arrows(angle=90, code=3, length=0, lwd=5,
            x0=tmp[1] + tmp[2],
            x1=tmp[1] - tmp[2],
            y0=yloc[r, c])
    }
}
mtext(1, text='elapsed time (sec)', line=2.5)
segments(x0=axs, y0=usr[3], y1=usr[4], lty=c(1, rep(2, length(axs) - 1)))

# mask for legend
rect(mean(usr[1:2]), mean(usr[3:4]), usr[2], usr[4], border=NA, col='white') 

# legned
nr = nrow(ary)
leg = c(rev(dimnames(ary)[[1]]), NA, 'the 95% confidence interval')
leg = sub('^(simple)', '\\1 (modified pl/pgsql)', leg)
leg = sub('^(st)', '\\1 (built-in pl/pgsql)', leg)
legend(
    adj=c(1, 1/2),
    border=c(rep(1, nr+1), NA),
    bty='n',
    legend=leg,
    fill=c(rev(cols), NA, NA),
    lwd=c(rep(NA, nr+1), 5),
    x=xlim[1] + diff(xlim) * 0.85,
    x.intersp=c(rep(-3.5, nr), NA, -2.7),
    y=usr[4]
)


ピクセル範囲が狭い（所要時間が少ない）部分でエラーバーがほとんど見えないので、ピクセル範囲別に五つのグラフを作って並べると ↓ こんな感じ。凡例は省略しました。これでもやっぱり、信頼区間が狭いケースだとエラーバーがほとんどor全然見えません。
par(lend=1, mfrow=c(2,3), cex=1, plt=c(15, 96, 22, 99)/100)
for (i in rev(colnames(ary))) {
    mat = ary[,i,]
    bar = barplot(mat[,1], horiz=T, col=cols, xlim=c(0, max(rowSums(mat))))
    apply(cbind(mat, bar), 1, function(d) {
        arrows(angle=90, code=3, length=0, lwd=5,
            x0=d[1]-d[2], x1=d[1]+d[2], y0=d[3])
    })
    mtext(1, text=paste(i, 'pixel range, elapsed time (sec)'), line=2.5)
}


ピクセル範囲から所要時間を予測する回帰式（関数別）

絶対的な所要時間の傾向は分かったので、もう少し分析的に、関数別の回帰式（説明変数:ピクセル範囲、応答変数:所要時間）を作ってみます。もう棒グラフの時のような集計値のarrayは不要で、単純に全データをカテゴリの順序を揃えて読み込めばOK。↓
(dat2 = dbGetQuery(con, "
    select pixel_range, func_pre, sec
    from result_range4ma
    order by 1, 2
"))

    pixel_range func_pre        sec
1             0     plv8   2.215204
2             0     plv8   2.215203
3             0     plv8   2.246404
4             0     plv8   2.184004
5             0     plv8   2.230803
6             0   simple   5.444409
7             0   simple   5.397609
8             0   simple   5.428810
9             0   simple   5.382010
10            0   simple   5.444409
11            0      sql   7.425613
12            0      sql   7.690814
13            0      sql   7.488013
(...)


↓ まず線型（1次式）の回帰で、4つの関数別に実行しつつ結果を1画面にまとめてプロット。
par(lend=1, mfrow=c(2,2), cex=1, plt=c(18, 98, 20, 98)/100)
options(show.signif.stars=F)
ylm = c(0, max(dat2$sec))

for(f in unique(dat2$func_pre)) {
    d = subset(dat2, func_pre==f)
    attach(d, warn.conflicts=F)
    pxr = pixel_range
    lm = lm(sec ~ pxr)
    smry = summary(lm)
    str = paste(collapse='\n', c(
            capture.output(smry)[9:12],
            paste('r.squared', round(smry$r.squared, 4))))
    message('-----\n', f, '\n', str)

    plot(pxr, sec, xlab='pixel_range', ylim=ylm)
    abline(lm, col='red')
    text(adj=0, cex=2, lab=f, x=0, y=102.5)
    text(adj=0, cex=0.8, x=0, y=80, lab=str)
}
-----
plv8
Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   1.1993     0.3122   3.842 0.000833
pxr           3.1431     0.1275  24.660  < 2e-16
r.squared 0.9636
-----
simple
Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  -0.6520     1.8237  -0.358    0.724
pxr          16.4759     0.7445  22.129   <2e-16
r.squared 0.9551
-----
sql
Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   5.7390     0.5550   10.34 4.04e-10
pxr           5.2622     0.2266   23.23  < 2e-16
r.squared 0.9591
-----
st
Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  -0.6995     2.9250  -0.239    0.813
pxr          21.0264     1.1941  17.608 7.61e-15
r.squared 0.9309


回帰式の当てはまりは悪くないけど、プロットされた実測値を見ると非線型の方が妥当な気がします。方法が色々ある中、今回は簡単に説明変数（ピクセル範囲）を2乗して線型回帰に代入するだけ。ただし予測曲線のプロットで一手間増えます（ablineを使えないので、curveで描く）。↓
for(f in unique(dat2$func_pre)) {
    d = subset(dat2, func_pre==f)
    attach(d, warn.conflicts=F)
    pxr = pixel_range
    lm = lm(sec ~ I(pxr^2))
    smry = summary(lm)
    str = paste(collapse='\n', c(
            capture.output(smry)[9:12],
            paste('r.squared', round(smry$r.squared, 4))))
    message('-----\n', f, '\n', str)

    plot(pxr, sec, xlab='pixel_range', ylim=ylm)
    usr = par('usr')
    
    coef = smry$coefficients
    func = function(x) { coef[1,1] + x ^ 2 * coef[2,1] }
    
    curve(add=T, col='red', expr=func, from=usr[1], to=usr[2])
    text(adj=0, cex=2, lab=f, x=0, y=102.5)
    text(adj=0, cex=0.8, x=0, y=80, lab=str)
}
-----
plv8
Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  2.90126    0.12979   22.35   <2e-16
I(pxr^2)     0.76404    0.01543   49.53   <2e-16
r.squared 0.9907
-----
simple
Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  8.11989    0.53482   15.18 1.77e-13
I(pxr^2)     4.02996    0.06356   63.40  < 2e-16
r.squared 0.9943
-----
sql
Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  8.56346    0.19504   43.91   <2e-16
I(pxr^2)     1.28333    0.02318   55.37   <2e-16
r.squared 0.9926
-----
st
Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  10.2803     1.2068   8.519 1.44e-08
I(pxr^2)      5.1788     0.1434  36.109  < 2e-16
r.squared 0.9827


処理するピクセル数がピクセル範囲の2乗なので、この方が当てはまりが良くて当然。simpleとsqlを比べると、切片（ピクセル範囲ゼロでの所要時間）はsimple < sqlですが係数の大小が逆なので、範囲が増えるほどsimple > sqlになります。もし説明変数が連続値なら、回帰式どうしの交点から「優劣が逆転する境い目」を確定できます。


全体のまとめ

1月23日から続いた記事は、これで一区切り。まとめると

こんな感じかなぁ。最後の点はいずれ別の記事にします。