How to obtain coordinates of maximum density

user_12274156 0
Viewed 1400

How could I get the coordinates of a point in the space with the greatest density.

enter image description here

I have this code to generate a random point and density analyze from this point.

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)

Thank you for help. :)

gaussian scipy matplotlib python-3.x python
1 Answers

For starters, let me state the obvious by saying that sns.jointplot computes the kernel density on its own, so your kernel variable is as of yet unused.

Here's what sns.jointplot generated for me with a random sample: fancy 2d kernel density contour with marginal distributions

There's a nice maximum at around (7, 5.4).

Here's what your kernel corresponds to:

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')

output from original kernel's KDE: scattered sharp peaks around input points

This will clearly not do: the density contains peaks centered around your input points; you will never be able to get a similar estimate than what sns.jointplot gave you.

We can easily fix this: you just have to drop the custom bw_method argument in the call to 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')

This looks just the way you want it: fixed contour plot which looks like the one from sns.jointplot

Now you know that this kernel.pdf is a bivariate function for which you're looking for the maximum.

And to find the maximum you should probably use something from scipy.optimize, for instance scipy.optimize.minimize (the trick is to look at the negative of your function, which turns maxima into minima).

Since your function will probably have a few local maxima, finding the global maximum reliably is not trivial. I would either use the aforementioned minimize, but first use a sparse mesh over the relevant domain and find the best maximum candidate first, or use a heavy-weight solver such as differential_evolution which is a stochastic solver that's supposed to be good at finding the true global minimum of a function.

Root finding and minimization is always fickle business, so you will have to play around with your real data and available methods to find a reliable workflow that gives you your maximum.

user_50673110
Related