戯言日記

Rの話だと思ったら唐突にサバゲーが混じってくる何か。

ジャンクのマテバを手に入れたので修理

f:id:doubtpad:20211024174046j:plain

という訳で、つい買ってしまった。

f:id:doubtpad:20211024174134j:plain
マテバ

以前から見た目がめっちゃ好みでずっと欲しかったところに、ハンマーの動作不良が原因のジャンク扱いで安く出てたのを見つけて、1日悩んでから結局ポチった。
修理できなかったらそのまま飾る予定。

f:id:doubtpad:20211024175149j:plain
全体

f:id:doubtpad:20211024181521j:plain
銃身

f:id:doubtpad:20211024181610j:plain
スイングアウト状態

リボルバーなのにオートマチックみたいな見た目とか、銃身が下とか、好きなポイントが多すぎる。
それとヘビーウェイト仕様なのでめちゃくちゃ重い。手元のハンドガンの中で一番かも。

ジャンクとのことだったので動作を確認。
トリガーは動くのだが、確かにハンマーが連動していない。というか完全に起きた状態で戻ってこなくなっている。
あとグリップに若干だけどガタついてる。この辺も直るようなら何とかしたいかも。


ざっくり状態が把握できたので、とりあえず中の状態確認のためにばらしていく。
説明書に分解画像とかパーツリストとかが無かったのでググったのだが、意外と情報が少なくて苦労したので少し細かく載せておく。

まずは銃口のリングを回して外す。正ねじ方向なので注意。
これを外すと銃身側の外装が引っこ抜ける。

f:id:doubtpad:20211024181902j:plain
銃身側の外装

シリンダーは本体フレームの正面側からマイナスネジを外すと引っこ抜ける。銃身を外さなくても作業できるが、外していた方がドライバーは回しやすい。
あとマイナスネジが意外と長い上にちょっと干渉して引き抜くのに苦労するので、少しずつ様子見しながら進めるといいかも。

f:id:doubtpad:20211024182101j:plain
シリンダー

シリンダーの各パーツは特にビス止めなどされていないので簡単にバラせる。

f:id:doubtpad:20211024182433j:plain
シリンダー各部

グリップは下からマイナスネジを1本外すだけで抜き取れる。

f:id:doubtpad:20211025234043j:plain
グリップ下部

グリップを引き抜くとガスタンクが剥き出しになる。
何となく冬場は冷えそうな見た目←

f:id:doubtpad:20211025234140j:plain
ガスタンク

バレルを押さえているパーツは銃口側から外す。
バレル周りの取り出しには影響しないので後回しでいいかも。
なお、ネジの上にあるビスはホップの調整ビスで、これを回して調整する。バラしていない時は銃身の上から六角レンチを差し込んで調整できる仕様。

f:id:doubtpad:20211025083753j:plain
バレルの安定用パーツ?

バレルはフレーム側にプラスネジ2本と六角ネジ1本で止められている。
これも銃口側からドライバーで外せる。

f:id:doubtpad:20211025083919j:plain
バレル部分

バレル後部にはスプリングとパーツがある。
本来はシリンダー内の弾頭とバレルに隙間が出来てしまうはずなのだが、この前後するパーツが微妙に前後することで隙間を埋めて上手く噛み合うようになっている。こういう「ものすごくしっかりと練られた構造」って見てるだけで楽しい。

f:id:doubtpad:20211025084014j:plain
バレル後部のスプリング

バレルそのものは根元の左側面から六角ネジを緩めることで外せる。

f:id:doubtpad:20211025084342j:plain
バレル周り

バレルに張り付いてるように見える黒いやつはホップ。
面ホップ仕様で、上から棒状のパーツで押すことでホップ調整されている。

f:id:doubtpad:20211025084508j:plain
面ホップ

トリガー周りとかはやらかすと怖いので分解はいったん保留して、ハンマー周りの不調を確認。
原因は思ったより簡単に見つかった。

f:id:doubtpad:20211025234912j:plain
ハンマー周り(通常状態)

