7

クラスを使用してscipy.stats.gaussian_kde、緯度と経度の情報で収集された離散データを滑らかにしようとしているので、最終的には、高密度がピークで低密度が谷である等高線図に似たものとして表示されます。

2 次元のデータセットをgaussian_kdeクラスに入れるのに苦労しています。私はそれが1次元データでどのように機能するかを理解するために遊んだので、2次元は次のようなものになると思いました:

from scipy import stats
from numpy import array
data = array([[1.1, 1.1],
              [1.2, 1.2],
              [1.3, 1.3]])
kde = stats.gaussian_kde(data)
kde.evaluate([1,2,3],[1,2,3])

これは、 に 3 点あることを意味し[1.1, 1.1], [1.2, 1.2], [1.3, 1.3]ます。x軸とy軸の幅1を使用して、1から3までを使用してカーネル密度を推定したいと考えています。

gaussian_kde を作成すると、次のエラーが表示され続けます。

raise LinAlgError("singular matrix")
numpy.linalg.linalg.LinAlgError: singular matrix

のソース コードをgaussian_kde調べてみると、データセットの意味と次元の計算方法がまったく異なることに気付きましたが、モジュールで多次元データがどのように機能するかを示すサンプル コードは見つかりませんでした。gaussian_kde多次元データを使用するいくつかのサンプル方法を教えてもらえますか?

4

4 に答える 4

7

この例はあなたが探しているもののようです：

import numpy as np
import scipy.stats as stats
from matplotlib.pyplot import imshow

# Create some dummy data
rvs = np.append(stats.norm.rvs(loc=2,scale=1,size=(2000,1)),
                stats.norm.rvs(loc=0,scale=3,size=(2000,1)),
                axis=1)

kde = stats.kde.gaussian_kde(rvs.T)

# Regular grid to evaluate kde upon
x_flat = np.r_[rvs[:,0].min():rvs[:,0].max():128j]
y_flat = np.r_[rvs[:,1].min():rvs[:,1].max():128j]
x,y = np.meshgrid(x_flat,y_flat)
grid_coords = np.append(x.reshape(-1,1),y.reshape(-1,1),axis=1)

z = kde(grid_coords.T)
z = z.reshape(128,128)

imshow(z,aspect=x_flat.ptp()/y_flat.ptp())

ここに画像の説明を入力してください

明らかに、軸は修正する必要があります。

を使用してデータの散布図を作成することもできます

scatter(rvs[:,0],rvs[:,1])

ここに画像の説明を入力してください

于 2011-05-25T14:55:52.200 に答える
4

カーネル密度推定と補間またはカーネル回帰を混同していると思います。ポイントのサンプルが多い場合、KDE ​​はポイントの分布を推定します。

どの補間が必要かはわかりませんが、scipy.interpolate のスプラインまたは rbf のいずれかがより適切になります。

1 次元のカーネル回帰が必要な場合は、scikits.statsmodels でいくつかの異なるカーネルを含むバージョンを見つけることができます。

更新：これが例です（これが必要な場合）

>>> data = 2 + 2*np.random.randn(2, 100)
>>> kde = stats.gaussian_kde(data)
>>> kde.evaluate(np.array([[1,2,3],[1,2,3]]))
array([ 0.02573917,  0.02470436,  0.03084282])

gaussian_kde は行に変数、列に観測値を持っているため、通常の統計とは方向が逆になっています。あなたの例では、3 つの点すべてが線上にあるため、完全な相関関係があります。それが特異行列の理由だと思います。

配列の向きを調整して小さなノイズを追加すると、例は機能しますが、(3,3) の近くにサンプル ポイントがないなど、非常に集中しているように見えます。

>>> data = np.array([[1.1, 1.1],
              [1.2, 1.2],
              [1.3, 1.3]]).T
>>> data = data + 0.01*np.random.randn(2,3)
>>> kde = stats.gaussian_kde(data)
>>> kde.evaluate(np.array([[1,2,3],[1,2,3]]))
array([  7.70204299e+000,   1.96813149e-044,   1.45796523e-251])
于 2010-11-09T00:28:43.233 に答える
1

gaussian_kdeSciPy マニュアルの2D データの操作方法の説明を理解するのが難しいことがわかりました。@endolith の例を補足するための説明を次に示します。コードをいくつかのステップに分割し、直感的ではない部分を説明するコメントを付けました。

まず、インポート:

import numpy as np
import scipy.stats as st
from matplotlib.pyplot import imshow, show

いくつかのダミー データを作成します。これらは、"X" および "Y" ポイント座標の 1 次元配列です。

np.random.seed(142)  # for reproducibility
x = st.norm.rvs(loc=2, scale=1, size=2000)
y = st.norm.rvs(loc=0, scale=3, size=2000)

2-D 密度推定の場合、gaussian_kdeオブジェクトは "X" および "Y" データセットを含む 2 行の配列で初期化する必要があります。NumPy の用語では、「それらを垂直に積み重ねる」:

xy = np.vstack((x, y))

したがって、「X」データは最初の行xy[0,:]にあり、「Y」データは 2 行目xy[1,:]にありxy.shape、(2, 2000). gaussian_kdeオブジェクトを作成します。

dens = st.gaussian_kde(xy)

推定された 2 次元密度 PDF を 2 次元グリッドで評価します。NumPy でこのようなグリッドを作成する方法は複数あります。ここでは、 @endolith の方法とは異なる (ただし機能的には同等の) アプローチを示します。

gx, gy = np.mgrid[x.min():x.max():128j, y.min():y.max():128j]
gxy = np.dstack((gx, gy)) # shape is (128, 128, 2)

gxyは 3 次元配列で、 の[i,j]番目の要素にgxyは、対応する "X" 値と "Y" 値の 2 要素リストが含まれます:gxy[i, j]の値は[ gx[i], gy[j] ]です。

各 2 次元グリッド ポイントでdens()(または同じことを)呼び出す必要があります。dens.pdf()NumPy には、この目的のために非常に洗練された関数があります。

z = np.apply_along_axis(dens, 2, gxy)

つまり、callable dens(同様に可能性があります) は3 次元配列の (3 番目の軸) にdens.pdf沿って呼び出され、値は 2 次元配列として返されます。唯一の不具合は、 の形状が期待したものではなく、 になることです。ドキュメントには次のように記載されていることに注意してください。axis=2gxyz(128,128,1)(128,128)

out [戻り値、LD] の形状は、軸方向の次元を除いて、arr の形状と同じです。この軸は削除され、func1d の戻り値の形状に等しい新しい次元に置き換えられます。したがって、func1d がスカラーを返す場合、out は arr よりも 1 つ少ない次元になります。

ほとんどdens()の場合、私が望んでいたスカラーではなく、長さ 1 のタプルが返されました。これは簡単に修正できるため、これ以上問題を調査しませんでした。

z = z.reshape(128, 128)

その後、画像を生成できます。

imshow(z, aspect=gx.ptp() / gy.ptp())
show()  # needed if you try this in PyCharm

これが画像です。( @endolith のバージョンも実装しており、これと見分けがつかない画像が得られたことに注意してください。)

上記のコマンドの出力

于 2020-10-23T12:13:17.690 に答える