Следующий код сравнивает вычисление объема пересечения либо с помощью scipy's dblquad, либо путем взятия среднего значения по сетке.
Примечания:
- Для 2D-случая (и только с 100 точками выборки) кажется, что дельта должна быть значительно больше 10%. В приведенном ниже коде используется 25%. При дельте 10% расчетные значения для
f1 и f2 составляют около 0.90, тогда как теоретически они должны быть 1.0. При дельте 25% эти значения составляют около 0.994. - Чтобы приблизить объем простым способом, среднее значение необходимо умножить на площадь (здесь
(xmax - xmin)*(ymax - ymin)). Кроме того, чем больше точек сетки учитывается, тем лучше аппроксимация. В приведенном ниже коде используются точки сетки 1000x1000. - В Scipy есть несколько специальных функций для вычисления интеграла, например
scipy.integrate.dblquad. Это намного медленнее, чем «простой» метод, но немного точнее. Точность по умолчанию не работает, поэтому приведенный ниже код значительно снижает эту точность. (dblquad выводит два числа: приблизительный интеграл и указание ошибки. Чтобы получить только интеграл, в коде используется dblquad()[0].) - Тот же подход можно использовать для других измерений. Для «простого» метода создайте сетку большего размера (
xx, yy, zz = np.mgrid[xmin:xmax:100j, ymin:ymax:100j, zmin:zmax:100j]). Обратите внимание, что деление на 1000 в каждом измерении приведет к созданию сетки, которая слишком велика для работы. - При использовании
scipy.integrate, dblquad необходимо заменить на tplquad для 3-х измерений или nquad для размеров N. Это, вероятно, также будет довольно медленным, поэтому необходимо дополнительно снизить точность.
import numpy as np
import pandas as pd
import scipy.stats as st
from scipy.integrate import dblquad
df1 = pd.DataFrame({'x':np.random.uniform(0, 1, 100), 'y':np.random.uniform(0, 1, 100)})
df2 = pd.DataFrame({'x':np.random.uniform(0, 1, 100), 'y':np.random.uniform(0, 1, 100)})
# Extract x and y
x1 = df1['x']
y1 = df1['y']
x2 = df2['x']
y2 = df2['y']
# Define the borders
deltaX = (np.max([x1, x2]) - np.min([x1, x2])) / 4
deltaY = (np.max([y1, y2]) - np.min([y1, y2])) / 4
xmin = np.min([x1, x2]) - deltaX
xmax = np.max([x1, x2]) + deltaX
ymin = np.min([y1, y2]) - deltaY
ymax = np.max([y1, y2]) + deltaY
# fit a gaussian kernel using scipy’s gaussian_kde method
kernel1 = st.gaussian_kde(np.vstack([x1, y1]))
kernel2 = st.gaussian_kde(np.vstack([x2, y2]))
print('volumes via scipy`s dblquad (volume):')
print(' volume_f1 =', dblquad(lambda y, x: kernel1((x, y)), xmin, xmax, ymin, ymax, epsabs=1e-4, epsrel=1e-4)[0])
print(' volume_f2 =', dblquad(lambda y, x: kernel2((x, y)), xmin, xmax, ymin, ymax, epsabs=1e-4, epsrel=1e-4)[0])
print(' volume_intersection =',
dblquad(lambda y, x: np.minimum(kernel1((x, y)), kernel2((x, y))), xmin, xmax, ymin, ymax, epsabs=1e-4, epsrel=1e-4)[0])
В качестве альтернативы можно вычислить среднее значение по сетке точек и умножить результат на площадь сетки. Обратите внимание, что np.mgrid намного быстрее, чем создание списка с помощью itertools.
# Create meshgrid
xx, yy = np.mgrid[xmin:xmax:1000j, ymin:ymax:1000j]
positions = np.vstack([xx.ravel(), yy.ravel()])
f1 = np.reshape(kernel1(positions).T, xx.shape)
f2 = np.reshape(kernel2(positions).T, xx.shape)
intersection = np.minimum(f1, f2)
print('volumes via the mean value multiplied by the area:')
print(' volume_f1 =', np.sum(f1) / f1.size * ((xmax - xmin)*(ymax - ymin)))
print(' volume_f2 =', np.sum(f2) / f2.size * ((xmax - xmin)*(ymax - ymin)))
print(' volume_intersection =', np.sum(intersection) / intersection.size * ((xmax - xmin)*(ymax - ymin)))
Пример вывода:
volumes via scipy`s dblquad (volume):
volume_f1 = 0.9946974276169385
volume_f2 = 0.9928998852123891
volume_intersection = 0.9046421634401607
volumes via the mean value multiplied by the area:
volume_f1 = 0.9927873844924111
volume_f2 = 0.9910132867915901
volume_intersection = 0.9028999384136771