# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# %matplotlib inline
import pdb


# ## 【問題1】中心点の初期値を決める
# 入力されたデータから $ K $ 個の中心点 $ μ1 $ から $ μ_K $ の初期値を決めるコードを作成してください。  
# $ K $ は指定したクラスタ数です。  
# 最もシンプルな初期値の決定方法は、データ点 $ X_n $ の中からランダムに $ K $ 個選ぶことです。今回はこれを実装してください。  
# K-meansの学習結果は中心点 $ μ $ の初期値に影響を受けます。  
# そのため、学習時には複数個の初期値で計算を行います。

# ## 【問題2】SSEを求める関数の作成
# クラスタ内誤差平方和（SSE, Sum of Squared Errors）を計算する関数を作成してください。  
# K-meansはこのSSEを最小化する $ r_{nk} $ と $ μ_k $ を求めることが目的となります。  
# 複数個の初期値で計算したクラスタリング結果から、どれを最終的に採用するかを決める際にこのSSEを求める関数を使用します。
# $$
# SSE = \sum^{N}_{n=1} \sum^{K}{k=1}r_{nk} |X_n - \mu_k|^2
# $$

# ## 【問題3】クラスタへの割り当て
# 全てのデータ点 $X_n $ を最も近い中心点 $ μ_k $ に割り当てるコードを作成してください。  
# K-menasにおける **近い** とは点と点のユークリッド距離が小さくなることです。ユークリッド距離とはピタゴラスの定理（三平方の定理）で求められるものですが、ベクトル $ p, q $ に対しては以下の数式で表現できます。
# $$
# |q - p| = \sqrt{(q-p)\cdot(q-p)}
# $$

# ## 【問題4】中心点の移動
# 中心点 $ μ_k $ を $ k $ 番目のクラスタに割り当てられる全てのデータ点 $ X_n $ の平均値（重心）に移動するコードを作成してください。  
# クラスタの割り当て $ r_{nk} $ を固定した上でSSEを最小化していることになりますが、SSE自体を求める必要はありません。

# ## 【問題5】繰り返し
# クラスタへの割り当てと中心点の移動を繰り返すことで学習を行うコードを作成してください。  
# 次のどれかの条件を満たしたらこれを止めます。
# - 繰り返しが指定した回数に到達する。（雛形では引数max_iter）
# - 中心点 $ μ_k $ と重心の差が指定した許容値以下になる。（雛形では引数tol）

# ## 【問題6】異なる初期値で計算
# 問題5の繰り返しを、指定した回数初期値を変更して計算するコードを作成してください。回数は雛形では引数n_initで指定します。  
# 複数計算した中でSSEが最小になったものが最終的な学習結果となるため、その時の中心点 $ μ_k $ を保持しておきます。

# ## 【問題7】推定
# 新たなデータに対してどのクラスタに所属するかの推定の際は、学習時の1つ目のフェーズと同じことを行います。  
# - データ点 $ X_n $ と学習で決まった中心点 $ μ_k $ に対し、SSEを最小化する $ r_{nk} $  を選ぶ。  
#
# 具体的には次のことです。
#
# - データ点 $ X_n $ を最も近い $ μ_k $ に割り当てる。

