[PYTHON] (Suite) Essayez d'autres fonctions de distance avec kmeans dans Scikit-learn

Que présenter dans cet article

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.

Try1 utilise la méthode cdist

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

Écrire avec Try2 Cython

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

Try3 Utiliser d'autres fonctions de distance

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.

41094b24-3bcb-0d04-93ab-edc54eb59501.png

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

J'ai essayé de mesurer

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.

Résumé

Pour utiliser d'autres fonctions de distance avec les k-means de scikit-learn,

Recommended Posts

(Suite) Essayez d'autres fonctions de distance avec kmeans dans Scikit-learn
kmeans ++ avec scikit-learn
Définissez votre propre fonction de distance avec k-means de scikit-learn
Essayez d'utiliser scikit-learn (1) - Clustering K par méthode moyenne
Identifiez les valeurs aberrantes avec le classificateur de forêt aléatoire de scikit-learn
SVM essayant l'apprentissage automatique avec scikit-learn
Essayez de vous connecter à qiita avec Python
Essayez SVM avec scikit-learn sur Jupyter Notebook
Essayez d'utiliser Python avec Google Cloud Functions
Regrouper les écoles représentatives à l'été 2016 avec scikit-learn
Remplissez les valeurs manquantes avec Scikit-learn impute
[Suite] Essayez l'accès au registre PLC avec Python
Essayez de travailler avec Mongo en Python sur Mac
Obtenez les distances de latitude et de longitude en mètres avec QGIS
Essayez de convertir des vidéos en temps réel avec OpenCV