t-SNEのt-分布

  • 次元削減の1手法のtSNE
  • 高次元の標本間の遠近情報を
    • ある点を中心として別の点が観察される確率として評価する
    • その確率を、観測標本全体で重み付き標準化する
    • さらに、その標準化確率行列が非対称なので対称化する
  • 低次元に埋め込むとき
    • 全く同様に低次元での対称化標準化確率行列をもたらす低次元座標を推定してもよいが
    • それをやると、それほどよい視覚化結果にならない
    • 低次元空間が狭いから
    • それを解決するために、裾野が広い分布(t-分布)を使うのがt-SNEのt
    • 裾野が広い分布は、狭い空間の遠い場所に標本を配置する確率を上げてくれるので、高次元の広々スペースを低次元に押し込めるときの場所の確保法として適当
  • 実際のtSNE低次元座標推定では、実利をとるためにいくつかのアルゴリズム上の工夫がされている
  • 参考:こちら
  • Rで書いて手続きを整理したのが以下のRmdフォーマットファイル

デンドログラムを二次元平面の上に描く

  • 二次元平面上に標本があり、それらの間に階層的クラスタリング構造が得られたとする
  • 二次元平面上に標本を配置しつつ、その平面から立ち上がるようにデンドログラムを描きたいかもしれない
library(FactoMineR)
# Compute PCA with ncp = 3
res.pca <- PCA(USArrests, ncp = 3, graph = FALSE)
# Compute hierarchical clustering on principal components
res.hcpc <- HCPC(res.pca, graph = FALSE)
plot(res.hcpc, choice = "3D.map")

f:id:ryamada22:20181019151527j:plain

球面調和関数係数も一般化逆行列で推定

  • 球面に場があったとする
  • それを球面調和関数で分解する
  • 球面調和関数は直交基底関数なので、個々の関数が場を説明する成分を積分して計算してもよいが、SPHARMというツール(こちら)では、一般化逆行列を使っている(この論文のメソッドセクション)
  • どういうことかというと、球面に多数の点があって、その点にある値v(\theta,\phi)があるとする。\theta,\phiは球面上の点の位置を角度で表したもの
  • 他方、この点には、球面調和関数ごとに値がある。もちろんそれらは\theta,\phiの関数としても表せる
  • さて。今、たくさんの点があったとき、個々の球面調和関数の重みがw_l^mだとしよう。l,mは球面調和関数のID番号のようなものである
  • この問題を次のような推定問題として考えることにする
  • ある点V_iには、従属変数の値y_iがあるが、この従属変数の値が、球面場の値とする
  • また、この点V_iには、説明変数の列X_iがある。この説明変数の値を、このV_iという球面上の点が持つ、球面調和関数の値とする
  • このようにすることで、標本数=球面上の点の数、説明変数の数=用いる球面調和関数の数、という多変量解析問題とみなせる
  • さて。説明変数(球面場)を、説明変数の一次線形和で表すというモデルにすることにする(球面調和関数同士は独立だから、このモデルで全く問題ない)
  • 『ある意味での最もよい推定値』は線形代数として求めても良い
  • Y = Xbを解くだけ
  • b=X^{-1}Yとすることにして、X^{-1}にムーアペンローズ一般化逆行列を使いましょう、と。
  • しかも、球面場の値(Y)を観測する点を固定しておくことにすれば、ムーアペンローズ逆行列は、1度だけ計算しておけば、後は、その逆行列を観察値ベクトルにかけるだけで、bは推定される
  • RのMASSパッケージのginv()関数はsvd()特異値分解をしているし、複素行列にも対応していて、以下に示すように、複素行列のときに、ユニタリをとっていたりする。
  • 一般化逆行列(Wiki)
> ginv
function (X, tol = sqrt(.Machine$double.eps)) 
{
    if (length(dim(X)) > 2L || !(is.numeric(X) || is.complex(X))) 
        stop("'X' must be a numeric or complex matrix")
    if (!is.matrix(X)) 
        X <- as.matrix(X)
    Xsvd <- svd(X)
    if (is.complex(X)) 
        Xsvd$u <- Conj(Xsvd$u)
    Positive <- Xsvd$d > max(tol * Xsvd$d[1L], 0)
    if (all(Positive)) 
        Xsvd$v %*% (1/Xsvd$d * t(Xsvd$u))
    else if (!any(Positive)) 
        array(0, dim(X)[2L:1L])
    else Xsvd$v[, Positive, drop = FALSE] %*% ((1/Xsvd$d[Positive]) * 
        t(Xsvd$u[, Positive, drop = FALSE]))
}
<bytecode: 0x00000000332756f0>
<environment: namespace:MASS>
  • ちなみに、(Z^t Z)^{-1} Z^tとして、正方行列の逆行列を計算するプロセスを使っても求められる。ginv()よりそのほうが速かった
