Macでスクリーンショットの保存先をショートカットキーで振り分ける

ツイッターで尋ねられて、自分でも気になったので設定してみました。
Macでスクリーンショットを撮るときに、ショートカットキーを工夫して、保存先を振り分ける設定です。
Macのスクリーンショットは、じつはいろいろオプションがあって、ウィンドウ単位で撮影するときに影をつけるかどうかとかも変えられるのですが、私は基本的には、

shift+command+3（全画面を撮影）
shift+command+4（インタラクティブモードで範囲を指定して撮影。space）

の2パターンしか使いません。
で、仮に用途に応じて2つのフォルダに振り分けたいとすれば、計4つのショートカットがあればいいことになります。振り分けるフォルダは、仮に、
"~Desktop/ScreenShot1"
"~Desktop/ScreenShot2"
だとしておきます。

どうやってショートカットをつくるかというと、

1. 処理自体は、シェルスクリプトで書く
2. Automatorで、そのシェルスクリプトを呼び出すService（Auromator上ではQuick Actionと呼ばれるのでややこしいが）を作っておく
3. System Preferences（日本語環境の場合は「環境設定」） > Keybord > Shortcuts > Servicesで、さっきつくったサービスを呼び出すショートカットを設定する*1

という手順になります。
サービスってのは、たとえば↓の画像の場合だとFinderの歯車マークの中に並んでいる、様々な処理のことで、私は正確に理解してないですが「Mac上のいろいろな場面で、ワンクリックで呼び出せる処理」ぐらいに思っています。

さて設定手順ですが、まずAutomatorで、New DocumentをQuick Actionとして立ち上げます。

左のメニューの、Utilitiesから、Run Shell Scriptを、右のワークフローのところにドラッグします。その上にある設定欄は、一行目のinputを設定するところを、no inputにしときます。

次にスクリーンショットを撮って名前を付けて保存する処理をbashで書きますが（上の画像ではもう書いてありますが）、スクリーンショットを撮るコマンドである"screencapture"の細かいオプション設定は下記を参考にしてください。
https://do-zan.com/mac-terminal-screencapture/

ここでは、4パターンのうち、インタラクティブモードで撮影し、ファイルをデスクトップの"ScreenShot1"というフォルダに保存するパターンだけやります。ファイル名は、頭に"sc"をつけて、その後に"20200710140230"（2020年07月10日14時02分30秒）みたいに日付と時刻を並べて、拡張子".png"をつけることにします。
1行目は時刻を取得して格納する処理です。
"-i"ってのが、インタラクティブモードを表すオプション設定。
その後ろは、要するに保存するファイルのパスですね。${current_time}でさっき取得した時刻をくっつけています。

current_time=`date +"%Y%m%d%H%M%S"`
screencapture -i ~/Desktop/ScreenShot1/sc${current_time}.png

これをAutomatorのスクリプトの欄に書いて、たとえば"ScreenShot(Interactive-ScreenShot1)"という名前で保存します。

次に、System Preferences（日本語環境の場合は「システム環境設定」） > Keybord > Shortcuts > Servicesに進むと、一覧の下のほうに、さっきつくった"ScreenShot(Interactive-ScreenShot1)"がありますので、

その右のAdd Shortcutっていうボタンをクリックして、例えば
shift+command+7
と押します。すると、下記のように、ショートカットが登録されます。
（なお、shift+command+6のショートカットは、タッチバーを撮影する機能に割り当てられてて使えませんでした。）

さっきつくったサービス（クイックアクション）を修正したい場合は、
~/Library/Services
という場所にファイルが保存されてるので、それをAutomatorで開いて編集すればいいです。（さっきショートカットを設定した画面から右クリックでもその場所が開けます）

で、さっそく
shift+command+7
を押してみたときに、保存先のフォルダが事前に作成されてないとエラーになります。
あと、場合によっては、スクリーンショットを撮る権限を、Automatorに与えないといけないです。↓のように、システム環境設定の、Security & PrivacyのScreen Recordingの権限を、Automatorに与えます。
これで完了です。

なお新しいサービス（Quick Action）をAutomatorで作成したのに、ショートカット設定画面の一覧に出てこない場合は、システム環境設定のウィンドウをいったん落として開き直すと、出てきたりします。
あと、このやりかただと、まずないとは思いますが1秒以内に2回撮影すると、1回目のは上書きで消えてしまいます。

*1:私はMacの言語設定を英語にしてるのですが、英語の勉強のためでも、カッコつけてるわけでも、ましてや英語のほうが理解しやすいからでもなく、分からないことを調べるときに英語の機能名称とかでググったほうが圧倒的に情報が多いからです。一時期、なぜか英語設定にしてないと無料版がインストールできないソフトとかがあったってのもありますが。

Pythonでよく忘れる、よく間違える書き方

• Pandasのデータフレームは、=でコピーしようとすると、コピーじゃなくて参照渡しになるので、コピーしたつもりのdfを処理すると元のdfも処理されてしまう。df2 = df1.copy()とするのを忘れないように。
• Pandasで要素がNaNかどうかを判定させようとする時、要素を指定してその後ろに.isnull()メソッドを使うんじゃなくて、データフレームに.isnull()を実行してそのあとで要素を指定する。

df.isnull().iloc[1,1]

• リストのappendは、直接編集するメソッドなので、=で代入しなおさなくてよい。（リンク）
• リストをソートする時も、.sort()をつけるだけで、代入はしなおさなくてよい。
• Pandasから.to_csvするとき、「quoting=csv.QUOTE_NONNUMERIC」という引数をつけると文字列のフィールドは""をつけることができるが、これは事前にimport csvしてないといけない。
• リストの要素を置換する時、置換対照表を辞書にしておいて内包表記を使うと簡単にかけたりする。

X = ['あ','い','あ','あ','う','え','お','あ']
Y = {'あ':'ア','お':'オ'}
X2 = [Y[x] if x in Y.keys() else x for x in X]
print(X2)

