Как получить координаты максимальной плотности

Question

Как я могу получить координаты точки в пространстве с наибольшей плотностью.

У меня есть этот код для генерации случайной точки и анализа плотности из этой точки.

import numpy as np
from scipy import stats
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

def random_data(N):
    # Generate some random data.
    return np.random.uniform(0., 10., N)

x_data = random_data(50)
y_data = random_data(50)

kernel = stats.gaussian_kde(np.vstack([x_data, y_data]), bw_method=0.05)

b = plt.plot(x_data, y_data, 'ro')
df = pd.DataFrame({"x":x_data,"y":y_data})
p = sns.jointplot(data=df,x='x', y='y',kind='kde')

plt.show(p)

Спасибо за помощь. :)

Answer 1

Для начала позвольте мне сформулировать очевидное, сказав, что sns.jointplot вычисляет плотность ядра самостоятельно, поэтому ваша переменная kernel еще не используется.

Вот что sns.jointplot сгенерировано дляя со случайной выборкой:

Есть хороший максимум около (7, 5,4).

Вот то, что соответствует вашему kernel:

x,y = np.mgrid[:10:100j, :10:100j]  # 100 x 100 grid for plotting
z = kernel.pdf(np.array([x.ravel(),y.ravel()])).reshape(x.shape)

fig,ax = plt.subplots()
ax.contourf(x, y, z, levels=10)
ax.axis('scaled')

Это явно не подходит: плотность содержит пики, центрированные вокруг ваших входных точек;вы никогда не сможете получить оценку, аналогичную той, которую дал вам sns.jointplot.

Мы легко можем это исправить: вам просто нужно сбросить пользовательский аргумент bw_method в вызове gaussian_kde:

kernel = stats.gaussian_kde(np.vstack([x_data, y_data]))

x,y = np.mgrid[:10:100j, :10:100j]  # 100 x 100 grid for plotting
z = kernel.pdf(np.array([x.ravel(),y.ravel()])).reshape(x.shape)

fig,ax = plt.subplots()
ax.contourf(x, y, z, levels=10)
ax.axis('scaled')

Это выглядит так, как вы этого хотите:

Теперь вы знаете, что эта kernel.pdf является двумерной функцией, для которой вы ищетедля максимума.

И чтобы найти максимум, вы, вероятно, должны использовать что-то из scipy.optimize, например scipy.optimize.minimize (хитрость заключается в том, чтобы посмотреть наотрицательно для вашей функции, которая превращает максимумы в минимумы).

Поскольку ваша функция, вероятно, будет иметь несколько локальных максимумов, надежное нахождение глобального максимума не является тривиальным. Я бы либо использовал вышеупомянутый minimize, но сначала использовал бы разреженную сетку над соответствующим доменом и сначала нашел бы наилучшего максимального кандидата, или использовал бы решатель с большим весом, такой как differential_evolution, который является стохастическим решателем, который должен хорошо находить истинный глобальный минимум функции.

Поиск корня и минимизация - это всегда непростое дело, поэтому вам придется поиграться с вашими реальными данными и доступнымиметоды, чтобы найти надежный рабочий процесс, который дает вам максимум.