基本状態はこういった形で、ハンマーを起こし始めるとハンマー側のツメがシリンダー・トリガー側のツメを引っかける形で引っ張ってくれる。

f:id:doubtpad:20211025235143j:plain
ハンマーを動かし始めると

f:id:doubtpad:20211025235233j:plain
前方のパーツを引っかけて動かす

この時に、ハンマーが後ろに行き過ぎてしまうと前パーツのツメとの噛み合わせが外れてしまい、ハンマーが戻らなくなる。
一応「カチッ」と音がして固定される位置があり、さらにグリップがハンマーに当たることで一定以上後ろに行かないようには対策されているので、そこまで頻発する故障ではなさそう。
ただし、それ以上行かないように固定される訳ではないので、グリップが緩んでいる状態でハンマーを引き過ぎると外れるっぽい。
確かに最初の確認でグリップにがたつきがあったので、もしかするとそれが悪さをしていたのかも。

f:id:doubtpad:20211025235445j:plain
噛み合わせが外れている

なお、トリガーを引くとバネの力でハンマーが戻る場合もあるのだが、ダメな時はダメっぽい。
まぁグリップの取り外しは簡単にできるので、原因が分かってさえいれば気にはならないレベル。

あと、フレーム側のネジも緩んでいた。
不調には影響してないと思うけど、気になるので締めなおしておく。

f:id:doubtpad:20211025234312j:plain
ここのネジが緩んでいた

たぶん直ったと思うので組み直す。
ついでに各部のクリーニングもやっておく。

f:id:doubtpad:20211026000925j:plain
修理後

f:id:doubtpad:20211026001207j:plain
シリンダー側から

f:id:doubtpad:20211026001232j:plain
スイングアウト

f:id:doubtpad:20211026001304j:plain
縦向き

実射してみたところ、ホップ強めだとそもそも弾が出なくて笑った。
ちょっと弱めた状態なら問題なく撃てたので、しっかり調整するようにした方が良さそう。


という訳で、さくっと修理できてしまったのでめちゃくちゃお得感あった。
強いて言えば、やっぱりグリップがガタついてるのが気になr

f:id:doubtpad:20211026001644j:plain
あっ(察し)

いや、なんとなくグリップ振ったらカタカタ言うから、覗き込んでみたら奥のこのパーツにがっつりヒビが入ってた。
このパーツを挟む形でグリップと本体をネジ止めしているので、これが無いとグリップが留められない。
ガタつきの原因が確定した代わりに、ガチの修理案件になった感。


とりあえず動作そのものは問題なくなったので良かった。
グリップ周りは今後修理するつもりでいるので、満足するレベルで仕上がったら書くかもしれない。

pngの保存形式が変わったと思ってたらggsave()の保存形式が変わっていた

去年の後半に作ったShinyアプリを久しぶりに動かしたら、処理の途中でいきなり落ちるとかいう恐ろしい目にあった。

具体的にはggplot2で作図したグラフをggsave()でpng保存し、それをEBImageパッケージを使って読込と解析をするアプリなのだが、その途中でアプリが落ちてしまった。
なお、一連の動作を全て纏めて走らせる形にしていたので、自作関数やらモジュール分割やらで原因まで辿り着くのに時間がかかって大変だった。

原因だが、直接の原因はpngのチャンネル数がアプリを組んだ時と変わっていたことで、以前は3チャンネル(RGB)だったのがいつの間にか透過レイヤー込みの4チャンネルになっていた。
で、その根本的な原因がggplot2パッケージに含まれているggsave()のアップデートだったっぽい。

チャンネル数の変化

チャンネル数の確認はEBImageパッケージ1で実施。
これに関しては見られれば何でもいい2

library(tidyverse)
library(EBImage)

g <- iris %>%
  ggplot(aes(Sepal.Length, Sepal.Width, colour=Species)) +
  geom_point() +
  theme_void()

ggsave(filename = "dat.png", plot = g)
readImage("dat.png") %>% channel("rgb")


