【Illustrator JavaScript】カーネル密度推定の分布を描く

Adobe Illustrator で可視化シリーズです。
前回の記事では確率密度関数の区間に色を着けました。

今回は良く知られた関数形で与えられる確率密度関数ではなく、データから推定された分布を描きたいと思います。

1. 確率密度の推定

データから確率密度を推定する方法として、カーネル密度推定があります。
Python の Scipy.stats モジュールでカーネル密度推定ができるので、今回も Scipy の力に頼ります。

Python スクリプト

import numpy as np
from scipy import stats

# 平均と標準偏差
mu = 0
sigma = 1

# 正規乱数を発生
r_norm = stats.norm.rvs(loc=0, scale=1, size=1000)

# カーネル密度推定
norm_density = stats.kde.gaussian_kde(r_norm)


# 分布全体の確率密度を算出
prob_x = np.linspace(-4, 4, 100)
prob_y = norm_density(prob_x)
prob_xy = np.stack([prob_x, prob_y], axis=1)
np.savetxt('norm_prob.csv', prob_xy, delimiter=',')


# 2.5%点から97.5%点までの確率密度を算出
fill_x_min = np.percentile(r_norm, 2.5)
fill_x_max = np.percentile(r_norm, 97.5)
fill_x = np.linspace(fill_x_min, fill_x_max, 100)
fill_y = norm_density(fill_x)
fill_xy = np.stack([fill_x, fill_y], axis=1)

fill_xy = np.concatenate([
   np.array([[fill_xy[0, 0], 0]]),
   fill_xy,
   np.array([[fill_xy[-1, 0], 0]])
])

np.savetxt('norm_fill.csv', fill_xy, delimiter=',')

データとして正規乱数を 1000 個生成しました。ここは何らかの観測されたデータを使っても良いです。

2. 密度推定結果の描画

JavaScript

// 新規ドキュメント作成
app.documents.add(DocumentColorSpace.RGB, 800, 600);


// RGBカラーを設定し、カラーオブジェクトを返す
function setRGBColor(r, g, b) {
   var RGB = new RGBColor();
   RGB.red = r;
   RGB.green = g;
   RGB.blue = b;
   return RGB;
}


// x_margin
var x_margin = 400;

// y_margin
var y_margin = 200;

// 横方向の拡大率
var x_rate = 80;

// 縦方向の拡大率
var y_rate = 400;


// 塗り用テキストファイル読込
var fill_file = new File("C:\\Users\\mare_ism\\Documents\\norm_fill.csv");
fill_file.encoding = "UTF-8";
fill_file.open("r", "TEXT");
var fill_text = fill_file.read();

var fill_textRow = fill_text.split("\n");

// 確率密度関数の範囲を塗る
var fillObj = activeDocument.pathItems.add();
var fill_arrays = [];
for (var i = 0; i < fill_textRow.length; i++) {
   if (fill_textRow[i].length != 0) {
       var fill_rowItems = fill_textRow[i].split(",");
       var fill_array = [fill_rowItems[0] * x_rate + x_margin,
                         fill_rowItems[1] * y_rate + y_margin];
       fill_arrays.push(fill_array);
  }
}

fillObj.setEntirePath(fill_arrays);
fillObj.filled = true; // 塗りを「あり」にする
fillObj.fillColor = setRGBColor(105, 189, 131);
fillObj.stroked = false; // 線を表示しない


// x軸
var xAxisPos = [[-4, 0], [4, 0]];
var xAxisObj = activeDocument.pathItems.add();
xAxisObj.setEntirePath([
   [xAxisPos[0][0] * x_rate + x_margin, xAxisPos[0][1] * y_rate + y_margin],
   [xAxisPos[1][0] * x_rate + x_margin, xAxisPos[1][1] * y_rate + y_margin]
]);
xAxisObj.stroked = true;
xAxisObj.strokeWidth = 0.8;
xAxisObj.strokeColor = setRGBColor(0, 0, 0);

// y軸
var yAxisPos = [[-3.5, -0.1], [-3.5, 0.5]];
var yAxisObj = activeDocument.pathItems.add();
yAxisObj.setEntirePath([
   [yAxisPos[0][0] * x_rate + x_margin, yAxisPos[0][1] * y_rate + y_margin],
   [yAxisPos[1][0] * x_rate + x_margin, yAxisPos[1][1] * y_rate + y_margin]
]);
yAxisObj.stroked = true;
yAxisObj.strokeWidth = 0.8;
yAxisObj.strokeColor = setRGBColor(5, 5, 5);


// 曲線用テキストファイル読込
var prob_file = new File("C:\\Users\\mare_ism\\Documents\\norm_prob.csv");
prob_file.encoding = "UTF-8";
prob_file.open("r", "TEXT");
var prob_text = prob_file.read();

var prob_textRow = prob_text.split("\n");

// 確率密度関数の曲線を描画
var lineObj = activeDocument.pathItems.add();
var prob_arrays = [];
for (var i = 0; i < prob_textRow.length; i++) {
   if (prob_textRow[i].length != 0) {
       var prob_rowItems = prob_textRow[i].split(",");
       var prob_array = [prob_rowItems[0] * x_rate + x_margin,
                         prob_rowItems[1] * y_rate + y_margin];
       prob_arrays.push(prob_array);
  }
}

lineObj.setEntirePath(prob_arrays);
lineObj.filled = false; // 塗りを「なし」にする
lineObj.stroked = true; // 線を表示する
lineObj.strokeWidth = 2;
lineObj.strokeColor = setRGBColor(0, 135, 60);

こちらは前回のものと同じです。

実行結果

カーネル密度推定で得られた分布の 2.5% 点から 97.5% 点までの 95% の範囲に色を着けることができました。

事後分布の MCMC サンプルを使って、ベイズ信用区間を描画するのに使ったりできそうですね。

3. 注意事項

パーセンタイル値を算出するのに元のデータを使っていますが、これはカーネル密度推定で得られた分布のパーセンタイル値と一致する保証はありません。

サンプルサイズが大きければそれほど大きくはずれないでしょうが、サンプルサイズが小さい場合は注意が必要です。

これはサンプルサイズを 30 にして再実行した結果です。
左右の色が着いていない部分がおそらくは 2.5% 以上になっている結果が得られています。