library(MASS)
library(Matrix)
n <- 500
Z <- matrix(rnorm(n^2),n,n) + 1i * matrix(rnorm(n^2),n,n)
Zinv1 <- ginv(Z)
Zinv2 <- solve(t(Z)%*%Z) %*% t(Z)
range(Re(Zinv1-Zinv2))
range(Im(Zinv1-Zinv2))
> system.time(
+     for(i in 1:10){
+         Zinv1 <- ginv(Z)
+     }
+ )
   ユーザ   システム       経過  
     21.08       0.16      21.24 
> 
> system.time(
+     for(i in 1:10){
+         Zinv2 <- solve(t(Z)%*%Z) %*% t(Z)
+     }
+ )
   ユーザ   システム       経過  
     13.42       0.09      13.51 

ノンパラメトリック・ベイズ実践編2 回帰

  • まず、単純に線形回帰
s <- rnorm(n,0,1)
head(s)
m <- matrix(s,ncol=2)
head(m)
lm_res <- lm(m[,2] ~ m[,1])
lm_res
plot(m)

abline(lm_res)
x1 <- x2 <- seq(from=-2,to=2,length=100)

x12 <- expand.grid(x1,x2)

k <- function(p,l,x1,x2){
	exp(-2 * (sin(pi / p * sqrt((x1-x2)^2) / l) ^ 2))
}
d <- function(x1,x2){
	sqrt((x1-x2)^2)
}
krbf <- function(l,x1,x2){
	exp(-1 / 2 *d(x1 / l, x2 / l)^2)
}
p <- 2
l <- 1

k.val <- rep(0,length(x12[,1]))
for(i in 1:length(k.val)){
	k.val[i] <- krbf(l,x12[i,1],x12[i,2])
}

image(matrix(k.val,length(x1),length(x2)))

k(x_i, x_j) = exp(-1 / 2 d(x_i / length_scale, x_j / length_scale)^2)
  • 参加メンバーの方のパイソンコード(ipynb)も貼らせていただきます
import numpy as np
import matplotlib.pyplot as plt

x1 = np.linspace(-2,2,100)
x2 = np.linspace(-2,2,100)

x1_, x2_ = np.meshgrid(x1, x2)

def d(x,y):
    return(np.sqrt((x-y)**2))

def sin_kernel(x,y,periodicity,length_scale):
    return(np.exp(-2 * (np.sin(np.pi / periodicity * d(x, y)) / length_scale) ** 2))

out_sin = sin_kernel(x1_,x2_,1,1)
plt.pcolormesh(x1_, x2_, out_sin, cmap="spring")
plt.colorbar()

def rbf_kernel(x,y,l):
    return(np.exp((-1/2) * d(x/l,y/l) **2 ))

out_rbf = rbf_kernel(x1_,x2_,1)

plt.pcolormesh(x1_,x2_,out_rbf,cmap="spring")
plt.colorbar()

ノンパラメトリック・ベイズ実践編

  • Non-parametric bayesian clustering
    • Data set simulation
n <- 1000
d <- 4
k <- 5
s <- sample(1:k,n,replace=TRUE)

m <- matrix(rnorm(d*k,sd=6),k,d)

X <- matrix(0,n,d)
for(i in 1:k){
  tmp <- which(s == i)
  r <- matrix(rnorm(d * length(tmp)))
  X[tmp,] <- r
  for(j in 1:d){
    X[tmp,j] <- X[tmp,j] + m[i,j]
  }
}
plot(as.data.frame(X))[f:id:ryamada22:20180922131657j:plain]

f:id:ryamada22:20180922131657j:plain
point cloud

    • 乱数を使ってデータセットを作っているので、モノクロプロットとカラープロットとでポイントクラウドが異なっていることに注意
  • まず、ノンパラだけどベイズでない方法
library(clues)
out <- clues(X)
plot(out)

f:id:ryamada22:20180922133645j:plain

machine-learning.hatenablog.com
を使っているらしく、速い!

# coding: utf-8

# In[1]:


# In[2]:


from sklearn.mixture import BayesianGaussianMixture
from sklearn.datasets import *
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd


# In[3]:


plt.style.use("seaborn")
np.random.seed(42)


# In[4]:


X,y = make_blobs(1000,n_features=4, centers = 5)


# In[13]:


dpgmm = BayesianGaussianMixture(n_components=1000)
dpgmm.fit(X)
y_pred = dpgmm.predict(X)


# In[14]:


X_df = pd.DataFrame(X,columns=["var1","var2","var3","var4",])
X_df["label"] = y
X_df["label"] = X_df["label"].astype("str")
X_df["prediction"] = y_pred
X_df["prediction"] = X_df["prediction"].astype("str")


