Dans l'article précédent (http://qiita.com/Kensuke-Mitsuzawa/items/67c17f6bda8a1097740c), j'ai présenté que kmeans essaie des fonctions de distance autres que la distance euclidienne.
Vous pouvez également utiliser Euclidean. Cependant, il existe un cas où une autre fonction de distance est un meilleur regroupement. Alors pourquoi ne pas essayer? C'est le contenu.
Au fait, la dernière fois, je me suis retrouvé avec la conclusion: "J'ai essayé d'utiliser la fonction de distance du coefficient de corrélation de Pearson, mais bon!".
Ce n'est pas pratique. J'ai donc essayé diverses choses.
scipy a une méthode appelée cdist.
C'est une méthode très utile même si elle implémente une échelle de similarité sur une échelle de distance.
De plus, si vous n'avez pas d'échelle implémentée, vous pouvez la définir vous-même.
Définissons donc notre propre fonction de distance pour le coefficient de corrélation de Pearson.
#Fonction de distance utilisant cdist
def special_operation(X, Y):
coeficient = scipy.stats.pearsonr(X, Y)[0]
#coeficient = numpy.corrcoef(X, Y)[0][1]
if coeficient < 0:
return abs(coeficient)
else:
return 1 - coeficient
def special_pearsonr_corrcoef(X, Y):
return cdist(X, Y, special_operation)
special_operation ()
est une fonction qui prend les vecteurs X et Y et renvoie un float.
Si vous enveloppez cette fonction avec cdist ()
, vous aurez une méthode qui calcule complètement la distance entre les données en un rien de temps.
Remplacez la fonction euclidienne par le patch monkey comme auparavant.
# pearson distance with numpy cdist
def new_euclidean_distances(X, Y=None, Y_norm_squared=None, squared=False):
return special_pearsonr_corrcoef(X, Y)
from sklearn.cluster import k_means_
k_means_.euclidean_distances = new_euclidean_distances
kmeans_model = KMeans(n_clusters=3, random_state=10, init='random').fit(features)
print(kmeans_model.labels_)
elapsed_time = time.time() - start
print ("Special Coeficient k-means:{0}".format(elapsed_time))
print('='*40)
Les résultats seront présentés plus tard
Réimplémentation avec Cython pour des vitesses encore plus rapides.
J'ai pensé: "N'est-ce pas encore lent si j'appelle scipy même avec cython?", Alors je copie et colle la partie pertinente de scipy dans C.
Je fais aussi des définitions de type pour faire autant que je peux.
Compilez-le et, une fois terminé, modifiez-le à nouveau.
# -*- coding: utf-8 -*-
#cython: boundscheck=False
#cython: wraparound=False
__author__ = 'kensuke-mi'
import numpy as np
cimport numpy as np
cimport cython
from itertools import combinations
from scipy.stats import pearsonr
import time
import scipy.special as special
from scipy.special import betainc
ctypedef np.float64_t FLOAT_t
def _chk_asarray(np.ndarray[FLOAT_t, ndim=1] a, axis):
if axis is None:
a = np.ravel(a)
outaxis = 0
else:
a = np.asarray(a)
outaxis = axis
if a.ndim == 0:
a = np.atleast_1d(a)
return a, outaxis
def ss(np.ndarray[FLOAT_t, ndim=1] a, int axis=0):
"""
Squares each element of the input array, and returns the sum(s) of that.
Parameters
----------
a : array_like
Input array.
axis : int or None, optional
Axis along which to calculate. Default is 0. If None, compute over
the whole array `a`.
Returns
-------
ss : ndarray
The sum along the given axis for (a**2).
See also
--------
square_of_sums : The square(s) of the sum(s) (the opposite of `ss`).
Examples
--------
>>> from scipy import stats
>>> a = np.array([1., 2., 5.])
>>> stats.ss(a)
30.0
And calculating along an axis:
>>> b = np.array([[1., 2., 5.], [2., 5., 6.]])
>>> stats.ss(b, axis=1)
array([ 30., 65.])
"""
cdef np.ndarray[FLOAT_t, ndim=1] new_a
cdef int new_axis
new_a, new_axis = _chk_asarray(a, axis)
return np.sum(a*a, axis)
def betai(FLOAT_t a, FLOAT_t b, FLOAT_t x):
"""
Returns the incomplete beta function.
I_x(a,b) = 1/B(a,b)*(Integral(0,x) of t^(a-1)(1-t)^(b-1) dt)
where a,b>0 and B(a,b) = G(a)*G(b)/(G(a+b)) where G(a) is the gamma
function of a.
The standard broadcasting rules apply to a, b, and x.
Parameters
----------
a : array_like or float > 0
b : array_like or float > 0
x : array_like or float
x will be clipped to be no greater than 1.0 .
Returns
-------
betai : ndarray
Incomplete beta function.
"""
cdef FLOAT_t typed_x = np.asarray(x)
cdef FLOAT_t betai_x = np.where(typed_x < 1.0, typed_x, 1.0) # if x > 1 then return 1.0
return betainc(a, b, betai_x)
def pearsonr(np.ndarray[FLOAT_t, ndim=1] x, np.ndarray[FLOAT_t, ndim=1] y):
"""
Calculates a Pearson correlation coefficient and the p-value for testing
non-correlation.
The Pearson correlation coefficient measures the linear relationship
between two datasets. Strictly speaking, Pearson's correlation requires
that each dataset be normally distributed. Like other correlation
coefficients, this one varies between -1 and +1 with 0 implying no
correlation. Correlations of -1 or +1 imply an exact linear
relationship. Positive correlations imply that as x increases, so does
y. Negative correlations imply that as x increases, y decreases.
The p-value roughly indicates the probability of an uncorrelated system
producing datasets that have a Pearson correlation at least as extreme
as the one computed from these datasets. The p-values are not entirely
reliable but are probably reasonable for datasets larger than 500 or so.
Parameters
----------
x : (N,) array_like
Input
y : (N,) array_like
Input
Returns
-------
(Pearson's correlation coefficient,
2-tailed p-value)
References
----------
http://www.statsoft.com/textbook/glosp.html#Pearson%20Correlation
"""
# x and y should have same length.
cdef np.ndarray[FLOAT_t, ndim=1] typed_x = np.asarray(x)
cdef np.ndarray[FLOAT_t, ndim=1] typed_y = np.asarray(y)
cdef int n = len(typed_x)
cdef float mx = typed_x.mean()
cdef float my = typed_y.mean()
cdef np.ndarray[FLOAT_t, ndim=1] xm, ym
cdef FLOAT_t r_den, r_num, r, prob, t_squared
cdef int df
xm, ym = typed_x - mx, typed_y - my
r_num = np.add.reduce(xm * ym)
r_den = np.sqrt(ss(xm) * ss(ym))
r = r_num / r_den
# Presumably, if abs(r) > 1, then it is only some small artifact of floating
# point arithmetic.
r = max(min(r, 1.0), -1.0)
df = n - 2
if abs(r) == 1.0:
prob = 0.0
else:
t_squared = r**2 * (df / ((1.0 - r) * (1.0 + r)))
prob = betai(0.5*df, 0.5, df/(df+t_squared))
return r, prob
cdef special_pearsonr(np.ndarray[FLOAT_t, ndim=1] X, np.ndarray[FLOAT_t, ndim=1] Y):
cdef float pearson_coefficient
pearsonr_value = pearsonr(X, Y)
pearson_coefficient = pearsonr_value[0]
if pearson_coefficient < 0:
return abs(pearson_coefficient)
else:
return 1 - pearson_coefficient
def copy_inverse_index(row_col_data_tuple):
return (row_col_data_tuple[1], row_col_data_tuple[0], row_col_data_tuple[2])
def sub_pearsonr(np.ndarray X, int row_index_1, int row_index_2):
cdef float pearsonr_distance = special_pearsonr(X[row_index_1], X[row_index_2])
return (row_index_1, row_index_2, pearsonr_distance)
cdef XY_both_2dim(np.ndarray[FLOAT_t, ndim=2] X, np.ndarray[FLOAT_t, ndim=2] Y):
cdef int row = X.shape[0]
cdef int col = Y.shape[0]
cdef np.ndarray[FLOAT_t, ndim=2] pearsonr_divergence_set = np.zeros([row, col])
cdef float pearson_value = 0.0
for x_index in xrange(row):
x_array = X[x_index]
for y_index in xrange(col):
y_array = Y[y_index]
pearson_value = special_pearsonr(x_array, y_array)
pearsonr_divergence_set[x_index, y_index] = pearson_value
return pearsonr_divergence_set
def pearsonr_distances(X, Y, Y_norm_squared=None, squared=False):
#start = time.time()
cdef np.ndarray[FLOAT_t, ndim=2] pearsonr_divergence_set = XY_both_2dim(X, Y)
#elapsed_time = time.time() - start
#print ("Pearsonr(Cython version) k-means:{0}".format(elapsed_time))
return pearsonr_divergence_set
Comme je l'ai présenté la dernière fois, lorsque le clustering K-means est effectué à l'aide de diverses fonctions de distance, le résultat est une table.
En regardant cela, la similarité cos et la similitude Jackard (fonction de distance à partir de) donnent de bons résultats.
Il n'est pas nécessaire de s'en tenir au coefficient de corrélation de Pearson. Si vous obtenez de bons résultats, vous pouvez l'essayer.
Ces deux sont préparés après avoir été implémentés en tant que fonction de distance dans cdist
de scipy
, vous pouvez donc les utiliser immédiatement en écrivant simplement les paramètres.
distance cos
def new_euclidean_distances(X, Y=None, Y_norm_squared=None, squared=False):
return cdist(X, Y, 'cosine')
Distance de Jackard
def new_euclidean_distances(X, Y=None, Y_norm_squared=None, squared=False):
return cdist(X, Y, 'jaccard')
Je n'ai mesuré que le temps de traitement. Je n'ai pas vérifié la précision du clustering car je n'utilise pas de données étiquetées. (Parce que le but de cette fois est juste de me motiver à utiliser la méthode indiquée dans l'article)
Les données sont les mêmes que dans l'article précédent. Regardons-les dans l'ordre.
La première est la distance euclidienne normale. Les nombres correspondent à l'étiquette de clustering et à l'heure d'exécution.
[1 1 0 1 1 2 1 1 2 1 1 1 2 1 1 0 2 1 0 0 2 0 2 0 1 1 1 0 2 2 2 2 2 2 0 2 2]
Normal k-means:0.00853204727173
========================================
distance cos
[2 2 0 2 2 1 2 2 1 0 2 2 1 0 2 0 1 2 0 0 1 0 1 0 2 2 2 0 1 0 1 1 1 1 0 1 1]
Cosine k-means:0.0154190063477
========================================
Distance de Jackard
[0 0 0 0 1 0 0 0 2 0 0 0 1 0 0 0 0 0 2 2 0 2 0 2 0 0 0 0 2 2 0 0 0 0 0 0 0]
jaccard k-means:0.0656130313873
Coefficient de corrélation de Pearson utilisant cdist
[0 0 0 1 1 2 1 0 2 1 0 2 2 0 0 1 1 0 1 0 2 0 2 0 2 2 1 1 2 1 1 2 2 2 1 2 2]
Special Coeficient k-means:9.39525008202
========================================
coefficient de corrélation de Pearson cythonisé
[0 0 0 1 1 2 1 0 2 1 0 2 2 0 0 1 1 0 1 0 2 0 2 0 2 2 1 1 2 1 1 2 2 2 1 2 2]
Pearsonr k-means:9.40655493736
========================================
De toute évidence, le temps d'exécution du coefficient de corrélation de Pearson est inutile. La vitesse est
Euclidienne <Jackard <cos <(Mur insurmontable) <Pearson
était.
Le résultat de la distance du coefficient de Jackard est assez différent des autres étiquettes de regroupement ...
Peut-être que cette fois, les données d'entrée ont une petite quantité de fonctionnalités (une situation qui est très différente de celle du regroupement de documents), de sorte que le résultat peut être considérablement différent des autres.
À ce stade, je voudrais profiter d'une autre occasion pour le vérifier en profondeur.
Pour utiliser d'autres fonctions de distance avec les k-means de scikit-learn,
Recommended Posts