# +
class ScratchKMeans():

    def __init__(self, n_clusters=2, n_init=3, max_iter=30, tol=1e-5):
        # ハイパーパラメータを属性として記録
        self.n_clusters = n_clusters
        self.n_init = n_init
        self.max_iter = max_iter
        self.tol = tol
        
        self.X = None
        self.m = None
        self.mu, self.best_mu = None, None
        self.clusters, self.best_clusters = None, None
        self.sse, self.best_sse = None, None
        
    def fit(self, X):
        self.X = X
        self.m = len(self.X)
        #self.clusters = np.random.choice(self.n_clusters, self.m)
        self.best_sse = np.inf
        
        for n in range(self.n_init):
            self.sse_list = [np.inf]
            self.mu = X[np.random.choice(self.m, size=self.n_clusters, replace=True)].copy()
            count = 0
            while count < self.max_iter:
                self.clusters = self.asign_cluster(self.X)
                for i in range(self.n_clusters):
                    self.mu[i] = np.mean(self.X[self.clusters==i], axis=0)
                if self.sse_list[-1] - self.calc_sse() < self.tol:
                    break
                self.sse_list.append(self.calc_sse())
                count += 1
            
            if self.calc_sse() < self.best_sse:
                self.best_sse = self.calc_sse()
                self.best_mu = self.mu
                self.best_clusters = self.clusters

    
    def asign_cluster(self, X, predict=False):
        m = len(X)
        tmp = np.empty((m, self.n_clusters))
        if predict:
            mu = self.best_mu
        else:
            mu = self.mu
        for i in range(self.n_clusters):
            for j in range(m):
                tmp[j,i] = np.linalg.norm(X[j]-mu[i])
        return np.argmin(tmp, axis=1)   
    
    
    def calc_sse(self):
        sse = 0
        for i in range(self.n_clusters):
            sse += np.sum((self.X[self.clusters==i] - self.mu[i])**2)
        return sse
        
        
    def predict(self, X_test):
        y_pred = self.asign_cluster(X_test, predict=True)
        return y_pred
    

# -

from sklearn.datasets import make_blobs
X, _ = make_blobs(n_samples=100, n_features=2, centers=4, cluster_std=0.5, shuffle=True, random_state=0)

kmeans = ScratchKMeans(n_clusters=4)
kmeans.fit(X)

kmeans.sse_list

plt.scatter(X[:,0], X[:,1])
plt.scatter(kmeans.best_mu[:,0], kmeans.best_mu[:,1], marker='o')

# ## 【問題8】エルボー法の実装
# K-meansにおいて中心点の個数 $ k $ は人為的に決める必要がありますが、その判断方法としてエルボー法があります。  
# エルボー法により中心点の個数 $ k $ を検討するためのコードを作成してください。  
# SSEを縦軸、クラスタ数 $ k $ を横軸にした折れ線グラフを表示します。  
# このためにK-meansのインスタンスからはSSEを取得できるようにしてください。

n_elbow = 7
score = np.zeros(n_elbow)
for i in range(n_elbow):
    kmeans = ScratchKMeans(n_clusters=i+1)
    kmeans.fit(X)
    score[i] = kmeans.best_sse
plt.plot(np.arange(n_elbow)+1, score, '-o')
plt.xlabel('n_clusters')
plt.ylabel('SSE')
plt.show();


# ## 【問題9】シルエット図
# $ k $ を決める助けになるものとしてシルエット図もあります。  
# シルエット図を表示するコードを作成してください。  
# Matplotlibを使い図を表示するコードは用意しました。  
# n番目のデータ点に対するシルエット係数 $ s_n $ は以下の数式でを求められます。  
# この数式を実装し、シルエット図を確認してください。
# $$
# s_n = \frac{b_n-a_n}{max(b_n,a_n)}
# $$

def calc_silhouette(X, mu, clusters):
    m = len(X)
    silhouette_vals = np.zeros(m)
    for i in range(m):
        a = np.mean(np.sum((X[clusters==clusters[i]]-X[i])**2, axis=1))
        mu_dist = np.sum((mu-mu[clusters[i]])**2, axis=1)
        idx = np.argsort(mu_dist)[1]
        b = np.mean(np.sum((X[clusters==idx]-X[i])**2, axis=1))
        silhouette_vals[i] = (b - a) / max(b, a)
    return silhouette_vals


from matplotlib import cm
def plot_silhouette(silhouette_vals, silhouette_avg, y_km, cluster_labels, n_clusters):
    y_ax_lower, y_ax_upper = 0, 0
    yticks = []
    for i, c in enumerate(cluster_labels):
        c_silhouette_vals = silhouette_vals[y_km == c]
        c_silhouette_vals.sort()
        y_ax_upper += len(c_silhouette_vals)
        color = cm.jet(i / n_clusters)
        plt.barh(range(y_ax_lower, y_ax_upper), c_silhouette_vals, height=1.0, edgecolor='none', color=color)
        yticks.append((y_ax_lower + y_ax_upper) / 2)
        y_ax_lower += len(c_silhouette_vals)
    plt.axvline(silhouette_avg, color="red", linestyle="--")
    plt.yticks(yticks, cluster_labels + 1)
    plt.ylabel('Cluster')
    plt.xlabel('Silhouette coefficient')
    plt.show()