• listをユニークにするには、NumPyなら.unique()が使えるが、setにしてからlistに戻せばユニークにはなる。順序は崩れるので注意。

unique_list = list(set(original_list))

• 数値をゼロ埋め（ゼロパディング）した文字列がほしい時は、数字を文字列にしてから、たとえば.zfill(3)とつけるとゼロ埋め3桁の数字（文字列だけど）になる。
• Rばかりやってると間違えるのだが、Pythonでは代入はオブジェクトのコピーを意味せず、名前が増えるだけとなる。つまり、代入した先の変数名に対して編集すると、元のオブジェクトが編集されてしまう。

x1 = [1,2,3,4,5]
x2 = x1
x2.remove(3)
print(x1)
print(x2)
[1, 2, 4, 5]
[1, 2, 4, 5]

• andとorを組み合わせると、andが優先（先にくくられる）なので、A and B or C and Dだと(A and B) or (C and D)の意味になる。明示的に優先させたければ()でくくる。
• nanとnanは==マッチしないので、nanかどうか判定するにはmath.isnan()を使う必要がある。
• object型の列は、中に複数種類のデータ型が含まれている可能性がある。特に、pandasで文字列型の列を読み込むとobject型になるが、nanを含んでいる場合、そのnanが別の型になってたりとか。
• torch.tensorの次元や長さについて。

# 次元数の確認（何重の入れ子になってるか）
x.dim()
x.ndimension()
x.ndim

# 各次元ごとのサイズを確認（外側の次元から順番に要素数が表示される）
x.size()  # numpy.arrayだとこれでバラしたときの要素数が返る？ので注意
x.shape

# 要素数（次元をバラしたときの要素数）
x.numel()
x.nelement()

• 文字列の出力を、下に並べていくんじゃなくて上書きさせたいときは、キャリッジ・リターンで先頭に戻す動きをすればよい。（前の出力のほうが長くてケツのほうが残ったりするのを避けるには、スペースで埋める。）

for i in range(100):
    print(f'{i}週目です', end='\r')
    time.sleep(5)

• データフレームを複数の条件でソートするとき、Rと違って、高階層側（外側）の条件を先に書く。以下のようにすると、全体としてはmatch_id順に並んでいて、同じmatch_idの中でplayer名順に並ぶようになる。Rとかエクセルだと内側から先にソートするので注意。

d2 = d2.sort_values(by=['match_id','player'], ascending=True)

• pandasのデータフレームをソートしたあと、reset_index(drop=True)しておかないと、順番は代わってるのにインデックスのラベルはもとのままになっているので、あとでおかしなことになる。

d2 = d2.sort_values(by=['match_id','player'], ascending=True)
d2 = d2.reset_index(drop=True)

Rでジニ係数を算出

めちゃくちゃ簡単なしょうもない内容ですが、あとで個人的に使うので、ジニ係数を出す簡単な関数をメモしておきます。
ここでは例として、都道府県の人口データのジニ係数を出してみます。
データを小さい順に並べて、都道府県数の累積（これは単順番を表すindexと同じになる）の累積構成比をx軸、人口の累積構成比をy軸にとったのがローレンツ曲線で、そっから三角形とスイカ形の面積比を求めるわけですね。

> # 北海道から沖縄県までの人口（千人単位）が順番に入ったデータ
> print(pop)
 [1]  5286  1263  1241  2316   981  1090  1864  2877  1946
[10]  1952  7330  6255 13822  9177  2246  1050  1143   774
[19]   817  2063  1997  3659  7537  1791  1412  2591  8813
[28]  5484  1339   935   560   680  1898  2817  1370   736
[37]   962  1352   706  5107   819  1341  1757  1144  1081
[46]  1614  1448
> 
> # ジニ係数を取得する関数
> get_gini <- function(v){
+   v.sorted <- v[order(v, decreasing=F)]
+   x <- 1:length(v)
+   y <- cumsum(v.sorted)/sum(v.sorted)
+   AB <- ((length(v))/2)  # 三角形の高さは1なので2乗しない
+   B <- sum(y)
+   gini <- (AB-B)/AB
+   return(gini)
+ }
> 
> # グラフを描く関数
> plot_gini <- function(v, col="black", lty=1, lwd=1){
+   v.sorted <- v[order(v, decreasing=F)]
+   x <- c(0, (1:length(v))/length(v))
+   y <- c(0, cumsum(v.sorted)/sum(v.sorted))
+   plot(x=x, y=y, type="l", col=col, lty=lty, lwd=lwd)
+   par(new=T)
+   plot(x=x, y=x, type="l", col="black", axes=F, ann = F)
+ }
> 
> # ジニ係数を出す
> get_gini(pop)
[1] 0.4415188
> 
> # 図示する
> plot_gini(pop)

斜め45度の線と赤い弧に囲まれた部分の面積が、三角形に占める割合がジニ係数で、1に近づくほど偏りが大きいことになりますね。

ちなみに、1920年、1980年、2018年を比較すると人口ジニ係数は0.248, 0.396, 0.442とどんどん上がっており、図示すると以下のような感じでした。

> # pop1920, pop1980, pop2018はそれぞれの年の都道府県人口
> plot_gini(pop1920, col="#00CC00", lty="dashed")
> par(new=T)
> plot_gini(pop1980, col="#0000CC", lty="dashed")
> par(new=T)
> plot_gini(pop2018, col="#CC0000", lwd=1.5)
> legend(
+   "bottomright",
+   legend=c("1920年", "1980年", "2018年"),
+   lty=c("dashed", "dashed", "solid"),
+   lwd=c(1,1,1.5),
+   col=c("#00CC00", "#0000CC", "#CC0000")
+   )

肺炎死者数の月別データで季節調整を試してみる（RからX-13ARIMA-SEATSを利用）