# In[15]:


sns.pairplot(X_df, vars = ["var1","var2","var3","var4",],hue = "label")
plt.show()


# In[16]:


sns.pairplot(X_df, vars = ["var1","var2","var3","var4",],hue = "prediction")
plt.show()

f:id:ryamada22:20180922135544p:plain

install.packages("CRPClustering")
library(CRPClustering)
result <- crp_gibbs(as.matrix(X), mu=c(0,0,0,0), sigma_table=14, alpha
=0.3, ro_0 =0.1, burn_in=40, iteration =200)
crp_graph_2d(as.matrix(X), result)
  • 出力
[1] "Center coordinates of Clusters"
  cluster     mu_k.1      mu_k.2       mu_k.3      mu_k.4
1      33 -1.7320941  -3.9238190  -0.02343495 -0.64889908
2       8 -5.4245621  -5.9558809  -2.48482642  4.09809942
3       5 -5.7573034  -1.9062441 -11.68433652 -0.60945428
4      44  1.8258171  -0.3260943  -0.48906236  1.01262547
5       9 -0.5258223 -11.0697677 -13.77974419  2.86106750
6      30 -4.2536066  -3.3689614   3.33196781 -2.17071835
7      29 -7.5764670  -3.0255271  -4.04623234  6.15198519
8      23  3.9127942  -1.5834261  -4.18575549  0.08683102
9      34 -2.9826239  -3.4543564   1.57405714 -1.00104252
[1] "The Number of Total Clusters  :  9"
[1] "Entropy All :  2.40723270175617"
[1] "z_result has numbers of clusters."
  • じゃあ、クラスタ正規分布からずれているとどうなるか、とやってみるも、Gaussian mixtureはひとまとめにする傾向があるのに対し、中華料理店過程MCMCは分解しすぎの傾向があり
library(MASS)
dat <- rbind(mvrnorm(50, c(0,0), diag(c(1,1))),
       mvrnorm(50, c(2,2), diag(c(1,1))),
       mvrnorm(50, c(7,7), diag(c(1,1))),
       mvrnorm(50, c(8,9), diag(c(1,1))))

plot(dat,col=c(rep("red",50),rep("blue",50),rep("yellow",50),rep("black",50)))

result <- crp_gibbs(as.matrix(dat), mu=c(0,0), sigma_table=14, alpha =0.3, ro_0 =0.1, burn_in=40, iteration =200)
  • ちょっと戻って、クラスタ数推定を(ノンパラベイズではなく)やってみよう
    • クラスタ数を決め打ちで適用して、その結果を比較して、最適クラスタ数の推定値を求めるとともにその根拠を示す
    • kmeansクラスタリングをやってBICで判定、がやっぱりいいかな…
library(mclust)
Mclust(X) -> out
plot(out)
  • 今日のまとめとして高次元データ用のクラスタリング法のレビューを読んでみる(こちら)
    • parsimonious gaussian mixture model をRのpgmm パッケージでやってみる
    • 階層的クラスタリングをしつつ、どこで切るのがよいか教えてくれる
> str(gaelle.bclust)
List of 14
 $ merge          : num [1:54, 1:2] -25 1 -12 -21 -8 5 2 4 -32 9 ...
 $ height         : num [1:54] 0 14.3 26 37.8 49.6 ...
 $ order          : num [1:55] 2 3 45 46 44 42 41 40 47 39 ...
 $ logposterior   : num [1:54] -2931 -2916 -2905 -2893 -2881 ...
 $ clust.number   : int [1:54] 54 53 52 51 50 49 48 47 46 45 ...
 $ cut            : num 741
 $ optim.alloc    : int [1:55] 1 1 1 2 2 2 2 2 2 2 ...
 $ optim.clustno  : int 7
 $ data           : num [1:55, 1:43] 0.33321 -0.00103 -0.2521 -0.17771 -0.43585 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:55] "ColWT.1" "ColWT.3" "ColWT.2" "d172.1" ...
  .. ..$ : chr [1:43] "maltose.MX1" "L.ascorbic" "glutamic.3" "raffinose2" ...
 $ repno          : int [1:55] 1 1 1 1 1 1 1 1 1 1 ...
 $ transformed.par: num [1:6] -1.84 -0.99 1.63 0.08 -0.16 -1.68
 $ labels         : chr [1:55] "ColWT.1" "ColWT.3" "ColWT.2" "d172.1" ...
 $ effect.family  : chr "gaussian"
 $ var.select     : logi TRUE
 - attr(*, "class")= chr [1:2] "bclustvs" "hclust"

インフォグラフィックとデータ視覚化

[データ視覚化][インフォグラフィック]