以前の結果がこれ。
frames.totalがチャンネル数である。

> readImage("dat.png") %>% channel("rgb")
Image 
  colorMode    : Color 
  storage.mode : double 
  dim          : 640 360 3 
  frames.total : 3 
  frames.render: 1 

imageData(object)[1:5,1:6,1]
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    1    1    1    1    1
[2,]    1    1    1    1    1    1
[3,]    1    1    1    1    1    1
[4,]    1    1    1    1    1    1
[5,]    1    1    1    1    1    1


それがこうなっていた。

> readImage("dat.png") %>% channel("rgb")
Image 
  colorMode    : Color 
  storage.mode : double 
  dim          : 640 360 4 
  frames.total : 4 
  frames.render: 1 

imageData(object)[1:5,1:6,1]
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    1    1    1    1    1
[2,]    1    1    1    1    1    1
[3,]    1    1    1    1    1    1
[4,]    1    1    1    1    1    1
[5,]    1    1    1    1    1    1


Deviceの問題かと思ったのだが、RStudioのGlobal Option → General → Graphicsで設定を切り替えても変わらないので違うっぽい。
さらに、Rネイティブなpng()ではRGBの3チャンネルで問題なく出力できた。

ggsave("test_gg.png", g)
png("test_windows.png", type="windows"); g; dev.off()
png("test_cairo.png", type="cairo"); g; dev.off()
png("test_cairo-png.png", type="cairo-png"); g; dev.off()

c("test_gg.png", "test_windows.png", "test_cairo.png", "test_cairo-png.png") %>% 
  enframe(name=NULL) %>% 
  mutate(frame = value %>% map_dbl( ~.x %>% readImage() %>% numberOfFrames()) )
# A tibble: 4 x 2
  value              frame
  <chr>              <dbl>
1 test_gg.png            4
2 test_windows.png       3
3 test_cairo.png         3
4 test_cairo-png.png     3


ggsave()だけチャンネル数が増えているので、ggsave()かggplot2そのものが怪しい。
その前提でコードを眺めていて、何となくtheme_void()を外したら3チャンネルで出力された(出力結果は割愛)。
theme_void()はグラフの要素を全て消す設定だが、もしかして各要素を透過して見えないようにしているだけってこと……?

という訳で、theme_void()の中身やらggsave()のリファレンスやらを確認したところ、どうも背景色の設定が影響しているっぽい。

ggsave()で背景色を指定するbgオプションはデフォルト設定(bg=NULL)だとplot.backgroundのfillに指定した色を利用するのだが、theme_void()では色を設定しておらず3、その場合はtheme()の仕様としてNULLになる4

よって、ggplot2側でplot.backgroundの設定をするか、ggsave()側でgb="white"などのような設定をすれば透過レイヤ―が無くなるため、チャンネル数が3つで確定する。

g <- iris %>%
  ggplot(aes(Sepal.Length, Sepal.Width, colour=Species)) +
  geom_point() +
  theme_void() +
  theme(plot.background = element_rect(fill="white"))

ggsave("test_gg_bg.png", g, bg="white")


挙動が変わった理由

そもそもggsave()自体は冒頭で言及したShinyアプリで以前から使っていたので、何かアップデートに絡んで挙動が変わったんだろうと予想。
その方向性で情報を漁ってみたのだが、2020年1月のHadley神によるggsave()での保存時に関する提案をしたのがたぶん出発点……だと思う。

github.com

その後、いくつかやり取りがあって、最終的に背景色をthemeから引っ張ってくる形でまとまったっぽい。

github.com

github.com

github.com

機能としては2021-06-16のggplot2 3.3.4でリリースされた(Fixesの8番目に記載有り)。
思っていたより最近のアプデだったらしい。

ggplot2.tidyverse.org

なお、嫌な予感がして仕事のログを確認したら、冒頭のアプリを動かしてエラーを起こしたのがまさに3.3.4のリリース日だった。
そんなタイムリーなことってある???