Rでの季節調整をやってみます。
季節調整は、arima関数のseasonal引数を指定したSARIMAモデルや、stl関数でできてそっちは簡単なのですが（参考1・参考2・参考3）、一般的にアメリカ商務省センサス局の「X-12-ARIMA」とか「X-13ARIMA-SEATS」とかが有名なので、それを使ってみようと思います。

1965年に開発された「X-11」は移動平均をつかった成分分解法で、これにRegARIMAモデルをつかって事前にデータ補完する機能などを加えたX-12-ARIMAが1996年に開発され、さらに2012年に、X-11とは別系統の技術として利用されてきたARIMAモデルベースの分解方法であるSEATSを選択肢として使えるようにしたX-13ARIMA-SEATSがリリースされた……ということらしい。なのでX-13はX-12が新しくなったというより、X-12とは別の技術も使えるようにしたパッケージという感じのようです。

X-12-ARIMAは、大きく分けて、異常値等を除去するための回帰変数の決定とARIMAモデルによる予測値を推計する「事前調整パート（RegARIMA）」、X-11による季節調整を行う「移動平均パート」、季節調整系列の適切性や安定性に関する診断を行う「事後診断パート」の３つのパートから構成されている。
https://www.stat.go.jp/training/2kenkyu/ihou/71/pdf/2-2-712.pdf

ということなんですが、要するに季節調整のコア部分にはX-11が使われ続けていて、X-12というのはそれに便利機能を付加したものというイメージです。

Rのパッケージだけではできないので、インストール方法からメモしておきます。具体的には、X-13ARIMA-SEATSをインストールして、それをRの{seasonal}パッケージからよびだして使う形になります。{x12}というパッケージからも使えるので追加的に後でやります。
（OSはMacで、X-code、commandline tools、homebrewがインストールされていることを前提とします）

まず、X13-ARIMAのコンパイルにはgfortranというコンパイラが必要らしいですが、これはgccに統合されているとのことなので、gccをインストールすればよい。homebrewで入ります。

$ brew install gcc

次に、センサス局のサイトにいって、X13のソースコードをダウンロードしてきます。
X-13ARIMA-SEATS Download Page (Unix/Linux Version) - X-13ARIMA-SEATS Seasonal Adjustment Program - US Census Bureau
このURLにある、「x13ashtmlsrc_V1.1_B39.tar.gz」という圧縮ファイルをダウンロードして、解凍するとフォルダが出てきますので、それを適当な場所に置きます。私は、
~/Library/R
に置きました。

で、この記事に書いてあるのですが、
Compiling X 13ARIMA SEATS from Source for OS X · christophsax/seasonal Wiki · GitHub
コンパイル前に、makefileの中身を少し書き換えないと、エラーが出ます。さっき解凍したフォルダの中に、
makefile.gf
というファイルがあるので、これをテキストエディタで開いて、293行目の"-static"を消します。

削除前：$(LINKER) -static -o $@ $(OBJS) $(LDMAP) $(LIBS) $(LDFLAGS)
削除後：$(LINKER) -o $@ $(OBJS) $(LDMAP) $(LIBS) $(LDFLAGS)

あとはターミナルからmakeすればいいです。

$ cd ~/Library/R/x13ashtmlsrc_V1.1_B39
$ make -f makefile.gf

コンパイルが終わると、フォルダの中に
x13asHTMLv11b39
という実行ファイルが出てきます。こいつがX13-ARIMAの処理をやってくれるのですが、Rの{seasonal}パッケージは、ディレクトリを指定してやるとその中にある「x13ashtml」というファイルを探すようになっているので、ファイル名を「x13ashtml」に変えておきます。

$ mv x13asHTMLv11b39 x13ashtml

ここまできたら、あとはRのほうから設定を行います。

> library(seasonal)  # パッケージを読み込み
> Sys.setenv(X13_PATH = "~/Library/R/x13ashtmlsrc_V1.1_B39")  # 実行ファイルの置き場を指定
> checkX13()
X-13 installation test:
  - X13_PATH correctly specified
  - binary executable file found
  - command line test run successful
  - command line test produced HTML output
  - seasonal test run successful
Congratulations! 'seasonal' should work fine!

上のように、checkX13()を実行して上記のようなメッセージが出たら、セットが成功していて、X13をRから呼び出せるようになっているということになります。インストールに不具合があった場合、ここでエラーが出ます。

季節調整ができるかどうか、肺炎の月別死者数のデータをつかって試してみます。
データは、以下のURLに載っている「人口動態統計月報（概数）」を取りました。平成17年1月から、死因別の死者数が載っていて、前年分も書いてあるので、平成16年1月から数字が手に入りました。
人口動態調査 結果の概要｜厚生労働省

なお、2017年1月から死因の分類変更があって、それまで肺炎に分類されていた人が一部、神経系や認知症系の方に分類されるようになったらしくて、内訳をみて揃えることができなかったので、2017年1月以降は1とするダミー変数も作っておきました。あとで外生変数として使います。

データを読み込んで、とりあえず原系列のグラフ描いてみます。

# CSVの読み込みとts型に変換
# CSVの1列目は年、2列目が肺炎の死者数、4列目に分類変更後ダミーを入力しといた
# （3列目は、ここで使いませんが、自殺者数を入れてある）
d <- read.csv('source/death.csv')
d <- ts(d, start=c(2004,1), end=c(2019,10), frequency = 12)

# 原系列のグラフ
plot(d[,2], xlab="", ylab="肺炎死者数", xaxp=c(2004, 2020, 8))

y軸の目盛りの省略が気持ち悪ければylimを設定すればいいと思います。
これに、季節調整をかけてみます。ひとまずダミー変数を使わずに、死者数の系列をそのまま使います。
X-12も含めて季節調整の原理については『入門 季節調整―基礎知識の理解から「X‐12‐ARIMA」の活用法まで』とかに詳しく書いてあります（ちゃんとよんでないけど）。