#plot_silhouette(silhouette_vals, silhouette_avg, y_km, cluster_labels, n_clusters)


# +
kmeans = ScratchKMeans(n_clusters=4)
kmeans.fit(X)

silhouette_vals = calc_silhouette(X, kmeans.best_mu, kmeans.best_clusters)
silhouette_avg = np.mean(silhouette_vals)
y_km = kmeans.best_clusters
cluster_labels = np.arange(kmeans.n_clusters)
n_clusters = kmeans.n_clusters
plot_silhouette(silhouette_vals, silhouette_avg, y_km, cluster_labels, n_clusters)
# -

# ## 【問題10】クラスタ数kの選択
# 最適だと思われるクラスタ数 k を選択してください。エルボー法やシルエット図を活用する他、人間が後々どう利用したいのかも選択に関わってきます。  
# 選択した上で、なぜそれを選んだか説明してください。

data = pd.read_csv('Wholesale customers data.csv')
data.head()

data.info()

from sklearn.decomposition import PCA
NUM = 6
pca = PCA(n_components=NUM)
pca = pca.fit(data)
var_exp = pca.explained_variance_ratio_
cum_var_exp = np.cumsum(var_exp)
print(cum_var_exp)
plt.bar(range(1,NUM+1), var_exp, alpha=0.5, align='center', label='individual explained variance')
plt.step(range(1,NUM+1), cum_var_exp, where='mid', label='cumulative expalined variance')
plt.ylabel('Explained variance ratio')
plt.xlabel('Principal components')
plt.hlines(0.7, 0, 6,  "blue", linestyles='dashed')
plt.legend(loc='best')
plt.grid()
plt.show()

# PCAの結果では次元4で元データの分散の95%を維持できている。

data_dropped = data.drop(columns=['Region', 'Channel']).values

# ### RegionとChannel以外のデータでクラスタリング

n = 10
score = np.zeros(n)
for i in range(n):
    kmeans = ScratchKMeans(n_clusters=i+1)
    kmeans.fit(data_dropped)
    score[i] = kmeans.best_sse
plt.plot(np.arange(n)+1, score, '-o')
plt.xlabel('n_clusters')
plt.ylabel('SSE')
plt.show();

# エルボー図ではk=5, 6, 7, 8がよさそう。

for i in [5, 6, 7, 8]:
    kmeans = ScratchKMeans(n_clusters=i)
    kmeans.fit(data_dropped)

    silhouette_vals = calc_silhouette(data_dropped, kmeans.best_mu, kmeans.best_clusters)
    silhouette_avg = np.mean(silhouette_vals)
    y_km = kmeans.best_clusters
    cluster_labels = np.arange(kmeans.n_clusters)
    n_clusters = kmeans.n_clusters
    plot_silhouette(silhouette_vals, silhouette_avg, y_km, cluster_labels, n_clusters)

# シルエット図ではおおきな違いはなさそう。

# ### PCA処理後のデータでクラスタリング

pca = PCA(n_components = 2)
pca = pca.fit(data)
data_pca = pca.transform(data)

n = 10
score = np.zeros(n)
for i in range(n):
    kmeans = ScratchKMeans(n_clusters=i+1)
    kmeans.fit(data_pca)
    score[i] = kmeans.best_sse
plt.plot(np.arange(n)+1, score, '-o')
plt.xlabel('n_clusters')
plt.ylabel('SSE')
plt.show();

# エルボー図ではk=5, 6, 7, 8がよさそう。

for i in [5, 6, 7, 8]:
    kmeans = ScratchKMeans(n_clusters=i)
    kmeans.fit(data_pca)

    silhouette_vals = calc_silhouette(data_pca, kmeans.best_mu, kmeans.best_clusters)
    silhouette_avg = np.mean(silhouette_vals)
    y_km = kmeans.best_clusters
    cluster_labels = np.arange(kmeans.n_clusters)
    n_clusters = kmeans.n_clusters
    plot_silhouette(silhouette_vals, silhouette_avg, y_km, cluster_labels, n_clusters)

