先日の技術書典６で「お台場計算尺」というサークル名でいくつか本を頒布した。その際、10分ごとに売れた部数を記録した。そのデータを時系列解析してみる。

まずはデータ。「お台場計算尺」は何百部も売れるようなすごいサークルではない（やっとこ100部程度）。「12:00」は 12:00 までの10分間に売れた数だが、開始直前に見に来てくれた人（サークル参加の人）が結局買ってくれたので、11:00 の開始時刻までに1冊売れたことにした。

total.csv

11:00,1
11:10,0
11:20,0
11:30,0
11:40,2
11:50,4
12:00,14
12:10,5
12:20,6
12:30,4
12:40,2
12:50,4
13:00,2
13:10,5
13:20,4
13:30,8
13:40,6
13:50,4
14:00,2
14:10,1
14:20,1
14:30,3
14:40,0
14:50,8
15:00,0
15:10,5
15:20,1
15:30,5
15:40,1
15:50,1
16:00,3
16:10,2
16:20,2
16:30,0
16:40,0
16:50,2
17:00,0

データの様子

で、Ｒに読み込んで解析。まずはいつもの呪文。

> library(tidyverse);    # 呪文だから唱えてるけど、今回は不要
> library(ggfortify);    # 同上、snippet で自動的に書き込む
> library(latticeExtra); # 移動平均 simpleSmoothTs を使う
> theme_set(theme_bw(base_family = "HiraKakuProN-W3")); # plot で日本語

> dat <- read.csv("total.csv", header=F);
> dat[,1] <- (1:length(dat[,1]) - 1)/6; # hour 単位
> colnames(dat) <- c("hour", "売り上げ");

読み込んだら何はなくともまずプロット。書典では tidyverse をバリバリ使う本を出したが、正直めんどうなので、古式ゆかしく楽ちんプロット。12時ちょうど（開会後1時間）のところに外れ値っぽいのがある。全体的には、上がってから下がっていく感じ（山が一つ）っぽく見えるような気がする。

> total <- dat[,2];
> time  <- dat[,1];
> plot(time, total, ylab="10分ごとの売上げ数");

で、累積（開始から各時間までに売れた総数）でもプロット。外れ値のおかげか、全体の売上げの半分くらいは入場が有料だった時間帯（< 2h）に売れている。11:00〜17:00 の最初の1/3である。横線が売り上げ数の1/2のライン。

plot(time, cumsum(total), ylab="売上げ累積");
lines(c(0,6), c(sum(total)/2, sum(total)/2));

10分ごとの売り上げから、細かい変動を除いてみる（移動平均）。ウィンドウを広くしていくと、プロファイルは上がって下がる感じ（山が一つ）のような気がしないでもない。まぁ先入観かもしれない。なお、始めも終わりも裾を引かない。となると、変動に当てはめるのは二次曲線とかかなぁ、とこの段階では想像する。

x <- simpleSmoothTs(total, width= 5); plot(time, total); lines(time, x);
x <- simpleSmoothTs(total, width= 7); plot(time, total); lines(time, x);
x <- simpleSmoothTs(total, width= 11); plot(time, total); lines(time, x);
x <- simpleSmoothTs(total, width= 15); plot(time, total); lines(time, x);

直線や曲線の当てはめ

まずは手探り

ここまでも手探りだったけど、それじゃぁモデリングしてみましょうか、ということで、まずはとりあえず直線を当てはめてみる。まぁ、はまってるような、はまってないような...

mod <- lm(total ~ time);
plot(time, total);lines(time, predict(mod));

次に二次曲線を当てはめてみる。二次曲線は非線形だから nls で当てはめる。頂点の場所が分かりやすいように、曲線の式は $y = a(x - b)^2 + c$ にする。上に凸なので $a < 0$、$b$ はピークの時刻、$c$ はピークの高さになる。プロットから $b = 2.5$、$c = 6$ くらいになるかな、などと予想しながらやってみる。すると、まぁ直線より少しはいいような...でもまだ、なんとも言えないかなぁ。

mod <- nls(total ~ a*(time - b)^2 + c, start = c(a = -1, b = 2, c = 0));
plot(time, total);lines(time, predict(mod));

この二次曲線モデル自体を見てみると、ｐ値も小さいし、収束も早い。ピークは 2.46h あたりで、10分間に4.1冊の割合というのが売上げ速度の（期待値の）最大値のようだ。

> summary(mod);

