中央絶対偏差と分位数：MATLAB vs R

はじめに

「あれ、、結果がぜんぜん合わない。そんな大胆なバグある？」

データサイエンスのための統計学入門という O'Reilly の本の内容を MATLAB で試しながら読んでいてドキッとした点を備忘録として書いておきます。

１章「探索的データ分析」を MATLAB で実装したものはこちら（GitHub: PSFDS-MATLAB）に置いてあります。

中央絶対偏差と分位数：MATLAB vs R

「1.4.3 例：州別人口の散らばりの推定」に出てくる中央絶対偏差や分位数を求める計算。MATLAB、R 双方に同じような名前の関数がありますが、デフォルトの状態（何も指定しないと）計算内容が異なります。

中央絶対偏差 mad

参考ページ

• R: https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/mad
• MATLAB: https://jp.mathworks.com/help/stats/mad.html

MATLAB の mad は中央絶対偏差・平均絶対偏差の両方を兼ねていますが、R は中央絶対偏差のみ。

y = mad(X,flag) は、平均絶対偏差 (flag = 0、既定) と中央絶対偏差 (flag = 1) のどちらを計算するかを指定します。

とある通り、MATLAB の mad 関数は何も指定しない場合は平均絶対偏差 mean(abs(X – mean(X))) を返し、flag = 1 とすると中央絶対偏差 median(abs(X – median(X))) を返します。

R の mad は中央絶対偏差を計算しますが、デフォルトで constant = 1.4826 が掛かった値が返ってくる点にも要注意です。

RCode(Display)

> x=c(3,5,8,10,12,12,15,17,19,20)
> mad(x)
[1] 6.6717

MATLAB の場合は以下の通りです。

Code

x = [3,5,8,10,12,12,15,17,19,20];
mad(x,1)

Output

ans = 4.5000

Code

constant = 1.4826;
constant*mad(x,1)

Output

ans = 6.6717

分位数 quantile

参照ページ

• R: https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/quantile
• MATLAB: https://jp.mathworks.com/help/stats/quantile.html

デフォルトで使用されるアルゴリズムが R と MATLAB で異なりますので注意が必要です。結論から言うと R での type = 5 （既定では type = 7）が MATLAB のデフォルトと一致します。type 別のアルゴリズムの挙動については上の公式ページか「分位数 パーセンタイル」（http://kj01.kgu.mydns.jp/fujimoto/tips/qptile.html）が参考になりました。

分位数の計算方法についてはMATLAB のドキュメンテーションページ「分位数と百分位数」も参考になります。

R での実行結果：

RCode(Display)

> quantile(x=c(3,5,8,10,12,12,15,17,19,20),probs=c(0,0.25,0.50,0.75,1),type=5)
  0%  25%  50%  75% 100% 
   3    8   12   17   20 

MATLAB だと以下の通り。

Code

x = [3,5,8,10,12,12,15,17,19,20];
p = 0:0.25:1;
quantile(x,p)

Output

ans = 1x5    
     3     8    12    17    20

四分位数間範囲 iqr

参照ページ

• R: https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/IQR
• MATLAB: https://jp.mathworks.com/help/stats/prob.normaldistribution.iqr.html

上と同じ理由で IQR の結果も異なります。type = 5 設定で R で実行すると以下の通り。

Code(Display)

> IQR(x=c(3,5,8,10,12,12,15,17,19,20),type=5)
[1] 9

そして MATLAB の実行結果と一致します。

Code

x = [3,5,8,10,12,12,15,17,19,20];
iqr(x)

Output

ans = 9

まとめ

分位数を計算する方法も沢山あるもんだなぁ・・と勉強になりました。
R は普段使わないので勘違いがあるかもしれません。間違いありましたらご指摘ください！コメントお待ちしています。