まずは推定して結果のサマリをみてみます。

> pneu.x13.1 <- seas(x=d[,2])
> summary(pneu.x13.1)

Call:
seas(x = d[, 2])

Coefficients:
                   Estimate Std. Error z value             Pr(>|z|)    
Constant           0.022586   0.003928   5.751      0.0000000088872 ***
AO2005.Mar         0.217141   0.032563   6.668      0.0000000000259 ***
LS2015.Feb        -0.114776   0.028185  -4.072      0.0000465558118 ***
LS2017.Jan        -0.246416   0.027709  -8.893 < 0.0000000000000002 ***
AR-Nonseasonal-01  0.696167   0.052020  13.383 < 0.0000000000000002 ***
MA-Seasonal-12     0.763174   0.053533  14.256 < 0.0000000000000002 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

SEATS adj.  ARIMA: (1 0 0)(0 1 1)  Obs.: 190  Transform: log
AICc:  2644, BIC:  2665  QS (no seasonality in final):    0  
Box-Ljung (no autocorr.): 30.71   Shapiro (normality): 0.9937  
Messages generated by X-13:
Warnings:
- At least one visually significant trading day peak has been found in one or more of the estimated
  spectra.

「coefficients」の中にある、「AO2005.Mar」というのは、「Additive Outlier」というタイプの異常値が2005年3月に検出されていて、そのダミー変数の係数が推定されてるわけですね。これはいわゆるスパイク的な変化のことでしょう。「LS2015.Feb」と「LS2017.Jan」は「Level Shift」というタイプの異常値が検出されていて、名前はそのまんまで、水準が階段のように上がったり下がったりする変化です。2017年1月のほうは、分類変更の影響です。なおこれらの他に、増減傾向の変化を表す「TC: Trend Change」というタイプの異常値もありますが、これは今回は検出されませんでした。
「ARIMA: (1 0 0)(0 1 1)」とあるのは、いわゆるARIMA(p,d,q)(P,D,Q)の次数ですね。pはARの次数、dは何階和分過程か（何階階差を取って分析するか）、qはMAの次数を表しています。小文字のほうは時系列方向のARIMAの次数で、大文字のほうは季節方向の次数です。季節方向というのは言い方が難しいですが、年を行、月を列とした表とつくった場合に縦方向に眺めるようなイメージで、たとえば2次のAR(2)だったらその部分はYtがY_t-12、Y_t-24で説明されるということです。
「Transform: log」は、デフォルトで自動変換が選択されてて、ここでは対数変換されたようです。対数変換されているばあい、いわゆる「乗法型季節調整」になり、原系列=トレンド成分*季節成分*不規則成分という分解がされていることになります。

グラフを描いてみます。

# グラフをかく
par(mfrow = c(4,1), oma=c(3,2,2,0), mar=c(0,2,0,2), mgp=c(1,1,0))
plot(pneu.x13.1$series$s10, xlab="", ylab="季節成分", xaxt="n", yaxt="n")  # 季節成分
plot(pneu.x13.1$series$s11, xlab="", ylab="季節調整値", xaxt="n", yaxt="n")  # 季節調整値
plot(pneu.x13.1$series$s12, xlab="", ylab="トレンド成分", xaxt="n", yaxt="n")  # トレンド成分
plot(pneu.x13.1$series$s13, xlab="", ylab="不規則成分", yaxt="n", xaxp=c(2004, 2020, 8))  # 不規則成分
par(default.par)

こんな感じ分解されましたが*1、季節調整値（2段目）やトレンド（3段目）をみると、2017年に急激に下がってますね。さきほどの「分類変更」の影響です。
なお、s10とs12とs13をかけると、原系列に戻せます。また、季節調整値というのは、原系列を季節成分で割ったものになります。（そのへんの話はこの資料をみると分かる。）
ここでは貼り付けませんが、季節成分の抽出結果のテーブル（行が年、列が月）を眺めていると、肺炎死者数はやはり2月や1月が多く温かい時期に減りますが、不思議なことに8月は一貫して7月や9月に比べて増える傾向がある。あと、冬場の肺炎死者数と、4月や9・10月といった春秋の肺炎死者数の格差が、年々小さくなっていました（冬場のウェイトが小さくなって春秋のウェイトが大きくなっている）。

ちなみに、推定結果のオブジェクトをそのままplotに投げると、以下のように、原系列と季節調整値を重ねて、かつ3種類のoutlier（異常値）を図示したものが出てきます。

plot(pneu.x13.1)

予測値と信頼区間もついたグラフを描くならこうします。（描画に凝るならggplot2のほうがいいかも）

### 予測値ありのグラフを描く
# 予測値の取り出し
# 1列目が予測、2列目がlower、3列目がupper
pneu.x13.f1 <- series(pneu.x13.1, c('forecast.forecasts')) 

# 元の時系列のプロット
plot(pneu.x13.1, xlim=c(min(time(d)), max(time(pneu.x13.f1))))

# 信頼区間の塗りつぶし
# 信頼区間下限を左からなぞり、上限を右からなぞる設定でxy座標を作成しpolygonを描画
# time(d)[nrow(d)]とd[nrow(d),2]を入れるのは、実績の最後と
# 予測の最初のあいだに、隙間を開けさせないため
x.ci <- c(time(d)[nrow(d)], time(pneu.x13.f1), rev(time(pneu.x13.f1)), time(d)[nrow(d)])
y.ci <- c(d[nrow(d),2], pneu.x13.f1[,2], rev(pneu.x13.f1[,3]), d[nrow(d),2])
polygon(x.ci, y.ci, col = '#ccccff', border = NA)

# 予測値のプロット
# 隙間を開けさせないために、実績値の最後の値と、予測値の系列をつないで新たなtsオブジェクトにする
pred <- ts(c(d[nrow(d),2], pneu.x13.f1[,1]), start=c(2019,10), end=c(2022,10), frequency=12)
lines(pred, col="#5555ff", lwd=1, lty=1)

