Построение квантилей вокруг средней гистограммы

Question

У меня есть N образцов, в каждом из которых есть n значений. Для каждого образца я строю гистограмму из n значений. В общей сложности у меня была бы N гистограмма. Теперь я хотел бы иметь график, который показывает среднее значение этих гистограмм и имеет 5-95% квантильной области. Пример такого графика может выглядеть следующим образом (пожалуйста, не обращайте внимания на пунктирную линию и черную линию, только на затененную область.)

Итак пока что я построил все N гистограмм друг на друге.

Я понял, что среднее значение этих гистограмм будет гистограммой всех [NX n] значения вместе. Пример кода для этого будет выглядеть следующим образом:

import numpy as np
import matplotlib.pyplot as plt

samples = []
N = 20
n = 250

for i in range(N):
    samples.append(np.random.normal(loc=np.random.rand(1,)[0]/5-0.1, scale=1., size=n))

values_all = None
for i in range(len(samples)):

    values = samples[i]
    print(values)

    weights = np.ones_like(values) / float(len(values))
    plt.hist(values, range=[-4, 4], density=False, histtype='step', color='red', bins=15, weights=weights)

    if values_all is None:
        values_all = values
    else:
        values_all = np.concatenate(([values_all, values]), axis=0)

weights = np.ones_like(values_all) / float(len(values_all))
plt.hist(values_all, range=[-4, 4], density=False, histtype='step', color='black', bins=15, weights=weights)
plt.show()

Буду признателен за любые предложения о том, как найти и построить квантили 5-95%.

Answer 1

Вы должны получить вероятности для каждого бина и найти квантили для них. Вот пример кода

import numpy as np
import matplotlib.pyplot as plt

# generate data
samples = []
N = 20
n = 250

for i in range(N):
    samples.append(np.random.normal(loc=np.random.rand(1,)[0]/5-0.1, scale=1., size=n))


prob_all = None
for i in range(len(samples)):

    values = samples[i]

    weights = np.ones_like(values) / float(len(values))
    n, bins, patches = plt.hist(values, range=[-4, 4], density=False, histtype='step', color='red', bins=15, weights=weights, alpha = 0.5)

    # concatanate bin probabilities
    if prob_all is None:
        prob_all = n.reshape(-1,1)
    else:
        prob_all = np.concatenate(([prob_all, n.reshape(-1,1)]), axis=1)

    plt.close() # don't plot previous histograms

# find quantiles for each bin
quant = np.quantile(prob_all, [0.05, 0.5, 0.95], axis=1)

# plot histogram from bins and probabilities
def plt_hist(bins, quant, clr, alph, lw):
    for j in range(len(n)):
        plt.plot([bins[j], bins[j + 1]], [quant[j], quant[j]], color=clr, linewidth=lw, alpha = alph)

    plt.plot([bins[0], bins[0]], [0., quant[0]], color=clr, linewidth=lw, alpha = alph)
    plt.plot([bins[len(n)], bins[len(n)]], [quant[len(n) - 1], 0.], color=clr, linewidth=lw, alpha = alph)
    for j in range(len(n) - 1):
        plt.plot([bins[j + 1], bins[j + 1]], [quant[j], quant[j + 1]], color=clr, linewidth=lw, alpha = alph)

fig, ax = plt.subplots()

ax.set_ylim([0., 0.3])
ax.set_xlim([-4.5, 4.5])

# plot 50% quantile (mean)
plt_hist(bins, quant[1], clr='blue', alph=1., lw=1.)

# shade between quantiles
for i in range(len(n)):
    x = np.arange(bins[i], bins[i+1], 0.0001)
    y1 = quant[0,i]
    y2 = quant[2,i]

    ax.fill_between(x, y1, y2, facecolor='red', alpha=0.4)

# boarder for shadings
plt_hist(bins, quant[0], clr='black', alph=.2, lw=1.)
plt_hist(bins, quant[2], clr='black', alph=.2, lw=1.)

plt.show()