アプデが最近だったこともあって、ここ最近も関連する話題が出ている模様。
気になる人は追ってみても面白いかも。

github.com

github.com

Enjoy!!


  1. Bioconductorに含まれている画像解析用のパッケージ。日本語の記事はそこまで多くないけど、比較的使いやすい印象。

  2. 画像解析ならPython使えよってツッコミはさておき、よく見かけるのはimagerパッケージとか。あとopencvパッケージも個人的に気になっている。更新止まったかもなぁとか思ってたら最新のPublishedが2021-05-04に来ているので、ちょこちょこ進めているのかも知れない。

  3. というかほとんどのテーマでplot全体の背景色は設定していないっぽい。グラフエリアの背景色は大体設定されているので混乱しそうになる。

  4. 今回見ているtheme_void()やよく使われるtheme_gray()などのtheme*()関数については、theme()で設定できる全てのパラメータがNULLになっているtheme(theme_all_null)に対して必要なところだけ上書きするような形で設定しているため、こういう挙動になる。……ということを今回初めてtheme()関連の仕組みを調べて知って、軸の設定やら何やらを弄った後でtheme*()を使うとリセットされるのはそういうことだったんだなっていうド素人ムーブをかましてました。

RのバージョンアップをRからやろうとして失敗したからゴリ押した話

少し前に、メジャーアップデートとしてR 4.1.0が来ていた。
地味に分量が多くて把握しきれていないが、Rネイティブなパイプ「|>」の実装とか、関数の「(x) x+1」表記とかが目立つ変更点かと思う。

せっかくならメジャーアップデートは対応しておきたいので、RユーザーらしくRからアップデートする。
と言っても、installrパッケージの中にあるupdateR()を使うだけなので楽。

……楽なんだけど、ちょっと油断すると割と面倒な落とし穴があるので書いておく。


やり方

ただ入れるだけなら、単純にこれでいい。

install.packages("installr")
library(installr)
updateR()


ちなみにRStudioから動かすとエラーが出るらしいので、R本体から実施した。

ただし、単純に関数だけで走らせてしまうと今まで使用していたパッケージをコピーできないという結構エグイ問題がある。
updateR(copy_packages = TRUE)で実行すれば勝手にやってくれるので、ちゃんと指定してから走らせると安心。

また、オプションを指定し忘れたとしても丁寧にポップアップ表示で「古いRからパッケージ移す?(意訳)」みたいな感じで訊いてくれるので、落ち着いて「Yes」を選択すればいい。
そのあと、古い方のRにパッケージを残すかなど訊かれるので、必要に応じて選択。
あとは勝手にやってくれる。


本題

さて。

自分の環境で、オプションは設定せずにポップアップでパッケージのコピーを実施したら、何故か失敗した。
エラーが出ていたので確認1

f:id:doubtpad:20210627152459p:plain
「権限がありません(意訳)」

どうやらフォルダを作ろうとしたけど権限がないから弾かれたっぽい。
そんなことある……?

悔しいので、あくまでもRで対応する。


旧バージョンからのパッケージのコピー

コピーとは書いたが、正確には「旧バージョンの持っているパッケージの情報を拾ってきてインストールする」という形。

まずは旧バージョンに入れていたパッケージの情報を拾ってくる。
パッケージは基本的に「Program Files」→「R」の各バージョンのフォルダ直下の「library」フォルダに保存されていると思うが、場合によっては「ドキュメント」→「R」→「win-library」の各バージョン直下に置いてある2
この中のパッケージ名を纏めて拾ってしまえばいい。

install.packages( list.files( path = "C:/Program Files/R/R-4.0.3/library", full.names = F ) )


なお、一括読み込みするためのコードで初期インストールされている関数以外を使用すると、「Rを再起動(≠RStudioの再起動)してからインストールします」と出るが、言う通りに再起動するとパッケージの読み込みが解除されるためコードが走らなくなる。
なのでbase関数で実行する必要がある点は注意。

library(fs)