次に、2017年1月以降とそれ以前を区別するダミー変数をつかった推定をしてみたいと思ったのですが、なんかエラーが出まくりました。
下記の記事に書いてあるように、外生変数を使う場合、将来予測用に余計なデータが求められて、ややこしいことになるらしい（X13の仕様により）。そこで、x11という引数もあわせて指定すると動くのですが、要するにこれは、成分分解をSEATSではなくX-11で行うということでしょうね。
NOTE not shown · Issue #200 · christophsax/seasonal · GitHub

とりあえずやってみます。
xregという引数に、外生変数を指定します。

> d_sadj.2 <- seas(x=d[,2], xreg=d[,4], 
+                forecast.maxlead = "0", series.span=", 2019.10", x11 = "")
> summary(d_sadj.2)

Call:
seas(x = d[, 2], xreg = d[, 4], forecast.maxlead = "0", series.span = ", 2019.10", 
    x11 = "")

Coefficients:
                   Estimate Std. Error z value             Pr(>|z|)    
xreg              -0.246421   0.027708  -8.894 < 0.0000000000000002 ***
Constant           0.022586   0.003928   5.750      0.0000000089088 ***
AO2005.Mar         0.217145   0.032563   6.668      0.0000000000259 ***
LS2015.Feb        -0.114767   0.028183  -4.072      0.0000465725741 ***
AR-Nonseasonal-01  0.696129   0.052021  13.382 < 0.0000000000000002 ***
MA-Seasonal-12     0.763082   0.053542  14.252 < 0.0000000000000002 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

X11 adj.  ARIMA: (1 0 0)(0 1 1)  Obs.: 190  Transform: log
AICc:  2644, BIC:  2665  QS (no seasonality in final):    0  
Box-Ljung (no autocorr.): 30.71   Shapiro (normality): 0.9937  
Messages generated by X-13:
Notes:
- Unable to test LS2017.Jan due to regression matrix singularity.
- Unable to test LS2017.Jan due to regression matrix singularity.

xregが説明変数として導入したダミー変数で、係数は先ほどのLS2017.Janとほぼ同じ値になってますね。そしてLS2017.Jan は消えています。（一番下のNotesに、重複するので削る旨が書いてあります。）
次にグラフです。

par(mfrow = c(4,1), oma=c(3,2,2,0), mar=c(0,2,0,2), mgp=c(1,1,0))
plot(d_sadj.2$series$d10, xlab="", ylab="季節成分", xaxt="n", yaxt="n")  # 季節成分
plot(d_sadj.2$series$d11, xlab="", ylab="季節調整値", xaxt="n", yaxt="n")  # 季節調整値
plot(d_sadj.2$series$d12, xlab="", ylab="トレンド成分", xaxt="n", yaxt="n")  # トレンド成分
plot(d_sadj.2$series$d13, xlab="", ylab="不規則成分", yaxt="n", xaxp=c(2004, 2020, 8))  # 不規則成分

seriesから取り出す値が「d10」とかになってるんですが、これはX-11で分解（decompositionのdかな？）された結果のデータです。SEATSでの分解結果は、さっきやったように、「s10」とかで取り出される。
トレンド成分をみると、2017年1月の分類変更は取り除かれている感じになっています。「季節調整値」だけだと、単なる季節調整の結果なので、分類変更の影響が出てしまいますが。
なお、上では示してませんが、ダミー変数にかかる係数の推定結果とかも格納されてて、将来予測を出すときはそれも使われることになります。

性能とかは判断できませんが、インストールも面倒だし、X13の仕様を詳しく見ないと使いこなせないので、ほかのもっと基本的なパッケージ・関数を使った方が楽な気はしますねｗ

ちなみに、↑でインストールしたX13のエンジンは、{x12}というパッケージからも使えます。
外生変数を使わないなら、{x12}のほうがplotの描き方が楽でいいと思いました。

library(x12)
x12path("~/Library/R/x13ashtmlsrc_V1.1_B39/x13ashtml")  # 実行ファイルまで記述

d_x12.1 <- x12(d[,2])

plot(d_x12.1, ylim=c(0, 14200),
     original=T, col_original="black",
     sa=T, col_sa="blue",
     trend=T, col_trend="red", lwd_trend=1.5)

ただ、さっきの分類変更を反映するダミー変数みたいに、回帰の変数を使った推定をしようとするのであれば、設定がいろいろややこしそうです。そもそもX12やX13の機能とシンタックスに関する基礎知識が必要になる模様。
X12やX13の基本が分かっていると、かなりいろいろな要因を季節調整時に考慮できるみたいですが、少し調べたぐらいでは使えそうになかったので、ひとまず放置しました。季節調整を凝らなければならない用事ができたときにまた勉強しようと思います……。

ちなみに自殺者数はこのようになってました。

全体としては減り続けていて、東日本大震災のあった2011年3月に大きく減少し、その後5月に激増しました。

参考

Rの{seasonal}パッケージの仕様
https://cran.r-project.org/web/packages/seasonal/seasonal.pdf

米センサス局のX-13ARIMA-SEATSの仕様
https://www.census.gov/ts/x13as/docX13AS.pdf

X12に関する日本語の説明サイト
http://www.cirje.e.u-tokyo.ac.jp/research/dp/2001/2001cj47.pdf
https://ecitizen.jp/x-12-arima/x-12-arima/

Rの{x12}パッケージの使い方の解説資料（テキストとスライド）
https://www.jstatsoft.org/article/view/v062i02
https://www.r-project.org/conferences/useR-2015/presentations/135.pdf

Rの{x12}パッケージの仕様
https://cran.r-project.org/web/packages/x12/x12.pdf

X-11、X-12-ARIMA、X-13ARIMA-SEATSの関係の解説
http://www.esri.go.jp/jp/archive/snaq/snaq150/snaq150d.pdf

