# Untitled

## 1865 days ago by Buyan

## functions and script for Clustering with Dirichlet Process(sort of probability process) ## written by DongHyuk.yi import numpy as np from numpy.random import multinomial as mnrand import matplotlib.pyplot as plt import matplotlib.cm as cm import warnings warnings.filterwarnings('ignore') # for warning off np.seterr(all='ignore') # for warning off ## function for 'PDF of Multivariate-normal distribution' # In sage, couldn't use 'mvnorm module in Scipy.stats', so define alternative function def mvnpdf(x,mu,cov): diff = mu - x diff = np.reshape(diff,(1,2)) diff_trans = np.matrix.transpose(diff) diff_trans = matrix(diff_trans) # for SAGE diff = matrix(diff) cov_inv = np.linalg.inv(cov) cov_inv = matrix(cov_inv) # for SAGE #multiply1 = np.matmul(diff, cov_inv) multiply1 = diff*cov_inv #high_part = np.exp(-0.5*np.matmul(multiply1, diff_trans)) high_part = -0.5 * multiply1 * diff_trans # for SAGE phi = np.pi #low_part = np.sqrt(det(2*phi*cov)) low_part = np.sqrt(det(matrix(2*phi*cov))) # For SAGE pdf = high_part/low_part return pdf ## function for conversion "sigma list" to "covariance matrix" def calc_cov(sigma): cov = np.diag((sigma)) return cov ## Main function def DP(data, alpha, sigma, init_cluster, iteration=10): ## data = raw data, alpha = concentration paramter, sigma = likely, variation of DP, init_cluster = initial number of cluster, iteration = number of whole process(default : 10000) print('function message : start') # set value with input data data = np.array(data); nr_data = data.shape dp_alpha = alpha; dp_cov = calc_cov(sigma); nr_cluster = init_cluster ## GMM with DP # 1. set arbitrary cluster for each data dp_cluster = np.random.randint(low=1, high=nr_cluster + 1, size=(nr_data)) count = 0 # calc initial centers of each cluster for ci in np.unique(dp_cluster): cluster_index = np.nonzero(dp_cluster == ci) center_pos = np.mean(data[cluster_index], axis=0) if count == 0: dp_center = center_pos else: dp_center = np.vstack([dp_center, center_pos]) count += 1 p = np.zeros(1) lik_list = np.zeros(1) for iter in range(iteration): print('%d iteration is doing' % iter) for i in range(0, nr_data): cur_point = data[i] p = np.zeros(nr_cluster + 1) lik_list = np.zeros(nr_cluster) for j in range(0, nr_cluster): ## Do Polya urn sequance(process) i_cluster = dp_cluster[i] # pick cluster index sequantly n_j = np.sum(np.nonzero(dp_cluster == (j + 1))) # count number of picked cluster if i_cluster == (j + 1): # if picked cluster is the same as current cluster n_j = (n_j - 1) # if picked cluster is the same as current cluster ## Doing Chinese restaurant process # calculate probability of previous cluster if np.isnan(dp_center[j].any()): lik = 0 else: lik = mvnpdf(cur_point, dp_center[j], dp_cov) cur_p = lik * n_j / (nr_data - 1 + dp_alpha) p[j] = cur_p lik_list[j] = lik ## Doing Chinese restaurant process(keep going) # calculate probability of making of new-cluster lik_new = np.max(lik_list) new_p = lik_new * dp_alpha / (nr_data - 1 + dp_alpha) p[j + 1] = new_p # (1) normalize all probability for clusters and (2) get index randomly(with multinomial distribution) p_norm = p / np.sum(p) nr_p = len(p_norm) z = np.nonzero(mnrand(1, np.reshape(p_norm, len(p_norm))) == 1) try: # current Deburging part, it will be fixed z = int(z) except: continue # decision wheter making new cluster or not if (z + 1) <= nr_cluster: dp_cluster[i] = z + 1 else: dp_cluster[i] = z + 1 nr_cluster += 1 # statement of making new-cluster # update information of cluster del dp_center for ci in range(nr_cluster): cluster_index = np.nonzero(dp_cluster == (ci + 1)) filter_data = data[cluster_index] try: center_pos = np.mean(filter_data, axis=0) except: continue if ci == 0: dp_center = center_pos else: dp_center = np.vstack((dp_center, center_pos)) print('number of current clusters : %d' % len(np.unique(dp_cluster))) print('function message : finish') print('\n') # result is (1) center of each cluster, (2) cluster index correspond to each data return dp_center, dp_cluster # end of 'DP' function ## execution function ## number of whole process : [length of data] x iteration # load data data = np.loadtxt(DATA + 'sample.csv', delimiter=',') # set 'hyper parameter' of Dirichlet Process alpha = 50 sigma = [0.5, 0.5] # set other value init_cluster = 4 # user can input any other integer bigger than 1(check code is needed) iteration = 100 # call 'DP' function # dp_center : center-point of each cluster # dp_cluster = cluster label of each data dp_center, dp_cluster = DP(data, alpha, sigma, init_cluster, iteration) ## print result of cluster (1st : raw data, 2nd : center of each cluster, 3rd : cluster label of each data) print('showing raw data(2-dimensional)'); print(data); print('\n'); print('\n') print('center point of each cluster'); print(dp_center); print('\n'); print('\n') print('given cluster of each data'); print(dp_cluster) point(data)
 WARNING: Output truncated! full_output.txt function message : start 0 iteration is doing number of current clusters : 4 1 iteration is doing number of current clusters : 4 2 iteration is doing number of current clusters : 4 3 iteration is doing number of current clusters : 4 4 iteration is doing number of current clusters : 4 5 iteration is doing number of current clusters : 4 6 iteration is doing number of current clusters : 4 7 iteration is doing number of current clusters : 4 8 iteration is doing number of current clusters : 4 9 iteration is doing number of current clusters : 4 10 iteration is doing number of current clusters : 4 11 iteration is doing number of current clusters : 4 12 iteration is doing number of current clusters : 4 13 iteration is doing number of current clusters : 4 14 iteration is doing number of current clusters : 4 15 iteration is doing number of current clusters : 4 16 iteration is doing number of current clusters : 5 17 iteration is doing number of current clusters : 5 18 iteration is doing number of current clusters : 4 19 iteration is doing number of current clusters : 4 20 iteration is doing number of current clusters : 4 21 iteration is doing number of current clusters : 4 22 iteration is doing number of current clusters : 4 23 iteration is doing number of current clusters : 4 24 iteration is doing number of current clusters : 4 25 iteration is doing number of current clusters : 4 26 iteration is doing number of current clusters : 4 27 iteration is doing number of current clusters : 4 28 iteration is doing number of current clusters : 4 ... [ 6.0713 7.6785 ] [ 7.5075 7.5547 ] [ 5.9073 8.6857 ] [ 5.5931 8.4311 ] [ 2.6977 5.9013 ] [ 4.2878 7.2245 ] [ 2.4675 5.8667 ] [ 2.5269 6.9344 ] [ 2.9182 6.9404 ] [ 2.2525 6.5458 ] [ 2.7593 7.9806 ] [ 2.4569 7.8856 ] [ 2.6884 7.2995 ] [ 3.3518 6.5285 ] [ 3.1219 6.7108 ] [ 1.5529 6.8987 ] [ 3.8122 7.4315 ] [ 2.3163 5.9903 ] [ 2.65 7.0086 ] [ 2.6037 5.469 ] [ 1.67 7.2917 ] [ 2.2802 7.7611 ] [ 2.9473 7.3843 ] [ 2.6626 6.5118 ] [ 2.9231 8.0081 ] [ 2.453 6.2581 ] [ 2.7578 7.5046 ] [ 3.7104 7.0077 ] [ 2.3597 8.2589 ] [ 1.4814 6.5148 ] [ 3.3948 6.5937 ] [ 2.1707 6.8331 ] [ 2.6787 6.9252 ] [ 2.8403 6.5456 ]] center point of each cluster [[ 5.39538976 5.67203684] [ 6.19154304 5.40182625] [ 5.56969799 5.85872367] [ 5.48545157 5.75329946] [ nan nan]] given cluster of each data [4 3 4 2 1 4 4 3 3 4 4 4 1 4 2 4 1 3 4 3 4 4 4 3 2 4 3 3 2 3 4 1 4 2 1 3 4 3 2 3 1 3 3 4 4 1 4 4 4 3 4 4 4 4 4 3 3 4 1 4 4 4 4 4 4 3 2 4 3 4 4 4 4 4 1 4 4 4 4 4 4 3 4 4 2 3 4 3 3 4 4 4 3 4 4 4 1 4 1 2 3 3 1 3 4 4 4 4 4 4 2 3 4 3 4 2 4 1 3 4 4 4 2 3 3 4 1 3 3 4 2 1 1 1 4 1 3 2 1 4 3 4 4 1 3 2 1 4 3 4 3 4 1 4 4 2 4 2 4 2 4 3 3 4 3 4 3 1 4 3 3 3 3 4 4 3 4 3 1 4 4 3 4 4 4 4 4 4 4 3 4 3 1 4 4 4 4 4 3 4 4 3 2 4 4 4 1 3 4 3 4 4 4 3 3 1 4 3 4 3 3 4 2 2 1 4 4 4 4 3 4 4 4 4 4 4 1 4 4 4 4 3 2 4 3 1 1 4 4 4 4 2 4 4 4 4 3 2 4 4 3 4 4 4 4 4 1 3 4 1 3 1 3 4 4 4 4 4 3 4 3 4 3 1 4 1 4 3 4 4 1 3 3 4 4 3 1 3 3 3] WARNING: Output truncated! full_output.txt function message : start 0 iteration is doing number of current clusters : 4 1 iteration is doing number of current clusters : 4 2 iteration is doing number of current clusters : 4 3 iteration is doing number of current clusters : 4 4 iteration is doing number of current clusters : 4 5 iteration is doing number of current clusters : 4 6 iteration is doing number of current clusters : 4 7 iteration is doing number of current clusters : 4 8 iteration is doing number of current clusters : 4 9 iteration is doing number of current clusters : 4 10 iteration is doing number of current clusters : 4 11 iteration is doing number of current clusters : 4 12 iteration is doing number of current clusters : 4 13 iteration is doing number of current clusters : 4 14 iteration is doing number of current clusters : 4 15 iteration is doing number of current clusters : 4 16 iteration is doing number of current clusters : 5 17 iteration is doing number of current clusters : 5 18 iteration is doing number of current clusters : 4 19 iteration is doing number of current clusters : 4 20 iteration is doing number of current clusters : 4 21 iteration is doing number of current clusters : 4 22 iteration is doing number of current clusters : 4 23 iteration is doing number of current clusters : 4 24 iteration is doing number of current clusters : 4 25 iteration is doing number of current clusters : 4 26 iteration is doing number of current clusters : 4 27 iteration is doing number of current clusters : 4 28 iteration is doing number of current clusters : 4 ... [ 6.0713 7.6785 ] [ 7.5075 7.5547 ] [ 5.9073 8.6857 ] [ 5.5931 8.4311 ] [ 2.6977 5.9013 ] [ 4.2878 7.2245 ] [ 2.4675 5.8667 ] [ 2.5269 6.9344 ] [ 2.9182 6.9404 ] [ 2.2525 6.5458 ] [ 2.7593 7.9806 ] [ 2.4569 7.8856 ] [ 2.6884 7.2995 ] [ 3.3518 6.5285 ] [ 3.1219 6.7108 ] [ 1.5529 6.8987 ] [ 3.8122 7.4315 ] [ 2.3163 5.9903 ] [ 2.65 7.0086 ] [ 2.6037 5.469 ] [ 1.67 7.2917 ] [ 2.2802 7.7611 ] [ 2.9473 7.3843 ] [ 2.6626 6.5118 ] [ 2.9231 8.0081 ] [ 2.453 6.2581 ] [ 2.7578 7.5046 ] [ 3.7104 7.0077 ] [ 2.3597 8.2589 ] [ 1.4814 6.5148 ] [ 3.3948 6.5937 ] [ 2.1707 6.8331 ] [ 2.6787 6.9252 ] [ 2.8403 6.5456 ]] center point of each cluster [[ 5.39538976 5.67203684] [ 6.19154304 5.40182625] [ 5.56969799 5.85872367] [ 5.48545157 5.75329946] [ nan nan]] given cluster of each data [4 3 4 2 1 4 4 3 3 4 4 4 1 4 2 4 1 3 4 3 4 4 4 3 2 4 3 3 2 3 4 1 4 2 1 3 4 3 2 3 1 3 3 4 4 1 4 4 4 3 4 4 4 4 4 3 3 4 1 4 4 4 4 4 4 3 2 4 3 4 4 4 4 4 1 4 4 4 4 4 4 3 4 4 2 3 4 3 3 4 4 4 3 4 4 4 1 4 1 2 3 3 1 3 4 4 4 4 4 4 2 3 4 3 4 2 4 1 3 4 4 4 2 3 3 4 1 3 3 4 2 1 1 1 4 1 3 2 1 4 3 4 4 1 3 2 1 4 3 4 3 4 1 4 4 2 4 2 4 2 4 3 3 4 3 4 3 1 4 3 3 3 3 4 4 3 4 3 1 4 4 3 4 4 4 4 4 4 4 3 4 3 1 4 4 4 4 4 3 4 4 3 2 4 4 4 1 3 4 3 4 4 4 3 3 1 4 3 4 3 3 4 2 2 1 4 4 4 4 3 4 4 4 4 4 4 1 4 4 4 4 3 2 4 3 1 1 4 4 4 4 2 4 4 4 4 3 2 4 4 3 4 4 4 4 4 1 3 4 1 3 1 3 4 4 4 4 4 3 4 3 4 3 1 4 1 4 3 4 4 1 3 3 4 4 3 1 3 3 3] 