Formula: total ~ a * (time - b)^2 + c

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
a  -0.3548     0.1602  -2.215   0.0335 *  
b   2.4608     0.4336   5.675 2.28e-06 ***
c   4.1456     0.6606   6.276 3.78e-07 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.756 on 34 degrees of freedom

Number of iterations to convergence: 3 
Achieved convergence tolerance: 9.315e-10

もっと高次だったらどうかな？ということで四次曲線を当てはめてみる。すると、当てはまりは悪くないし、また「上がる > プラトー > 下がる」というプロファイルはまぁ、売っていた時の感覚からも納得できる。最後の方はヒマだったし。

> mod <- nls(total ~ a*(time - b)^4 + c, start = c(a = -0.5, b = 2, c = 0));
> plot(time, total);lines(time, predict(mod));
> summary(mod);
Formula: total ~ a * (time - b)^4 + c

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
a -0.04681    0.01944  -2.408   0.0216 *  
b  2.84146    0.23405  12.141 6.52e-14 ***
c  3.78541    0.56799   6.665 1.20e-07 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.768 on 34 degrees of freedom

Number of iterations to convergence: 5 
Achieved convergence tolerance: 9.558e-08

最適モデル

二次より四次がいいなら、もっと高いといいのかもしれない。次数には最適値があるのかもしれない。というわけで探してみる。AIC を使って探すと六次式が最適（AIC 極小）になる。プロットを見ると、2h の前にピークがあってその後プラトー、そして最後に下がる。より感覚的に納得できる感じ。ｐ値は今一つだが、まぁこの辺りがとりあえずのゴールかな、と思う。もっと売り上げ数の多いサークルのデータで見てみたいところ。

> t <- time; # 長い変数名がイヤだから
> mod <- lm(total ~ 1 + t + I(t^2) + I(t^4) + I(t^6) + I(t^8) + I(t^10));
> mopt <- step(mod);
> plot(time, total);lines(time, predict(mopt));
> summary(mopt);

Call:
lm(formula = total ~ t + I(t^2) + I(t^4) + I(t^6))

Residuals:
    Min      1Q  Median      3Q     Max 
-3.1733 -1.4567 -0.5842  0.8173  9.5552 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept) -0.4237932  1.7057428  -0.248   0.8054  
t            7.4616040  2.9388969   2.539   0.0162 *
I(t^2)      -2.6696918  1.1822114  -2.258   0.0309 *
I(t^4)       0.0777012  0.0448673   1.732   0.0929 .
I(t^6)      -0.0010574  0.0007093  -1.491   0.1458  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.666 on 32 degrees of freedom
Multiple R-squared:  0.2726,	Adjusted R-squared:  0.1817 
F-statistic: 2.998 on 4 and 32 DF,  p-value: 0.03292

有名な時系列解析法ではどうか

まぁもういいかなと思いつつも、一応時系列解析的なこともやってみる。自己回帰とFFT。自己回帰は最適モデルの次数が低く、元データに見るべき周期性はないことを示している。FFT も周期成分（左端のパワー = sum(total) 以外の成分）には目立ったものはない。

> ar(total, order=16)

Call:
ar(x = total, order.max = 16)

Coefficients:
      1        2        3  
 0.2367   0.3650  -0.2760  

Order selected 3  sigma^2 estimated as  7.463
> plot(0:36, abs(fft(total)), type="h");

結論

技術書典６における、６時間の開催時間内の売り上げの推移は、多項式でモデリングするとスタートダッシュとプラトーのバイモーダルな感じになる。直線的あるいは上がって下がるだけ（山が一つ）と言った感じではない。また時間内には意味のある周期性（振動）はなさそう。

もっとマシな解析のためには

時間の単位が「分」なのは、それが人にとって感覚的に理解／把握しやすい長さだからだろう。その細かさのデータがあれば、もうちょっと分かりやすい解析ができるように思う。具体的には、かんたん後払いで400件（毎分1件）以上くらい売れていればよかったのかな、と思う（そのくらいあれば、10分間あたりの平均が10件以上になって、分布も見えたかもしれない）。今回は29件だった。理想は遠い。計算尺本は売れないのか（売るけど）。

なお、手でとった今回の記録をあとで「かんたん後払い」のデータと照らし合わせたところ、書き漏らしがあることが分かった。悲しい。昔はこんなミスはしなかった（と今は思うんだが）。年はとりたくないものだ...