# シルエット図ではおおきな違いはなさそう。

#Channelの値
data['Channel'].value_counts()

#Regionの値
data['Region'].value_counts()

# クラスタ数は6とする。  
# 理由は、Channelは2値、Regionは3値からなる特徴量であり、これらの特徴量を組み合わせて6値に分類できると考えられるため。  
# またデータは、Fresh, Milk, Grocery, Frozen, Detergents_Paper, Delicassenの6つの特徴量から成るのでそれぞれのクラスタがそれぞれの列に対応づけられることが期待できるため。

# ## 【問題11】既知のグループとの比較
# RegionとChannelはK-meansに使用しなかったが、これらは既知のグループです。  
# クラスタとこれらのグループがどう違うか確認し、説明してください。

kmeans = ScratchKMeans(n_clusters=6)
kmeans.fit(data_pca)
y_pred = kmeans.predict(data_pca)
# plt.scatter(data_pca[data['Region']==1, 0], data_pca[data['Region']==1, 1], label='Region=1')
# plt.scatter(data_pca[data['Region']==2, 0], data_pca[data['Region']==2, 1], label='Region=2')
# plt.scatter(data_pca[data['Region']==3, 0], data_pca[data['Region']==3, 1], label='Region=3')
plt.scatter(data_pca[kmeans.clusters==0, 0], data_pca[kmeans.clusters==0, 1], label='Cluster0', s=80, alpha=0.4)
plt.scatter(data_pca[kmeans.clusters==1, 0], data_pca[kmeans.clusters==1, 1], label='Cluster1', s=80, alpha=0.4)
plt.scatter(data_pca[kmeans.clusters==2, 0], data_pca[kmeans.clusters==2, 1], label='Cluster2', s=80, alpha=0.4)
plt.scatter(data_pca[kmeans.clusters==3, 0], data_pca[kmeans.clusters==3, 1], label='Cluster3', s=80, alpha=0.4)
plt.scatter(data_pca[kmeans.clusters==4, 0], data_pca[kmeans.clusters==4, 1], label='Cluster4', s=80, alpha=0.4)
plt.scatter(data_pca[kmeans.clusters==5, 0], data_pca[kmeans.clusters==5, 1], label='Cluster5', s=80, alpha=0.4)
#plt.scatter(kmeans.best_mu[:,0], kmeans.best_mu[:,1], marker='o', label='cluster', c='r', s=100)
plt.legend();

plt.scatter(data_pca[data['Channel']==1, 0], data_pca[data['Channel']==1, 1], label='Channel=1')
plt.scatter(data_pca[data['Channel']==2, 0], data_pca[data['Channel']==2, 1], label='Channel=2')
plt.scatter(kmeans.best_mu[:,0], kmeans.best_mu[:,1], marker='o', label='cluster', c='r', s=100)
plt.legend();

plt.scatter(data_pca[data['Region']==1, 0], data_pca[data['Region']==1, 1], label='Region=1')
plt.scatter(data_pca[data['Region']==2, 0], data_pca[data['Region']==2, 1], label='Region=2')
plt.scatter(data_pca[data['Region']==3, 0], data_pca[data['Region']==3, 1], label='Region=3')
plt.scatter(kmeans.best_mu[:,0], kmeans.best_mu[:,1], marker='o', label='cluster', c='r', s=100)
plt.legend();

# Channel==1はcluster0,2,5に、Channel==2はcluster1,3,4に属する傾向がみえた。  
# Regionについては傾向はみられなかった。

c1r1 = (data['Channel']==1) & (data['Region']==1)
c1r2 = (data['Channel']==1) & (data['Region']==2)
c1r3 = (data['Channel']==1) & (data['Region']==3)
c2r1 = (data['Channel']==2) & (data['Region']==1)
c2r2 = (data['Channel']==2) & (data['Region']==2)
c2r3 = (data['Channel']==2) & (data['Region']==3)