X-11のベースになっている移動平均法による加法型・乗法型の季節調整の原理については、下記資料のp.79あたりからをみると分かりやすい。X-11の問題点*2をX-12がどのように解決するかの説明も書かれてある。
https://www3.boj.or.jp/josa/past_release/chosa199605e.pdf

*1:目盛りの数字が重ならないように設定するのが面倒だったので、y軸の目盛りを消してしまってます。4つのグラフで縦方向のスケールが全然違うのでそこは注意が必要。

*2:前後7年移動平均を取るが末端から7年以内のところではデータがないので後方移動平均に変わっていく点や、妥当性の事後評価ができない点。

四半期GDP成長率の「年率換算」はどれぐらいブレるのか？

GDPは「四半期速報値」というのが作成されており、前の四半期との比（たとえば2019年7-9期と10-12月期の比）を4乗したものが、「年率換算」の成長率として報道などで使われている。最新の速報では、2019年10-12月期の実質GDPが、年率換算でマイナス7.1%と大幅な減少を記録したと報道されていた。

この、四半期成長率を4乗した値というのは、トレンドが安定していればもちろん年率の成長に一致するが、短期的な変動が4乗されるわけなので、けっこう大きく上下する。

以下は、ちょっとわかりにくいが、実質GDPの
-年間（暦年）の成長率［黒い線］
-四半期成長率を「年率換算」したもの（前期比の4乗）［赤い線］
-四半期GDPを4倍したものと、前年の年間GDPの比［青い線］
をグラフにしたものだ。四半期はいずれも季節調整済み。

データは、暦年GDPについては平成30年度の国民経済計算年次推計、四半期速報値は2019年10-12月期の一次速報値（実質・季節調整済み）（年率換算マイナス6.3%と報道されたやつ）を用いた。

年成長率（黒線）と、四半期成長率の年率換算値（赤線）は、けっこう乖離があるということがわかる。ちなみに、「年成長率」と「四半期成長率の年率換算」の相関係数をとると0.33ぐらいになる。

最新の速報では2019年10-12月期の年率換算成長率がマイナス7.1%になったということだが*1、リーマンショックの影響が出た2009年1-3月期は年率換算で-17.7%、今回と同じく消費増税があった2014年4-6月期は年率換算-7.4%だったので、今回はまあ、それらよりはマシということになる（それでも大きな落ち込みだが）。次の2020年1-3月はコロナの影響があるので、リーマン並かそれ以上いくのかも？

一方、「四半期GDPを4倍したもの」と「前の年のGDP」の比をとったのが上のグラフでいうと青い線なのだが（「前年同期比」ではない点に注意）、これは年成長率とあまり大きくは乖離していない。年成長率との相関係数も、0.93ぐらいになる。

もちろん、短期的の成長を瞬間風速的にみたい場合はもちろんあるだろうから、「四半期成長率の年率換算」に意味がないわけではない。というか、「年成長率」ですら、人によっては「短期的すぎる」指標かもしれないので、まあ要するにみる人の関心次第としか言いようがない。また、上のグラフでいう青い線のほうが優れた性質のデータだ、というわけでもない。赤と青は意味合いが異なるわけで、単に「年成長率との近さ」を知りたいなら、青い線のほうをみたほうが近くなるという話である。

重要なのは、四半期成長率の年率換算というのは、あくまで「その3ヶ月と前の3ヶ月のあいだの短期的な変動」を誇張したものなので、「年率換算」したからといって、「年スケールでのトレンド」を見通そうとしているわけでは全くないということ。まぁ当たり前だが。

上のグラフで確認したかったのは、「四半期スパンの変動と年スパンの変動はどれぐらい乖離しているのか」であって、四半期GDPをみることが役に立つか立たないかということではない。

ただ、年成長率とはあまり関係ないものを見てるのであり、4乗したところで年成長率に近くなってるわけでもないのだから、換算しなくていい気はするが。四半期の速報値から年スパンでの変化を推し量りたいのであれば、上のように前の年の年間GDPと比較するか、過去4四半期の平均をその前の4四半期の平均と比較するかすればいいんじゃないだろうか。

会社員時代、営業データを週次・月次・年次で報告する仕事をしてたのだが、その経験からすると、上のグラフの赤い線（四半期GDPの年率換算）のように大きくジグザグする指標は、直接モニタリングするのにはあまり向いていない可能性もある。時系列データをどういう単位でみるべきかというのは、考え方はいろいろあり、もちろん短期的なショックを検出することも大事ではある。しかし、大きく下がったあとに反動で大きく上がるような数字は、天候とか平日日数*2みたいな特殊要因に大きく左右されてるわけで、解釈がめんどくさい。というか「短期的な解釈」だけをして「トレンドとして解釈」してしまわないような注意がいる。トレンドはトレンドで別途確認するという心がけが必要だろう。

*1:上のグラフは、最新の速報の前の、「マイナス6.3%」になってる1次速報データで作っている

*2:これは四半期GDP統計では考慮済みらしいが

Macのパーティションがおかしくなった場合の対処

Mac OS（最近はmacOSというらしいが）のディスクユーティリティでパーティションを変更しようとするとき、よくわからない現象に直面することがたまによくあります。
下記の記事のように、一部のボリュームが消してくても消せない（マイナスボタンがない）状態になったりします。
mac GUIで結合できないパーティションを結合する – MY ROBOTICS

今回、外付けHDDのパーティションをいじろうとしたらそれに陥ったので、ターミナルからdiskutilコマンドで操作してみたのですが、それでも変なことがおきました。
具体的には、

$ diskutil list

