2017-01-10 30 views
4

カーネル密度推定の等高線図を描きたいと思います.KDEは、等高線プロットの塗りつぶし領域内に統合されています。Python - 等高線内の2Dカーネル密度推定の統合

例として、私は2DデータのKDEを計算想像:

data = np.random.multivariate_normal((0, 0), [[1, 1], [2, 0.7]], 100) 
x = data[:, 0] 
y = data[:, 1] 
xmin, xmax = min(x), max(x) 
ymin, ymax = min(y), max(y) 
xx, yy = np.mgrid[xmin:xmax:100j, ymin:ymax:100j] 
positions = np.vstack([xx.ravel(), yy.ravel()]) 
values = np.vstack([x, y]) 
kernel = st.gaussian_kde(values) 
f = np.reshape(kernel(positions).T, xx.shape) 

私はKDEのContourPlotはを描画する方法を知っています。

fig = plt.figure() 
ax = fig.gca() 
ax.set_xlim(xmin, xmax) 
ax.set_ylim(ymin, ymax) 
cfset = ax.contourf(xx, yy, f, cmap='Blues') 
cset = ax.contour(xx, yy, f, colors='k') 
plt.show() 

ただし、この等高線図は、塗りつぶされた各領域内の確率密度が何であるかを示しています。代わりに、私は塗りつぶし領域のそれぞれの中に入る総確率を示すプロットを望みます。

答えて

1

輪郭線が「単調」、つまり輪郭線内では、対応する輪郭レベルより上のピクセル値しか見つからない限り、次の点が正しいことに注意してください。また、密度がマルチピークになっている場合、別々のピークの対応する領域がまとめて表示されます。

これが真実/許容値である場合、ピクセルを値で並べ替えることで問題を解決できます。

私は「レベル」という変数にあなたが

ff = f.ravel() 
order = np.argsort(ff) 
fsorted = ff[order] 
F = np.cumsum(fsorted) 
# depending on how your density is normalised next line may be superfluous 
# also note that this is only correct for equal bins 
# and, finally, to be unimpeachably rigorous, this disregards the probability 
# mass outside the field of view, so it calculates probability condtional 
# on being in the field of view 
F /= F[-1] 
boundaries = fsorted.searchsorted(levels) 
new_levels = F[boundaries] 
ような何かを試みることができるあなたのプロットプログラムはその等高線のレベルを選択しますが、あなたがそれらを持っていると仮定すると(言って、昇順に)保存されたヒューリスティックれるかわかりません

これを使用できるようにするには、プロットプログラムで自由に輪郭ラベルを選択するか、少なくとも輪郭を配置するレベルを選択する必要があります。 1が実際に持っているしたい場合は、後者の場合には、kwarg「レベル」

# make a copy to avoid problems with in-place shuffling 
# i.e. overwriting positions whose original values are still to be read out 
F[order] = F.copy() 
F.shape = f.shape 
cset = ax.contour(xx, yy, F, levels=new_levels, colors='k') 

がありますと仮定すると、私は、最後に

それをより見やすくするために、あなたのコメントのいずれかから、次をコピーしましたcb = fig.colorbar(cfset、ax = ax)values = cb.values.copy()values [1:] - = values [: - 1] .copy(これは各塗りつぶし領域内の確率です。 )cb.set_ticklabels(値) - Laura

+0

ありがとうございます。答えは完全には機能しませんが、解決策に近づいています(少なくとも輪郭が単調であるこの特定のケースでは)。あなたはF [order] = 1 - Fを割り当てるときにあなたが何をするのかを明確にしてください。ありがとうございます – Laura

+0

私は問題を発見しました。この解は、ax.contourf(xx、yy、(1-F)、levels = new_levels)をプロットすると機能します。 – Laura

+0

素晴らしい!だからF [order] = 1-FをF [order] = Fで置き換えて、もう一方の行はそのままにしておきましょう。 F [order] = F行が実際に何をしているかについてのアドバイスが必要ですか? –

関連する問題