#これだと上手く走らない
xxx <- fs::dir_info( path = "C:/Users/name/documents/R/win-library/4.0.3", full.names = F )
install.packages(xxx)


パッケージの名前をメモ帳などにコピーしてからインストールしたいという奇特な人については、まず初めにパッケージ名を全て纏めて取得してからメモ帳か何かにコピペ。

library(tidyverse)

xxx <- list.files( path = "C:/Users/name/documents/R/win-library/4.0.3", full.names = F ) %>%
    str_flatten(collapse=",")


これを新しいバージョンのR上で、以下のように処理すれば完了する。
ただし、この方法でinstall.packages()に大量のパッケージ名を渡すとエラーを起こすことがあるため、場合によっては分割して渡すなど注意が必要。

xxx <- "..." #メモ帳に保存した文字列をコピペ
xxx1 <- unlist(strsplit(xxx, ",")) #分割してunlist()
install.packages(xxx1)


この手法については、自分のRで入っているパッケージを他のPCのRにインストールする場合にも地味に役立つため、知ってると活躍することがある、かもしれない。

Enjoy!!


  1. ちなみにこのエラー、最後のポップアップで「今開いてるRは閉じる?(意訳)」でYesを押していると見れなくなるのが最大のトラップだと思う。

  2. 管理者権限がないPCか、管理者として起動せずにRを使っていた場合だと有り得る。

data.tableにfurrr::future_map()を入れるとむしろ速度が低下するっぽい

data.tableに対する処理の一部が妙に重くて困ってたのだが、その原因を調べてたらfurrr::future_map()との食い合わせが最悪だったってことが分かったのでメモ。


data.tableとは

data.frameを発展させたデータ形式で、高速かつメモリ効率の良さが売り。

cran.r-project.org

このdata.table形式で解析をしている場合、こういう書き方がよく出てくる(例としてirisで書いてみる)。

#花弁サイズの単位を cm → mm に変更
data.table(iris)[, lapply(.SD, function(x) x*10), ,
                 .SDcols=Petal.Width:Petal.Length]


この時のlapply()について、無名関数1を使う場合はmap()に置き換えている人も多い。

data.table(iris)[, lapply(.SD, ~.x*10), ,
                 .SDcols=Petal.Width:Petal.Length]


furrr::future_map()とは

tidyverseユーザーにはお馴染みのpurrr::map()系の関数について、置き換えるだけでお手軽に並列処理を実装できる関数。
様々なmap()系の関数に対応していて、しかも「future_」を頭につけるだけで使えるので非常に便利。

cran.r-project.org


何が問題だったの?

「大きめのデータの解析をするつもりで」「並列処理を入れたかった」ので、こんな感じの処理を書いていた(本来はもう少し複雑)。

data.table(iris)[, future_map(.SD, ~.*10), ,
                 .SDcols=Petal.Width:Petal.Length]


が、どうにも処理が重い。というか時間がかなりかかる。

処理に条件分岐とかを入れていたので、その辺がおかしいのだろうと勘違いしてずっと調べていたのだが、future_map()で書いていたところをmap()と書き間違えた時にいきなり爆速になって気が付いた


とりあえずベンチマークを取る

上で書いた処理3種類について、microbenchmarkパッケージで処理速度を比較してみる。

library(microbenchmark)

#処理は上で書いていたものをそのまま使用
microbenchmark(
  "lapply" = data.table(iris)[, lapply(.SD, function(x) x*10), ,
                              .SDcols=Petal.Width:Petal.Length],
  "map" = data.table(iris)[, map(.SD, ~.*10), ,
                           .SDcols=Petal.Width:Petal.Length],
  "future_map" = data.table(iris)[, future_map(.SD, ~.*10), ,
                                  .SDcols=Petal.Width:Petal.Length],
  times = 100) -> res
res %>% autoplot()
res %>% print()


結果:

f:id:doubtpad:20210514005224p:plain