plt.scatter(data_pca[c1r1, 0], data_pca[c1r1, 1], label='c1r1')
plt.scatter(data_pca[c1r2, 0], data_pca[c1r2, 1], label='c1r2')
plt.scatter(data_pca[c1r3, 0], data_pca[c1r3, 1], label='c1r3')
plt.scatter(data_pca[c2r1, 0], data_pca[c2r1, 1], label='c2r1')
plt.scatter(data_pca[c2r2, 0], data_pca[c2r2, 1], label='c2r2')
plt.scatter(data_pca[c2r3, 0], data_pca[c2r3, 1], label='c2r3')
plt.legend();

# ## 【問題12】卸売業者にとって有益な情報
# 卸売業者にとって有益な情報を探し出し文章にまとめてください。

data[['Fresh','Milk','Grocery','Frozen','Detergents_Paper', 'Delicassen']].describe().T

data.loc[c1r1,['Fresh','Milk','Grocery','Frozen','Detergents_Paper', 'Delicassen']].describe().T

# Fresh, Frozen以外はあまり買わないグループ。

data.loc[c1r2,['Fresh','Milk','Grocery','Frozen','Detergents_Paper', 'Delicassen']].describe().T

# Frozenをよく買う、Freshは平均程度、Milk, Grocery, Detergents_Paper, Dilicassenはかわないグループ。

data.loc[c1r3,['Fresh','Milk','Grocery','Frozen','Detergents_Paper', 'Delicassen']].describe().T

# どの品も平均程度に買うグループ。一番多い。

data.loc[c2r1,['Fresh','Milk','Grocery','Frozen','Detergents_Paper', 'Delicassen']].describe().T

# Milk, Grocery, Detergents_Paperをよく買うグループ。

data.loc[c2r2,['Fresh','Milk','Grocery','Frozen','Detergents_Paper', 'Delicassen']].describe().T

# Grocery, Detergents_Paperをよく買うグループ。

data.loc[c2r3,['Fresh','Milk','Grocery','Frozen','Detergents_Paper', 'Delicassen']].describe().T

# Milk, Grocery, Detergents_Paperをよく買うグループ。  
# 全体の25％程度で2番目に多い。

# #### DBSCAN

from sklearn import cluster, mixture, datasets

dataset1 = datasets.make_blobs(n_samples=1000, random_state=10, centers=6, cluster_std=1.2)[0]
dataset2 = datasets.make_moons(n_samples=1000, noise=0.05)[0]


def cluster_plots(set1, set2, colours1 = 'gray', colours2 = 'gray', title1 = 'Dataset 1', title2 = 'Dataset 2'):
 
    fig,(ax1,ax2) = plt.subplots(1, 2)
    fig.set_size_inches(6, 3)
 
    ax1.set_title(title1,fontsize=14)
    ax1.set_xlim(min(set1[:,0]), max(set1[:,0]))
    ax1.set_ylim(min(set1[:,1]), max(set1[:,1]))
    ax1.scatter(set1[:, 0], set1[:, 1],s=8,lw=0,c= colours1)
 
    ax2.set_title(title2,fontsize=14)
    ax2.set_xlim(min(set2[:,0]), max(set2[:,0]))
    ax2.set_ylim(min(set2[:,1]), max(set2[:,1]))
    ax2.scatter(set2[:, 0], set2[:, 1],s=8,lw=0,c=colours2)
 
    fig.tight_layout()
    plt.show()


cluster_plots(dataset1,  dataset2)

dbscan_dataset1 = cluster.DBSCAN(eps=1, min_samples=5, metric='euclidean').fit_predict(dataset1)
dbscan_dataset2 = cluster.DBSCAN(eps=1, min_samples=5, metric='euclidean').fit_predict(dataset2)

cluster_plots(dataset1, dataset2, dbscan_dataset1, dbscan_dataset2)

dbscan_dataset1 = cluster.DBSCAN(eps=0.1, min_samples=5, metric='euclidean').fit_predict(dataset1)
dbscan_dataset2 = cluster.DBSCAN(eps=0.1, min_samples=5, metric='euclidean').fit_predict(dataset2)

cluster_plots(dataset1, dataset2, dbscan_dataset1, dbscan_dataset2)

np.sum(dbscan_dataset1==-1)

StarcoderdataJupyterScriptsDedupFiltered