%matplotlib inline 
import random
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from matplotlib.patches import Ellipse
import scipy as sc
from scipy import random, linalg, stats, special

data=pd.read_csv('MoGDE.csv')

data

data.shape

colordict = {'abnormal': 'red', 'normal': 'blue'}
color = list(colordict[Class] for Class in data.label)

plt.scatter(data.iloc[:,0],data.iloc[:,1],c=color)

X = data.drop('label',axis=1)
# used only normal
x = X.iloc[list(y=='normal'),]
y = data['label']

N_var = x.shape[1]  # Variable의 갯수: 2
N_gaussian = 3      # 몇개의 가우시간 분포를 사용할 것인가?
N = x.shape[0]      # 데이터의 관측치 갯수

# EM step을 위한 초기 Mu, Cov 값 크기 설정
_mu = np.empty([N_gaussian, N_var])
_cov = np.empty([N_gaussian, N_var, N])


print("initial mu shape =" ,_mu.shape)
print("initial cov shape =", _cov.shape)
print("initial weight shape =", _w.shape)

# 각 gaussian 분포의 가중치, 어떤 분포를 더 중요깊게 봐야하나? 
_w = np.repeat(1.0/N_gaussian,N_gaussian)

# initial w
_w

# starting values
_mu = np.empty([N_gaussian, N_var])
for g in range(N_gaussian):
    _mu[g, :] = np.random.random(N_var) * np.max(x, axis=0)

# initial mu
_mu

def PSM_init(n):
    M = np.matrix(np.random.rand(n,n))
    M = (M + M.T)/2
    M = M + np.eye(n)
    return M # variance를 임의로 조정 

_cov = np.array([PSM_init(N_var) for i in range(0,N_gaussian)])

# 3개 가우시안과 2x2 covariance
_cov

#EM Algorithm: Expectation
def Estep(x, w, mu, cov):
    # 관측치 x 가우시안의 수 = (10,3)
    # wj*p(xi)를 저장
    w_x = np.zeros((x.shape[0], mu.shape[0]))
    # 관측치 갯수마다 적용
    for i in range(x.shape[0]):
        wj_xi = np.zeros(mu.shape[0]) # n gaussian

        for g in range(mu.shape[0]):
            # wj*p(xi)
            wj_xi[g] = w[g] * sc.stats.multivariate_normal.pdf(x.iloc[i, :], mu[g, :], cov[g, :, :])
        for g in range(wj_xi.shape[0]):
            # nomalized
            w_x[i, g] = wj_xi[g] / np.sum(wj_xi)

    return w_x

w_x = Estep(x, _w, _mu, _cov)
pd.DataFrame(w_x,columns=['gaussian01','gaussian02','gaussian03'])

#EM Algorithm: Maximazation
def Mstep(w_x, x, mu, cov):

    N = x.shape[0]

    # the weigths update
    new_w = np.sum(w_x, axis=0) / N
     
    # (N, N_gaussian, N_var)
    df_mu = np.zeros((N, mu.shape[0], mu.shape[1]))
    # (N, N_gaussian, nxn covariance)
    df_cov = np.zeros((N, cov.shape[0], cov.shape[1], cov.shape[2]))

    # 새롭게 업데이트 될 mu, Sigma
    new_mu = np.zeros((mu.shape[0], mu.shape[1]))
    new_cov = np.zeros((cov.shape[0], cov.shape[1], cov.shape[2]))

# mean
    for i in range(N):
        # (N_gaussian x 1) x (1 x N_var) = (N_gaussian, N_var)
        df_mu[i, :, :] = np.outer(w_x[i, :], x.iloc[i, :])
    
    for g in range(N_gaussian):
        new_mu[g, :] = (1 / np.sum(w_x, axis=0)[g]) * np.sum(df_mu, axis=0)[g, :]

# sd
    for i in range(N):
        for g in range(N_gaussian):
            df_cov[i, g, :, :] = w_x[i, g] * np.outer((x.iloc[i, :] - new_mu[g, :]),
                                                                (x.iloc[i, :] - new_mu[g, :]))
            
    for g in range(N_gaussian):
        new_cov[g, :, :] = (1 / np.sum(w_x, axis=0)[g]) * np.sum(df_cov, axis=0)[g, :, :]

    return new_w, new_mu, new_cov

