Эти единичные векторы получаются в результате диагонализации тензора инерции, поэтому направление определено, но смысл вырожден (это будет важно позже).
Для этого я использую функцию плотности ядра scikit-learn (sklearn.neighbors.KernelDensity).
Для воспроизводимости я приведу пример векторов и весов (но последние не совсем важны для вопроса):
Код: Выделить всё
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import GridSearchCV
from sklearn.model_selection import LeaveOneOut
from sklearn.neighbors import KernelDensity
nn = np.array([[ 0.03402818, 0.39301999, 0.91890009],
[ 0.59206726, 0.80298496, 0.06834847],
[-0.92554092, 0.34986418, -0.144807 ],
[ 0.60773445, 0.79307238, -0.04117088],
[ 0.18615284, 0.51452414, 0.8370257 ],
[ 0.40394095, 0.22357499, 0.88704336],
[ 0.16339813, 0.48618174, 0.85844532],
[ 0.86968201, -0.24154248, 0.43047697],
[-0.50666593, -0.83360426, -0.21998538],
[-0.80594915, 0.46939866, -0.36070884],
[ 0.03807831, 0.40621088, 0.91298563],
[ 0.2755237 , 0.22552711, 0.93446467],
[ 0.02600335, 0.38454888, 0.92273831],
[-0.96124128, 0.18894847, -0.20078269],
[ 0.60268438, 0.79609915, 0.05475108],
[-0.94328141, 0.30247806, -0.13684734],
[ 0.43765748, 0.22873265, 0.86956156],
[ 0.17128545, 0.50535992, 0.84573793],
[ 0.39959667, 0.22321366, 0.88909964],
[ 0.88907925, -0.21332485, 0.40500691],
[-0.14413127, 0.18262637, 0.97256043],
[ 0.27327328, 0.0623681 , 0.95991246],
[-0.13116824, 0.18737418, 0.97349156],
[ 0.94248182, -0.22452222, 0.24762429],
[-0.43556423, 0.19004559, 0.87986731],
[ 0.94787827, -0.25737596, 0.1878414 ],
[-0.46914909, -0.00164206, -0.88311745],
[-0.05133885, 0.18431249, 0.98152597],
[ 0.43146423, 0.015442 , 0.90199787],
[-0.88801615, 0.2035849 , -0.41228692],
[ 0.29721996, 0.0539346 , 0.95328451],
[-0.14843932, 0.18290349, 0.97186011],
[ 0.26234948, 0.06603363, 0.96271092],
[ 0.94420917, -0.22516556, 0.2403529 ],
[-0.45932668, -0.00433692, -0.88825683]])
weights = np.array([2.75792219, 2.80864853, 1.59750041, 4.20455058, 3.24719617,
1.70647554, 3.74628944, 5.88323566, 2.60350748, 3.29570828,
3.60952358, 1.21608154, 3.13088853, 1.48962109, 3.21465176,
1.77468426, 1.50534637, 4.14162718, 1.54502792, 3.18413058,
3.64251358, 3.65384557, 8.13181441, 3.02696656, 4.58928026,
5.16233972, 3.38317826, 4.24850591, 4.63408708, 3.04477234,
4.77449734, 3.59965002, 4.01439879, 3.25976281, 3.66809761])
Код: Выделить всё
def _vec_to_coord(vec):
"""
Form a set of unit vectors,
to longitude and latitude
"""
longitude = np.arctan2(
vec[:, 1],
vec[:, 0],
)
latitude = np.arcsin(vec[:, 2])
return longitude, latitude
def _gen_uGrid():
"""
Generate a longitude-latitude grid, and the corresponding unit vectors.
"""
lon, lat = np.mgrid[-np.pi : np.pi : 100j, -np.pi / 2 : np.pi / 2 : 100j]
nn = np.array(
[
(np.cos(lon) * np.cos(lat)).ravel(),
(np.sin(lon) * np.cos(lat)).ravel(),
np.sin(lat).ravel(),
]
).T
return lon, lat, nn
- Проецируем свои векторы на unit-sphere и использовать метрику хаверсинуса.
Чтобы принять во внимание вырождение, я продублировал набор данных на другую половину сферы. - Используйте специальную метрику (которая на 1 меньше абсолютного значения скалярного произведения):
Код: Выделить всё
def custom_metric(x, y):
"""
Compute the distance between two unit vectors
"""
dist = 0
for ii in range(len(x)):
dist += x[ii] * y[ii]
return 1 - abs(dist)
Для первого случая (т.е. ядра с метрикой хаверсинуса):
Код: Выделить всё
longitude, latitude = _vec_to_coord(nn)
# wrap ol the vectors in one half of the sphere and duplicate them
msk = ~((longitude > -(np.pi / 2)) & (longitude < (np.pi / 2)))
longitude[msk] = longitude[msk] - np.sign(longitude[msk]) * np.pi
latitude[msk] = -latitude[msk]
# duplicate the sample to properly take into account the periodicity
duplicated_longitude = longitude - np.sign(longitude) * np.pi
duplicated_latitude = -latitude
latitude = np.concatenate([latitude, duplicated_latitude])
longitude = np.concatenate([longitude, duplicated_longitude])
sample_weights = np.concatenate([weights, weights])
# the haversine metric needs the latitude first
projected_coords = np.stack([latitude, longitude]).T
# estiamte the best band with via LeaveOneOut cross-validation
nbins = 50
bandwidths = 10 ** np.linspace(-2, 1, nbins)
grid = GridSearchCV(
KernelDensity(kernel="gaussian", metric="haversine"),
{"bandwidth": bandwidths},
cv=LeaveOneOut(),
n_jobs=6,
verbose=0,
)
grid.fit(projected_coords, sample_weight=sample_weights)
# Now the ploting
projection = "mollweide"
fig, ax = plt.subplots(subplot_kw=dict(projection=projection))
X, Y = np.mgrid[-np.pi : np.pi : 100j, -np.pi / 2 : np.pi / 2 : 100j]
xy = np.vstack([Y.ravel(), X.ravel()]).T # the latitud first
Z = np.exp(grid.best_estimator_.score_samples(xy)).reshape(X.shape)
img = ax.pcolormesh(X, Y, Z)
fig.colorbar(img, orientation="horizontal")
ax.plot(longitude, latitude, ".", color="r")
ax.plot(duplicated_longitude, duplicated_latitude, ".", color="b")
ax.grid()
ax.set_title("haversine")
Код: Выделить всё
sample_weights = weights
# estiamte the best band with via LeaveOneOut cross-validation
bandwidths = 10 ** np.linspace(-3, 0, nbins)
grid = GridSearchCV(
KernelDensity(
kernel="gaussian",
metric="pyfunc",
metric_params={"func": custom_distance},
),
{"bandwidth": bandwidths},
cv=LeaveOneOut(),
n_jobs=6,
verbose=0,
)
grid.fit(nn, sample_weight=weights)
X, Y, samples = _gen_uGrid()
Z = np.exp(grid.best_estimator_.score_samples(samples)).reshape(X.shape)
fig, ax = plt.subplots(subplot_kw={"projection": "mollweide"})
img = ax.pcolormesh(X, Y, Z)
fig.colorbar(img, orientation="horizontal")
ax.plot(*_vec_to_coord(nn), "r.")
ax.grid()
ax.set_title("Custom Metric")
В чем может быть причина такого масштабирования плотностей?
haversine
Пользовательский
Подробнее здесь: https://stackoverflow.com/questions/782 ... -learn-kde