## Unit: microseconds
##        expr    min       lq      mean   median       uq     max neval
##      lapply  750.3   848.80  1528.067  1110.35  1500.30 13169.6   100
##         map  792.3   879.75  1210.073   992.10  1264.75  4440.0   100
##  future_map 9952.7 12710.25 15790.534 14612.25 17863.30 56270.4   100

lapply()とmap()は同程度(若干map()の方が有利)なので好みで選んで良さそうだが、future_map()を使うと中央値で約15倍、平均でも10倍くらい遅い。
これだけ遅ければ、処理のボトルネックになるに決まってる。

原因

今さら思い出してみると「fread()は最初から並列化されてる」とかって話を見かけた記憶があったので、もしかしてdata.tableって全体的にそういう感じになってる……?

という訳で、公式ドキュメントやら何やら調べてみたところ、下記の資料が分かりやすかった。

www.john-ros.com

「16.4.2 Parallel Data Munging with data.table」から引用。

We now recall it to emphasize that various operations in data.table are done in parallel, using OpenMP. 意訳:もっかい言うけど2、data.tableはOpenMPを使って並列処理してるよ。

やはり内部で並列化していた。

というか読んでて気が付いたのだが、パッケージの関数にgetDTthreads()とsetDTthreads()という、そのものずばりな関数が入っていた。
やはりパッケージの関数は全部確認するべきだわ(遠い目)。

ちなみに、getDTthreads()を使うと現状のスレッド数の確認、setDTthreads()ではスレッド数の設定が可能。
ただし、スレッド数を設定可能な最大数にするとむしろ速度が低下する(上の資料でもそんな話が出ている)。
どうしても自力で設定しないといけない場合は最大スレッド数の半分とかにしておけば問題なさそうだけど、自分みたいな詳しくない人間は関数に任せた方が安心。

もう少し詳しい話が知りたい場合はこの辺も参考になりそう。
なお自分は斜め読みしか出来てないので今度読み直す。

github.com


たぶん今回の問題は、並列化している処理の内部で、さらに並列化しようとしたことで訳の分からない状態になっていたんじゃなかろうかと思う。
foreachとかでも同じようなことをしたら(試してはいないけど)遅くなる気がするので、data.tableの時はそのまま素直に処理を書いた方が楽で良さそう。

Enjoy!!


  1. 「~abs(fft(.x))2」みたいに書くアレ。Rを使い始めて日が浅いと、使い方を調べるのに名前が分からなくて意外と手間取ると思うの自分だけだろうか。

  2. なお、ここより前の本文中でOpenMPに関する記述はない(16.3.2でOpenMPI関連の話をしているので、それを指しているのかも)。

「MAD」が問題を抱えすぎている気がする

一昔前に某動画サイトでうpされていたアレの話とかではなくて。

データの標準化の話を調べていた時、「中央絶対偏差(MAD)」に行き当たったのが事の発端。

doubtpad.hatenablog.com

このMAD、バラつきに関する指標の1つで、標準偏差と比較すると外れ値に強い(ロバストな)性質がある。計算式はこんな感じ。

 \displaystyle
MAD = median|X_i-\overline{X}|

これ自体には何も問題はない。
問題なのは、同じくバラつきを表現する指標にMADと呼ばれるものが存在することである。

どういうことかというと、バラつきの指標には「Median Absolute Deviation(中央絶対偏差)」と「Mean Absolute Deviation(平均絶対偏差)」が存在しており、そのどちらもMADと略されるということである。ややこしいなオイ。

日本語でならまだ判断が付くものの、「MAD」だけでは他に情報がないかぎり判定は不可能。
ちなみに平均絶対偏差はこんな感じで算出する。

 \displaystyle
MAD = \frac{1}{n}\sum_{i=1}^{n}|x_i-m(X)|

ご覧の通り、まったくの別物。
訳分からなくなるんだから、「Average Absolute Deviation」とかにすればいいじゃんと思う。



en.wikipedia.org

あるのかよ。
とりあえず内容確認(概要から引用)。