/dev/disk2 (external, physical):
   #:                       TYPE NAME                    SIZE       IDENTIFIER
   0:      GUID_partition_scheme                        *4.0 TB     disk2
   1:                        EFI EFI                     209.7 MB   disk2s1
   2:                  Apple_HFS TM_MBP2016_Portable     1.0 TB     disk2s2
   3:                  Apple_HFS Untitled                699.9 GB   disk2s4
   4:                  Apple_HFS Strage_main             2.3 TB     disk2s5

というふうに、確認すると（diskutil listの出力結果はexternalのところだけ貼ってます）、disk2という外付けHDDの中にEFI（システム側で必要なやつ？）を入れて4つのパートがあって、足したら4TBになります。これのうち、disk2s5だけ残してdisk2s2とdisk2s4を削除し、結合して（EFIを含めて）4TB使えるようにしたかったわけです。要するにパーティションをやめて全部結合するということ。

そこで、

$ diskutil eraseVolume "Free Space" not disk2s2
$ diskutil eraseVolume "Free Space" not disk2s4

というふうにボリュームの消去をしてみたのですが、その結果どうなったかというと、

/dev/disk2 (external, physical):
   #:                       TYPE NAME                    SIZE       IDENTIFIER
   0:      GUID_partition_scheme                        *4.0 TB     disk2
   1:                        EFI EFI                     209.7 MB   disk2s1
   2:                  Apple_HFS Strage_main             2.3 TB     disk2s2

こういうふうに、不要なボリュームは消えてくれました。ディスク全体としては4.0TBあって、残存するボリュームの合計が2.3TB+209.7MBであり、残りがFree Spaceになったのだと思います。
で、現在2.3TBになっているdisk2s2をリサイズしてマックスまで使えるようにしたかったのですが、ディスクユーティリティで見てみると、

こんなふうに、Free Spaceは表示されなくて、全体で2.3TBしかないよということになってしまいました。
そこで、ネットでみた説明を参考にdiskutil resizeVolumeというコマンドをつかって、リサイズ後のサイズを「100%」にしてやってみたのですが、

$ diskutil resizeVolume disk2s2 100%
Resizing to full size (fit to fill)
Started partitioning on disk2s2 Strage_main
Error: -69743: The new size must be different than the existing size

こういう結果になりました。100%だともとのサイズと同一という意味であって、ディスクをマックスまで使うという意味じゃなかったんですねｗ
なので、具体的に数値で指定してやってみたのですが、たとえば3.9TBに設定したとしても、

$ diskutil resizeVolume /dev/disk2s2 3.9TB
Resizing to 3900000000000 bytes
Started partitioning on disk2s2 Strage_main
Error: -69519: The target disk is too small for this operation, or a gap is required in your partition map which is missing or too small, which is often caused by an attempt to grow a partition beyond the beginning of another partition or beyond the end of partition map usable space

となり、容量が足りませんと出ます。3.9TBにEFIの分を足しても超過するわけではないはずなのですが・・・。
結局、もう調べるのがめんどくさくなったので、
Macのディスクユーティリティーコマンドで外部ストレージをフォーマット - ytooyamaのブログ
を参考に、論理ボリュームではなく物理ディスク全体を消去することにしました。データが全部きえますが、別のディスクにコピーしてあるので後で戻します。

$ diskutil eraseDisk JHFS+ Untitled /dev/disk2

そしたらまぁ、

こんなふうに、4TB使える状態にはなりました（画像はボリューム名をUntitledからStrage_mainに戻したあとのものです）。
Free Spaceと結合して論理ボリュームのサイズをでかくする方法は、後でわかったら追記します。

対数尤度が正で、AICが負になるケース

初歩的な話ですが、今日学生から訊かれたので例を考えてみました。
AICがマイナスの値になってもいいんだっけ？みたいな話です。

AICがマイナスに

統計モデルの最尤推定をする場合、尤度は1より小さい場合が多く、したがって対数尤度が負である場合が多いとしましょう*1。特に離散変数では必ずそうなります。その場合、AICは「-2×(対数尤度-パラメータ数)」なので、正になるでしょう。
でも、対数尤度が正でAICが負になることというのも、連続変数ではあり得ます。その場合AICは、負で絶対値が大きいほどよい、ということになる。

誤差の分布は正規分布だと考えているとして、たとえば、

• モデルの当てはまりがめちゃめちゃいい
• そもそも数値の単位が小さい

というような場合に、誤差の標準偏差（分散）の値が小さくなって、結果的にAICが負になりやすくなりますね。
　
　

正規分布の確率密度

統計モデルの誤差が平均ゼロ、標準偏差σの正規分布に従うとした場合、確率密度は誤差ゼロのときに最大になるわけですが、誤差ゼロのときの確率密度は、σが小さいと大きくなります。
Rでは、dnorm()という関数で正規分布上の確率密度が簡単に計算できますが、平均ゼロの正規分布上でゼロという値の確率密度を、標準偏差が1、0.1、0.01の3パターンで計算してみると、

> # sdは標準偏差
> dnorm(0, mean=0, sd=1)
[1] 0.3989423
> dnorm(0, mean=0, sd=0.1)
[1] 3.989423
> dnorm(0, mean=0, sd=0.01)
[1] 39.89423

というふうに、標準偏差に反比例して大きくなっていきます。正規分布の確率密度関数でいうと、も（平均）も0で、右のところがになるので、まぁ当たり前です。

平均値に一致する点（誤差ゼロの点）でなくても、たとえばプラス方向に1標準偏差分行ったところの値も、同じようになります。この場合はexpのところがexp(-0.5)に固定になりますからね。

> dnorm(1, mean=0, sd=1)
[1] 0.2419707
> dnorm(0.1, mean=0, sd=0.1)
[1] 2.419707
> dnorm(0.01, mean=0, sd=0.01)
[1] 24.19707

標準偏差（Sigma）が0.5の場合と0.1の場合について確率密度関数のグラフを書いてみると次のようになります。

x <- seq(from=-1, to=1, by=0.01)
Density <- dnorm(x, mean=0, sd=0.5)
plot(x, Density, type='l', col='blue', ylim=c(0,4))

