ヒストグラムからの乱数

Question

scipy/numpyを使用してヒストグラムを作成すると、2つの配列があります.1つはビン数、もう1つはビンエッジです。ヒストグラムを使用して確率分布関数を表すと、その分布から乱数を効率的に生成するにはどうすればよいですか？

Answer 1

それはnp.random.choiceがオピーオーンの答え@で何をするか、おそらくですが、あなたは一様乱数に基づいて選択し、その後、正規化された累積密度関数を構築することができます。

from __future__ import division 
import numpy as np 
import matplotlib.pyplot as plt 

data = np.random.normal(size=1000) 
hist, bins = np.histogram(data, bins=50) 

bin_midpoints = bins[:-1] + np.diff(bins)/2 
cdf = np.cumsum(hist) 
cdf = cdf/cdf[-1] 
values = np.random.rand(10000) 
value_bins = np.searchsorted(cdf, values) 
random_from_cdf = bin_midpoints[value_bins] 

plt.subplot(121) 
plt.hist(data, 50) 
plt.subplot(122) 
plt.hist(random_from_cdf, 50) 
plt.show() 

---

2Dケースは、次のように実行できます。

data = np.column_stack((np.random.normal(scale=10, size=1000), 
         np.random.normal(scale=20, size=1000))) 
x, y = data.T       
hist, x_bins, y_bins = np.histogram2d(x, y, bins=(50, 50)) 
x_bin_midpoints = x_bins[:-1] + np.diff(x_bins)/2 
y_bin_midpoints = y_bins[:-1] + np.diff(y_bins)/2 
cdf = np.cumsum(hist.ravel()) 
cdf = cdf/cdf[-1] 

values = np.random.rand(10000) 
value_bins = np.searchsorted(cdf, values) 
x_idx, y_idx = np.unravel_index(value_bins, 
           (len(x_bin_midpoints), 
           len(y_bin_midpoints))) 
random_from_cdf = np.column_stack((x_bin_midpoints[x_idx], 
            y_bin_midpoints[y_idx])) 
new_x, new_y = random_from_cdf.T 

plt.subplot(121, aspect='equal') 
plt.hist2d(x, y, bins=(50, 50)) 
plt.subplot(122, aspect='equal') 
plt.hist2d(new_x, new_y, bins=(50, 50)) 
plt.show() 

Answer 2

おそらくこのようなものです。ヒストグラムのカウントをウェイトとして使用し、このウェイトに基づいてインデックスの値を選択します。

import numpy as np 

initial=np.random.rand(1000) 
values,indices=np.histogram(initial,bins=20) 
values=values.astype(np.float32) 
weights=values/np.sum(values) 

#Below, 5 is the dimension of the returned array. 
new_random=np.random.choice(indices[1:],5,p=weights) 
print new_random 

#[ 0.55141614 0.30226256 0.25243184 0.90023117 0.55141614] 

Answer 3

@Jaimeソリューションは素晴らしいですが、あなたは、ヒストグラムのKDE（カーネル密度推定）を使用して検討すべきです。ヒストグラムよりも統計を行うのがなぜ難しいのか、なぜkdeを代わりに使うのかがわかりますhere

scipyからkdeを使用する方法を示す@Jaimeのコードを編集しました。それはほぼ同じように見えますが、ヒストグラムジェネレータをよりよくキャプチャします。

from __future__ import division 
import numpy as np 
import matplotlib.pyplot as plt 
from scipy.stats import gaussian_kde 

def run(): 
    data = np.random.normal(size=1000) 
    hist, bins = np.histogram(data, bins=50) 

    x_grid = np.linspace(min(data), max(data), 1000) 
    kdepdf = kde(data, x_grid, bandwidth=0.1) 
    random_from_kde = generate_rand_from_pdf(kdepdf, x_grid) 

    bin_midpoints = bins[:-1] + np.diff(bins)/2 
    random_from_cdf = generate_rand_from_pdf(hist, bin_midpoints) 

    plt.subplot(121) 
    plt.hist(data, 50, normed=True, alpha=0.5, label='hist') 
    plt.plot(x_grid, kdepdf, color='r', alpha=0.5, lw=3, label='kde') 
    plt.legend() 
    plt.subplot(122) 
    plt.hist(random_from_cdf, 50, alpha=0.5, label='from hist') 
    plt.hist(random_from_kde, 50, alpha=0.5, label='from kde') 
    plt.legend() 
    plt.show() 