new_w, new_mu, new_cov = Mstep(w_x, x, _mu, _cov)
print('new_w=',new_w)
print('new_mu=',new_mu)
print('new_cov=',new_cov)

# initialization
iteration = 100
N_gaussian = 2

_w = np.repeat(1.0/N_gaussian,N_gaussian)

_cov = np.array([PSM_init(N_var)/3 for i in range(0,N_gaussian)])

_mu = np.empty([N_gaussian, N_var])

for g in range(N_gaussian):
    _mu[g, :] = np.random.random(N_var) * np.max(x, axis=0)

for k in range(iteration):

# E step
    w_x = Estep(x, _w, _mu, _cov) # 주어진 모수로 w*p(x)를 계산한

# M step
    new_w, new_mu, new_cov = Mstep(w_x, x, _mu, _cov)

# log likelihood를 계산 
    logL_i = np.zeros((x.shape[0]))

    for i in range(x.shape[0]):

        L = np.zeros(N_gaussian)

        for g in range(N_gaussian):
            L[g] = new_w[g] * sc.stats.multivariate_normal.pdf(x.iloc[i, :], new_mu[g, :],
                                                                        new_cov[g, :, :])

        logL_i[i] = np.log(np.sum(L))
        
    if(k%20==0):print(k,'th logL=',logL)
    
    logL = np.sum(logL_i)
    

print ("mu=", new_mu)
print("sigma=",new_cov)
print("weight=",new_w)

# visualization

plotsize = 8
sizeMean = 30
text_size = 16
axis_font = {'fontname': 'Arial', 'size': '24'}
Title_font = {'fontname': 'Arial', 'size': '28'}
color = ['b', 'r', 'g', 'c', 'm', 'y', 'k']

fig = plt.figure()

ax = fig.add_subplot(1, 1, 1)
ax.plot(X.iloc[np.array(data.label=='normal'), 0], X.iloc[np.array(data.label=='normal'), 1], 'k.', markersize=sizeMean, color=color[0])
ax.plot(X.iloc[np.array(data.label=='abnormal'), 0], X.iloc[np.array(data.label=='abnormal'), 1], 'k.', markersize=sizeMean,color=color[1])


for i in range(N_gaussian):
    # the sd with ellipses
    # central point of the error ellipse
    pos = [new_mu[i, 0], new_mu[i, 1]]

    # for the angle we need the eigenvectors of the covariance matrix
    w, ve = np.linalg.eig(new_cov[i, 0:2, 0:2])

    # We pick the largest eigen value
    order = w.argsort()[::-1]
    w = w[order]
    ve = ve[:, order]

    # we compute the angle towards the eigen vector with the largest eigen value
    thetaO = np.degrees(np.arctan(ve[1, 0] / ve[0, 0]))

    # Compute the width and height of the ellipse based on the eigen values (ie the length of the vectors)
    width, height = 2 * np.sqrt(w)

    # making the ellipse
    ellip = Ellipse(xy=pos, width=width, height=height, angle=thetaO)
    ellip.set_alpha(0.5)
    ellip.set_facecolor(color[i])

    ax.add_artist(ellip)

for label in (ax.get_xticklabels() + ax.get_yticklabels()):
    label.set_fontname('Arial')
    label.set_fontsize(text_size)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.xaxis.set_ticks_position('bottom')
ax.yaxis.set_ticks_position('left')
ax.set_ylabel('X2', **axis_font)
ax.set_xlabel('X1', **axis_font)
ax.set_title('Novelty example', y=1.08, **Title_font)
ax.figure.set_size_inches(plotsize, plotsize)

plt.show()

score = list()
for i in range(10):
    novelty = [new_w[g] * sc.stats.multivariate_normal.pdf(
        X.iloc[i, :], new_mu[g, :], new_cov[g, :, :]) for g in range(N_gaussian)]
    score.append(1/np.sum(novelty))

pd.concat([X,data['label'],pd.DataFrame(np.log(score), columns=['Novelty score'])],axis=1)