ノンパラメトリック・ベイズ実践編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"

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

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

ノンパラ・ベイズ 夏休み集中セミナーメモ

  • 9月1日
  • パラとノンパラの基礎概念
    • 資料1『Parametric vs Nonparametric Models』
    • パラは有限個パラメタ、ノンパラは無限個パラメタのモデル
    • 無限個パラメタのモデルとはどういうことかをわかることが大事
    • ノンパラベイズは無限個パラメタを想定しつつ、実際には無限要素を扱わないで処理する仕組み
    • 関数近似、分類、クラスタリングなどに関して、パラとノンパラの違いをディスカッション
    • ノンパラ関数近似であるガウシアン過程による非線形回帰をディスカッション
    • ガウシアン過程回帰の補助資料
  • ノンパラ検定
  • 9月17日
  • チュートリアル
    • まずは、5分でチュートリアルから何を読み取るかを決める
      • Non-parametric Bayesian approach (NPB)の手法としての特徴を理解し
      • 基本的な方法例2つを学び
      • その基礎概念を数式レベルで理解する
    • 次にイントロダクションを短時間でまとめる
    • Mixture models and clustering
      • Finite mixture modeling
      • Chinese restaurant process
        • クラスタ数を不定にするために、クラスタ数を無限大にする
        • 前項にて、クラス多数固定の場合には、クラスタ別に割り当てられる事前確率を用いていたので、クラスタ数が無限大の場合のクラスタ別事前確率が生成できれば、前の手法を拡張できる
        • クラスタ数が無限大の場合のその事前確率分布を定める、確率過程のモデルがchinese restaurant process
        • その分布を与えるにあたり、exchangeabilityという概念を用いている
      • Chinese restaurant process mixture models
        • クラスタ数は無限だが、有限サンプルに観測されるクラス多数はサンプル数で頭打ちされており、その頭打ちされたクラスタ数の出現確率は計算できるから、それに基づけば観測クラスタ数の事後分布も含めてベイズ推定は回る
    • Latent variable model and Indian buffet process
      • 多変量解析で次元削減がしたい
      • 0,1でできた行列(列数不定)を生成する確率過程であるIndian buffet processを使うことで、削減後次元数を減らしつつ、その削減次元数を不定にしながら、説明モデルの事後分布を生成することができる

ノンパラ・ベイズ 夏休み集中セミナーメモ0

  • 予定
    • 9月1日(パラとノンパラの基礎概念。ノンパラ検定)
    • 9月17日(ノンパラ・ベイズの短いチュートリアル)
    • 9月22日(長文資料のつまみ食い。R・パイソンで遊ぶ、その1)
    • 9月24日(長文資料のつまみ食い。R・パイソンで遊ぶ、その2)
  • 参加者
    • A(統計遺伝学分野 院生): 多次元雲状データセットを説明するためのノンパラ手法を研究中
    • B(外部病院研修医 医学部生時代の分野研究参加者):タイプ数不定の頻度推定〜ディリクレ過程
    • C(京大病院研修医 医学部生時代の分野セミナー参加者):ビッグデータ型のメディカルサイエンスに興味がある=探索型データマイニングとしてのノンパラベイズに親和性
    • D(統計遺伝学分野ポスドク):軌跡解析にスムージング。スムージングとしてのノンパラアプローチ
    • E(近大産婦人科 分野提供講義受講者): オミクスを中心にデータ解析全般に興味があり、ノンパラベイズも〜
    • F(呼吸器内科院生 共同学位課程院生): オミクスを含めた量的生物学・医療ビッグデータ学徒
    • G(統計遺伝学分野 教員): 分野のテーマ指導のためを含めた複数の理由で参加