def kde(x, x_grid, bandwidth=0.2, **kwargs): 
    """Kernel Density Estimation with Scipy""" 
    kde = gaussian_kde(x, bw_method=bandwidth/x.std(ddof=1), **kwargs) 
    return kde.evaluate(x_grid) 


def generate_rand_from_pdf(pdf, x_grid): 
    cdf = np.cumsum(pdf) 
    cdf = cdf/cdf[-1] 
    values = np.random.rand(1000) 
    value_bins = np.searchsorted(cdf, values) 
    random_from_cdf = x_grid[value_bins] 
    return random_from_cdf 

Answer 4

私はOPと同じ問題を抱えていたし、私はこの問題への私のアプローチを共有したいと思います。

Jaime answerとNoam Peled answer私はKernel Density Estimation (KDE)を使用して2D問題の解決法を構築しました。

Frist、いくつかのランダムなデータを生成してから、そのProbability Density Function (PDF)をKDEから計算してみましょう。私はそれにexample available in SciPyを使用します。

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

def measure(n): 
    "Measurement model, return two coupled measurements." 
    m1 = np.random.normal(size=n) 
    m2 = np.random.normal(scale=0.5, size=n) 
    return m1+m2, m1-m2 

m1, m2 = measure(2000) 
xmin = m1.min() 
xmax = m1.max() 
ymin = m2.min() 
ymax = m2.max() 

X, Y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j] 
positions = np.vstack([X.ravel(), Y.ravel()]) 
values = np.vstack([m1, m2]) 
kernel = stats.gaussian_kde(values) 
Z = np.reshape(kernel(positions).T, X.shape) 

fig, ax = plt.subplots() 
ax.imshow(np.rot90(Z), cmap=plt.cm.gist_earth_r, 
      extent=[xmin, xmax, ymin, ymax]) 
ax.plot(m1, m2, 'k.', markersize=2) 
ax.set_xlim([xmin, xmax]) 
ax.set_ylim([ymin, ymax]) 

、プロットは次のとおりです。

今

、我々は変数ZあるKDE、から入手したPDFからランダムデータを取得します。

# Generate the bins for each axis 
x_bins = np.linspace(xmin, xmax, Z.shape[0]+1) 
y_bins = np.linspace(ymin, ymax, Z.shape[1]+1) 

# Find the middle point for each bin 
x_bin_midpoints = x_bins[:-1] + np.diff(x_bins)/2 
y_bin_midpoints = y_bins[:-1] + np.diff(y_bins)/2 

# Calculate the Cumulative Distribution Function(CDF)from the PDF 
cdf = np.cumsum(Z.ravel()) 
cdf = cdf/cdf[-1] # Normalização 

# Create random data 
values = np.random.rand(10000) 

# Find the data position 
value_bins = np.searchsorted(cdf, values) 
x_idx, y_idx = np.unravel_index(value_bins, 
           (len(x_bin_midpoints), 
           len(y_bin_midpoints))) 

# Create the new data 
new_data = np.column_stack((x_bin_midpoints[x_idx], 
          y_bin_midpoints[y_idx])) 
new_x, new_y = new_data.T 

そして、この新しいデータとプロットからKDEを計算することができます。

kernel = stats.gaussian_kde(new_data.T) 
new_Z = np.reshape(kernel(positions).T, X.shape) 

fig, ax = plt.subplots() 
ax.imshow(np.rot90(new_Z), cmap=plt.cm.gist_earth_r, 
      extent=[xmin, xmax, ymin, ymax]) 
ax.plot(new_x, new_y, 'k.', markersize=2) 
ax.set_xlim([xmin, xmax]) 
ax.set_ylim([ymin, ymax])