Density <- dnorm(x, mean=0, sd=0.1)
par(new=T)
plot(x, Density, type='l', col='red', ylim=c(0,4))

text(x=0.1, y=3, labels='Sigma=0.1', col='red', pos=4)
text(x=0.2, y=1, labels='Sigma=0.5', col='blue', pos=4)

標準偏差0.5の場合は、xがどんな値であっても確率密度が1より小さいので、それを掛け算してモデルの尤度が計算されるわけですから、1より小さい値でどんどん小さくなっていきます。当てはまりが良くないモデルだと、最尤法で計算する場合の母標準偏差の推定値も比較的大きくなって、こういう横に広がった分布が想定されることになります。
一方、標準偏差0.1の場合は、平均値まわりの割と広い範囲で確率密度が1を超えてますね。
　
　

当てはまりが良いとAICがマイナスになりやすくなる

切片ゼロ、傾き1（α=0、β=1）の線形モデルに、平均ゼロ、標準偏差sigmaの正規分布に従う乱数を加えてデータをつくって、単回帰分析を行い、そのモデルの決定係数、対数尤度、AICを計算してみます。まずsigma=1としてみます。

> # データを生成
> sigma <- 1
> x <- seq(from=0.1, to=10, by=0.1)
> y <- x + rnorm(n=100, mean=0, sd=sigma)
> 
> # 回帰モデルを推定して、対数尤度とAICを計算する
> lm1 <- lm(y~x)
> r2 <- summary(lm1)$r.squared
> lm1.loglik <- logLik(lm1)
> lm1.aic <- AIC(lm1)
> print(r2)
[1] 0.896859
> print(lm1.loglik)  # logLikの結果はlog Likという特殊なクラスになってる
'log Lik.' -136.302 (df=3)
> print(lm1.aic)
[1] 278.604
> # 対数尤度とAICの関係
> -2*(as.numeric(lm1.loglik)-3)
[1] 278.604

対数尤度（log Lik）は負で、AICが正ですね。df=3となっているのは、単回帰モデルなのでパラメータが「切片」「傾き」「誤差の標準偏差」の3つなのを表してるんでしょう。
以下のようにグラフにしてみます。

title <- paste('Sigma=',as.character(sigma),
               ', R2=',as.character(round(r2,2)),
               ', LogLik=', as.character(round(lm1.loglik,2)),
               ', AIC=', as.character(round(lm1.aic,2)),
               sep='')

plot(x, y, pch=20, main=title)  # 散布図
abline(lm1, col='red', lwd=2)  # 推定された回帰直線を引く

同じやり方で、sigma=0.5、0.2、-0.1の場合も出してみると、次のようになります。

誤差分散（標準偏差）が小さくて当てはまりがめちゃめちゃいい場合は、AICが負になることが分かりますね。
ちなみに、標準偏差を0.01から2まで少しずつ変えていって、その都度乱数でデータを生成して回帰モデルを当てはめ、AICを計算するという作業を繰り返してその結果をグラフにすると、次のように対数の形なります。

get_aic <- function(sigma){
  x <- seq(from=0.1, to=10, by=0.1)
  y <- x + rnorm(n=100, mean=0, sd=sigma)
  model.1 <- lm(y ~ x)
  return(AIC(model.1))
}

aics <- c()
sigmas <- seq(from=0.01, to=2, by=0.01)
for (i in sigmas){
  aic <- get_aic(i)
  aics <- c(aics, aic)
}

plot(x=sigmas, y=aics, type='l',xlab='Sigma', ylab='AIC')

　
　

単純に、小さい数値を扱ってるとAICが負になりやすくなる

ところで、もう一つ、単にデータの値が小さい数字になるような単位だと、AICが負になりやすくなります。
これも同様に、誤差の分散が値としては小さくなるからですが、以下のように試してみます。
まず、誤差の標準偏差2のデータに単回帰モデルを当てはめて決定係数、対数尤度、AICを出してみて、元データを10で割ってから同じことをしてみます。

# 標準偏差2でやってみる
sigma <- 2
x1 <- seq(from=0.1, to=10, by=0.1)
y1 <- x1 + rnorm(n=100, mean=0, sd=sigma)
lm1 <- lm(y1~x1)
lm1.r2 <- summary(lm1)$r.squared
lm1.loglik <- logLik(lm1)
lm1.aic <- AIC(lm1)
title <- paste('Sigma=',as.character(sigma),
               ', R2=',as.character(round(lm1.r2,2)),
               ', LogLik=', as.character(round(lm1.loglik,2)),
               ', AIC=', as.character(round(lm1.aic,2)),
               sep='')
plot(x1, y1, pch=20, main=title)  # 散布図
abline(lm1, col='red', lwd=2)  # 推定された回帰直線を引く

次に、単にデータを10で割り算してから同じことをやってみます。
Sigmaは、自分で指定してないですが、lmの結果の中に入っている推定値を取り出して見ておきます。（それもだいたい10分の1になります。）

# データを10分の1に
x2 <- x1/10
y2 <- y1/10
lm2 <- lm(y2~x2)
lm2.r2 <- summary(lm2)$r.squared
lm2.sigma.est <- summary(lm2)$sigma  # σは推定値を取り出す
lm2.loglik <- logLik(lm2)
lm2.aic <- AIC(lm2)
title <- paste('Sigma=',as.character(round(lm2.sigma.est,3)),
               ', R2=',as.character(round(lm2.r2,2)),
               ', LogLik=', as.character(round(lm2.loglik,2)),
               ', AIC=', as.character(round(lm2.aic,2)),
               sep='')
plot(x2, y2, pch=20, main=title)
abline(lm2, col='red', lwd=2)

当てはまり具合（たとえばR2）は変わってないけど、AICは負になりました。

*1:だからこそ、マイナス符号をつけるAICが慣例化してるんでしょう。