Как я могу извлечь контейнеры из объекта дистрибуции KDE? - PullRequest
Новая
       40

Как я могу извлечь контейнеры из объекта дистрибуции KDE?

2 голосов
/ 03 марта 2020

Фон

Я пытаюсь сгруппировать ~ 11 000 генов вида растения ( Arabidopsis thaliana ), по имеющимся у меня данным об их изменении уровней экспрессии в ответ на воздействие света.

Необработанные значения для каждого гена являются непрерывными случайными переменными, но я хочу sh, чтобы дискретизировать значения, скажем, для 20 дискретных классов.

Итак, вместо: 

change_in_expression = array([-2.2, -1.1, -1.2, ...,  0.6, -1. , -0.9])

У меня есть выходы класса: 

change_in_expression = array(["3_to_4","2_to_3","1_to_2",...])

Что я пробовал

Я строю распределение с seaborn's distplot(), который, я полагаю, использует KDE : 

import seaborn as sns

d = array([-2.2, -1.1, -1.2, ...,  0.6, -1. , -0.9]) # = change_in_expression

dis = sns.distplot(d, fit=stats.laplace, kde=False)

plt.title("Distribution of Differential Expression")
plt.xlabel("Log2FoldChange in expression")
plt.ylabel("Frequency")
plt.show()

enter image description here

И я знаю, что matplotlib.pyplot's hist() позволяет извлекать ячейки, когда настройки по умолчанию "автоматически" генерируют эти группировки ...

Резюме Вопрос

Вопрос в том, как я могу сгруппировать мои гены? Это более широкий вопрос, чем просто запрос seaborn версии matplotlib's hist() ..., поскольку seaborn's distplot использует KDE .

Не похоже, что я могу получить контейнеры из объекта ax, созданного seaborn, взглянув на методы, доступные из: 

dir(sns.distplot(d, fit=stats.laplace, kde=False)

Я думаю, один из способов был бы чтобы проверить внутренности distplot исходного кода Seaborn'а, выяснить, как они складывают данные перед построением графика ... но это далеко за пределы моего набора навыков единорога ...

1 Ответ

4 голосов
/ 03 марта 2020

Seaborn звонит pyplot.hist, что, в свою очередь, вызывает numpy.histogram. Таким образом, можно проверить, что seaborn использует в качестве аргумента bins, если ничего не указано. То есть a является данными, 

bins = min(_freedman_diaconis_bins(a), 50)

, где 1 

def _freedman_diaconis_bins(a):
    """Calculate number of hist bins using Freedman-Diaconis rule."""
    # From https://stats.stackexchange.com/questions/798/
    a = np.asarray(a)
    if len(a) < 2:
        return 1
    h = 2 * iqr(a) / (len(a) ** (1 / 3))
    # fall back to sqrt(a) bins if iqr is 0
    if h == 0:
        return int(np.sqrt(a.size))
    else:
        return int(np.ceil((a.max() - a.min()) / h))

и iqr 2 

def iqr(a):
    """Calculate the IQR for an array of numbers."""
    a = np.asarray(a)
    q1 = stats.scoreatpercentile(a, 25)
    q3 = stats.scoreatpercentile(a, 75)
    return q3 - q1

Все это должно быть примерно таким же, как 

bins = min(len(numpy.histogram_bin_edges(a, bins="fd")), 50)

или 

bins = 50 if len(numpy.histogram_bin_edges(a, bins="fd")) > 50 else "fd"

, затем передать bins в pyplot.hist, 

plt.hist(a, bins=bins)
...