MATLABで95%等確率楕円/信頼楕円の作成方法についてご紹介します。実際のコードも掲載しているのですぐに利用可能です。
あくまで素人が作成したコードですので、ご利用は自己責任でお願いします。
95%等確率楕円/信頼楕円の作成方法
コードの例
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 | %% 95%等確率楕円分布 p=0.95; fp = {'FontName', 'Arial','FontWeight','bold'} ; %データ生成 location=-5+10*rand(20,2); %95%等確率楕円用パラメータ means=mean(location); covariance=cov(location); [V,D]=eig(covariance); if D(1,1)>D(2,2) A(1,1)=D(1,1); A(1,2)=D(2,2); B(:,1)=V(:,1); B(:,2)=V(:,2); else A(1,2)=D(1,1); A(1,1)=D(2,2); B(:,2)=V(:,1); B(:,1)=V(:,2); end c=sqrt(icdf('Chisquare',p,2)); w=2*c*sqrt(A(1,1)); h=2*c*sqrt(A(1,2)); theta=atan(covariance(1,2)/covariance(1,1)); a=w/2;b=h/2; t=-pi:0.1:pi+0.1; n=length(t); x0=means(:,1);y0=means(:,2); c=[x0*ones(1,n);y0*ones(1,n)]; x1=[a*cos(t); b*sin(t)]; R=[cos(theta) -sin(theta); sin(theta) cos(theta)]; X=R*x1+c; %グラフ作成 figure(); scatter(location(1:20,1), location(1:20,2),60, 'filled', 'MarkerFaceColor','k'); hold on; fill(X(1,:),X(2,:),'k','FaceAlpha',0.1, 'LineStyle','none'); hold on %長軸の描画/式の表示 x2=-10:0.1:10; y2=tan(theta)*x2+(means(:,2)-(tan(theta)*means(:,1))); plot(x2,y2,'k:'); hold on str = {sprintf('y = %.2fx + %.2f',tan(theta),means(:,2)-(tan(theta)*means(:,1)))}; text(3,8,str,fp{:},'Fontsize', 10,'HorizontalAlignment', 'center') %軸などの設定 xlim([-10 10]) ylim([-10 10]) xticks(-10:5:10) yticks(-10:5:10) xlabel('x', fp{:}, 'fontsize', 18) ylabel('y', fp{:}, 'fontsize', 18) ax = gca; ax.XAxis(1).Color = [0 0 0]; daspect([1 1 1]) set(gca,'FontWeight','bold'); set(gca,'FontSize',14); set(gca,'linewidth',1.5); box on; title('95%等確率楕円',fp{:},'fontsize', 20) %画像として保存 cd figure saveas(gcf,'95%ellipse','png') cd .. |
ざっくりと本コードの手順は次の通りです。
- 乱数を使用して20個の座標を作成
- 95%等確率楕円用パラメータを計算
- 散布図・楕円を表示
- 楕円の長軸の式を表示
- 標準偏差や軸の設定
- 画像として保存
95%等確率楕円の作成に重要なのが③です。
この部分で何を行なっているかざっくり説明していきます。
コードの解説
1 2 3 4 5 6 7 8 9 10 11 12 13 | %95%等確率楕円用パラメータ means=mean(location); covariance=cov(location); [V,D]=eig(covariance); if D(1,1)>D(2,2) A(1,1)=D(1,1); A(1,2)=D(2,2); B(:,1)=V(:,1); B(:,2)=V(:,2); else A(1,2)=D(1,1); A(1,1)=D(2,2); B(:,2)=V(:,1); B(:,1)=V(:,2); end |
全データの中心から楕円の中心座標を求めます。
そして、固有値・固有ベクトルについてそれぞれ求めていきます。
また、データが中心からどの程度離れているかを表すマハラノビス距離(D)を求めます。
このマハラノビス距離の2乗はχ^2分布（カイ二乗分布）するという性質を持ちます。
そのため、内部のカイ二乗分布の累積確率が95%となる楕円を求めると95%等確率楕円が描けます。
1 2 3 4 5 | c=sqrt(icdf('Chisquare',p,2)); w=2*c*sqrt(A(1,1)); h=2*c*sqrt(A(1,2)); theta=atan(covariance(1,2)/covariance(1,1)); |
逆累積分布関数(icdf)を使ってカイ二乗分布(‘Chisquare’)の自由度2、p=0.95となるマハラノビス距離を求めます。
マハラノビス距離と固有値を利用すれば、95%等確率楕円の長軸/短軸の長さが求まります。
また、(共分散÷xの分散)で軸の傾きも求めることができます。
回転行列で軸の傾き分を回転させた楕円を描画することで95%等確率楕円の完成です。
追記
(共分散÷xの分散)での軸の傾きだと上手くフィットしない時があるようです。
原因は不明ですが現在調査中であり、原因が分かり次第追記します。
1 | theta=atan(((A(1,1)-A(1,2))/covariance(1,2))); |
なお、どのような計算を行っているのか理解できていませんが上記のコードで傾きが計算できるようです。
図を見た感じだとこのコードの方がフィットしている感じがあります。
むらしん計算ミスの原因と本コードの計算原理が分かる方おられましたら、コメントにてご教授いただけると幸いです。
出力結果
上の図が、上記コードの実行結果です。
乱数で作成した20個の座標に対して、95%等確率楕円を表示することができています。
今回はrand関数を使用してデータをグラフ作成を行いました。
実際に使用する際には、locationの部分に自分のデータを当てはめてみてください。
なお、使用する際は自己責任でお願いします。
さいごに
以上、MATLABで95%等確率楕円/信頼楕円の作成方法です。
自分なりに試行錯誤して作成したコードですが、間違いを含んでいる可能性があります。間違いの指摘やコメントがありましたら本記事のコメント欄にお願いします。
最後までお読みいただきありがとうございました。