7

次の問題に対する迅速でエレガントな Python/Scipy/Numpy ソリューションを知っていますか: 関連する値 w (すべて 1D 配列) を持つ x、y 座標のセットがあります。次に、x と y を 2D グリッド (サイズ BINSxBINS) にビン化し、各ビンの w 値の分位数 (中央値など) を計算します。これにより、最終的に、必要な分位数を持つ BINSxBINS 2D 配列が得られます。

これはネストされたループで簡単に実行できますが、もっと洗練された解決策があると確信しています。

ありがとう、マーク

4

4 に答える 4

5

これは私が思いついたものです。役に立つことを願っています。必ずしもループを使用するよりもクリーンまたは優れているとは限りませんが、より良いものに向けて始めることができるかもしれません。

import numpy as np
bins_x, bins_y = 1., 1.
x = np.array([1,1,2,2,3,3,3])
y = np.array([1,1,2,2,3,3,3])
w = np.array([1,2,3,4,5,6,7], 'float')

# You can get a bin number for each point like this
x = (x // bins_x).astype('int')
y = (y // bins_y).astype('int')
shape = [x.max()+1, y.max()+1]
bin = np.ravel_multi_index([x, y], shape)

# You could get the mean by doing something like:
mean = np.bincount(bin, w) / np.bincount(bin)

# Median is a bit harder
order = bin.argsort()
bin = bin[order]
w = w[order]
edges = (bin[1:] != bin[:-1]).nonzero()[0] + 1
med_index = (np.r_[0, edges] + np.r_[edges, len(w)]) // 2
median = w[med_index]

# But that's not quite right, so maybe
median2 = [np.median(i) for i in np.split(w, edges)]

numpy.histogram2d も見てください

于 2012-04-24T22:39:16.720 に答える
1

あなたのコードに感謝します。それに基づいて、私の問題の次の解決策を見つけました（コードのわずかな変更のみ）：

import numpy as np
BINS=10
boxsize=10.0
bins_x, bins_y = boxsize/BINS, boxsize/BINS
x = np.array([0,0,0,1,1,1,2,2,2,3,3,3])
y = np.array([0,0,0,1,1,1,2,2,2,3,3,3])
w = np.array([0,1,2,0,1,2,0,1,2,0,1,2], 'float')

# You can get a bin number for each point like this
x = (x // bins_x).astype('int')
y = (y // bins_y).astype('int')
shape = [BINS, BINS]
bin = np.ravel_multi_index([x, y], shape)


# Median 
order = bin.argsort()
bin = bin[order]
w = w[order]
edges = (bin[1:] != bin[:-1]).nonzero()[0] + 1
median = [np.median(i) for i in np.split(w, edges)]

#construct BINSxBINS matrix with median values
binvals=np.unique(bin)
medvals=np.zeros([BINS*BINS])
medvals[binvals]=median
medvals=medvals.reshape([BINS,BINS])

print medvals
于 2012-04-25T21:28:57.317 に答える
0

numpy/scipy を使用すると、次のようになります。

    import numpy as np
    import scipy.stats as stats

    x = np.random.uniform(0,200,100)
    y = np.random.uniform(0,200,100)
    w = np.random.uniform(1,10,100)

    h = np.histogram2d(x,y,bins=[10,10], weights=w,range=[[0,200],[0,200]])
    hist, bins_x, bins_y = h
    q = stats.mstats.mquantiles(hist,prob=[0.25, 0.5, 0.75])

    >>> q.round(2)
    array([ 512.8 ,  555.41,  592.73])

    q1 = np.where(hist<q[0],1,0)
    q2 = np.where(np.logical_and(q[0]<=hist,hist<q[1]),2,0)
    q3 = np.where(np.logical_and(q[1]<=hist,hist<=q[2]),3,0)
    q4 = np.where(q[2]<hist,4,0)

    >>>q1 + q2 + q3 + q4
    array([[4, 3, 4, 3, 1, 1, 4, 3, 1, 2],
   [1, 1, 4, 4, 2, 3, 1, 3, 3, 3],
   [2, 3, 3, 2, 2, 2, 3, 2, 4, 2],
   [2, 2, 3, 3, 3, 1, 2, 2, 1, 4],
   [1, 3, 1, 4, 2, 1, 3, 1, 1, 3],
   [4, 2, 2, 1, 2, 1, 3, 2, 1, 1],
   [4, 1, 1, 3, 1, 3, 4, 3, 2, 1],
   [4, 3, 1, 4, 4, 4, 1, 1, 2, 4],
   [2, 4, 4, 4, 3, 4, 2, 2, 2, 4],
   [2, 2, 4, 4, 3, 3, 1, 3, 4, 4]])

prob = [0.25, 0.5, 0.75] は分位数設定のデフォルト値です。変更することも、そのままにしておくこともできます。

于 2012-05-21T23:07:47.583 に答える