The average absolute deviation (AAD) of a data set is the average of the absolute deviations from a central point. It is a summary statistic of statistical dispersion or variability.

訳:AADは中心点からの絶対偏差の平均値で、分散の要約統計量です。

なるほど、平均絶対偏差だな。

In the general form, the central point can be a mean, median, mode, or the result of any other measure of central tendency or any reference value related to the given data set.

訳:中心点としては平均値や中央値、最頻値、その他にもデータセットから得られる中心なら何でも使えます。

なるほど、平均絶対偏差だな。

AAD includes the mean absolute deviation and the median absolute deviation (both abbreviated as MAD).

訳:AADはMADの両方を含みます。

なるほど分からん。
平均だと思って読み進めてたからマジで混乱したわ。

なお、直後の「Measures of dispersion」の項目では「平均値周りの平均絶対偏差と中央値周りの中央絶対偏差は両方ともMADで混乱する上に計算結果も全然違う」「偏差の算出方法と中心点の設定、この2つの組み合わせが複数あるからAADは一意に計算値を表すものではない」と解説されている。

つまり、ここでの「Average」とは「中心値」を示しており、AADについては「中心絶対偏差」とでも訳すのが比較するのにちょうどいいと思われる。
すると、2種類のMADが含まれるのも納得できる1

また、この2つの呼び分けについてだが、個人的にはIBMこのページのように「MAD」と「MeanAD」のような形で使い分けると、後で見直した時の混乱を最小限に抑えられて良さそう。

補足だが、この項目より、最後の処理として「平均値」と「中央値」(=MeanとMedian)のどちらを採用するかによって、どちらのMADになるか決まるとある。
どういうことかというと、中心点の設定方法は偏差の呼び方に影響しない
そのため、どういう処理を行っているのかについては解析内容をしっかりと確認する必要がある。


ちなみに、Rにはmad()という関数が最初から入っている。
これで得られる値は中央絶対偏差なのだが、「標準偏差の推定量としてのMAD」を計算しているため、

 \displaystyle
MAD = median|X_i-\overline{X}| \\
k = 1/qnorm(0.75) \approx 1.4826

を用いて、

 \displaystyle
\hat\sigma = k ・ MAD

である \hat\sigma を計算結果として返してくる2

また、median()とabs()を用いてMADの計算を自前で実装した場合の値と、mad()のデフォルトでの返り値を比較すると、近似地である k = 1.4826 を用いて計算していることが確認できる3

# MADの算出
iris[[1]] %>% mad()

# 値を比較すると、近似値を用いて計算した場合と一致している
iris[[1]] %>% { median( abs(. - median(.)) ) / qnorm(0.75) }
iris[[1]] %>% { median( abs(. - median(.)) ) * 1.4826 }


以上、あまり出番は無いかも知れないけど、ロバストな指標を用いている場面で「MAD」と見えたら、「中心値」と「中心点」に注意した方がいいというお話しでした。

Enjoy!!



~おまけ~

ここまでの話の流れがあった上で日本語版のWikipediaを見てみる。
が、そもそも日本語版だと平均絶対偏差などは単独での記載ページはない

で、一応は「偏差」の中の項目として記載されているのだが。

ja.wikipedia.org

そもそも中央絶対偏差については記述自体が存在しない。
サイト内検索をしてみると、「対数コーシー分布」のページのパラメータの推定で言及があるくらい。

かなりしっかりと記載されている英語版とは情報量が桁違いなので、日本語で読みたい場合はおとなしくWikipedia以外で情報収集した方がいい。


  1. ちなみに、AADのWikiをスクロールしていくと下の方に「Maximum absolute deviation」という項目を見つけて身構えてしまうが、厳密にはバラつきの指標ではないとのこと。

  2. 厳密(?)な意味でのMADが必要な場合はパラメータで「constant=1」を指定すれば得られる。

  3. mad()のパラメータに「constant=1.4826」とあるので、ぶっちゃけ敢えて検証するほどのことではないのだが、確認しないと何となく落ち着かなかった。