R에서 혼합 가우스 모델 구현 (깁스 샘플링)
16022 단어 R
기사의 목적
R에서 혼합 가우스 모델의 깁스 샘플링을 구현합니다.
For문을 가능한 한 많이 사용하지 않고 빨리 학습하도록 했습니다.
참고 : 비 파라 메트릭 베이즈 포인트 프로세스 및 통계적 기계 학습 수리
목차
No.
목차
1
모델 설명
2
데이터 및 라이브러리
3
구현
4
확인
1. 모델 설명
2. 데이터 및 라이브러리
데이터는 iris의 데이터 세트를 사용합니다.
X <- iris[,1:4]
D <- ncol(X)
N <- nrow(X)
library(mvtnorm)
library(MCMCpack)
library(cluster)
3. 구현
library(mvtnorm)
library(MCMCpack)
#(1)Kを求める
K <- 3
#(2)muを乱数で初期化
set.seed(100)
mu <- matrix(rep(apply(X,2,mean),each=K)+rnorm(K*D,0,0.1), nrow=K)
#(3)
a0 <- 1
b0 <- 1
alpha0 <- 1
#(4)
pi <- rep(1/K, K)
t <- 1
#(5)(ⅰ)-(ⅳ)を繰り返す
max.iter <- 5
a <- a0 + N*D/2
for(s in 1:max.iter){
#(ⅰ)zのサンプリング
tmp <- t(apply(mu, 1, function(x) dmvnorm(X, x, diag(D)/t)))*as.vector(pi)
z <- apply(tmp, 2, function(x) which.max(rmultinom(1,3,x)))
#(ⅱ)mu(µ)のサンプリング
n <- tapply(z, z, length)
x.k <- apply(X, 2, function(x) tapply(x, z, mean))
mu <- t(apply(cbind(x.k, n), 1,
function(x) rmvnorm(1, x[D+1]*x[1:(D)]/(x[D+1]+1), diag(D)/(t*(x[D+1]+1)))))
#(ⅲ)t(τ)のサンプリング
b <- b0 + sum(unlist(apply(X, 2, function(x) tapply(x, z, function(x) (x-mean(x))^2))))/2 +
sum(as.vector(n/(2*(1+n)))*((x.k%*%t(x.k))*diag(K)))
t <- rgamma(1, a, b)
#(ⅳ)pi(π)のサンプリング
alpha <- alpha0 + n
pi <- rdirichlet(1, alpha)
}
4. 확인
왼쪽이 정답이고 오른쪽이 구현의 결과입니다.
library(cluster)
par(mfrow=c(1,2))
clusplot(X, iris[,5], color=TRUE, shade=FALSE, labels=4, lines=0)
clusplot(X, z, color=TRUE, shade=FALSE, labels=4, lines=0)
Reference
이 문제에 관하여(R에서 혼합 가우스 모델 구현 (깁스 샘플링)), 우리는 이곳에서 더 많은 자료를 발견하고 링크를 클릭하여 보았다
https://qiita.com/Tatsuki-Oike/items/d0931a80676d1bf7c94b
텍스트를 자유롭게 공유하거나 복사할 수 있습니다.하지만 이 문서의 URL은 참조 URL로 남겨 두십시오.
우수한 개발자 콘텐츠 발견에 전념
(Collection and Share based on the CC Protocol.)
No.
목차
1
모델 설명
2
데이터 및 라이브러리
3
구현
4
확인
1. 모델 설명
2. 데이터 및 라이브러리
데이터는 iris의 데이터 세트를 사용합니다.
X <- iris[,1:4]
D <- ncol(X)
N <- nrow(X)
library(mvtnorm)
library(MCMCpack)
library(cluster)
3. 구현
library(mvtnorm)
library(MCMCpack)
#(1)Kを求める
K <- 3
#(2)muを乱数で初期化
set.seed(100)
mu <- matrix(rep(apply(X,2,mean),each=K)+rnorm(K*D,0,0.1), nrow=K)
#(3)
a0 <- 1
b0 <- 1
alpha0 <- 1
#(4)
pi <- rep(1/K, K)
t <- 1
#(5)(ⅰ)-(ⅳ)を繰り返す
max.iter <- 5
a <- a0 + N*D/2
for(s in 1:max.iter){
#(ⅰ)zのサンプリング
tmp <- t(apply(mu, 1, function(x) dmvnorm(X, x, diag(D)/t)))*as.vector(pi)
z <- apply(tmp, 2, function(x) which.max(rmultinom(1,3,x)))
#(ⅱ)mu(µ)のサンプリング
n <- tapply(z, z, length)
x.k <- apply(X, 2, function(x) tapply(x, z, mean))
mu <- t(apply(cbind(x.k, n), 1,
function(x) rmvnorm(1, x[D+1]*x[1:(D)]/(x[D+1]+1), diag(D)/(t*(x[D+1]+1)))))
#(ⅲ)t(τ)のサンプリング
b <- b0 + sum(unlist(apply(X, 2, function(x) tapply(x, z, function(x) (x-mean(x))^2))))/2 +
sum(as.vector(n/(2*(1+n)))*((x.k%*%t(x.k))*diag(K)))
t <- rgamma(1, a, b)
#(ⅳ)pi(π)のサンプリング
alpha <- alpha0 + n
pi <- rdirichlet(1, alpha)
}
4. 확인
왼쪽이 정답이고 오른쪽이 구현의 결과입니다.
library(cluster)
par(mfrow=c(1,2))
clusplot(X, iris[,5], color=TRUE, shade=FALSE, labels=4, lines=0)
clusplot(X, z, color=TRUE, shade=FALSE, labels=4, lines=0)
Reference
이 문제에 관하여(R에서 혼합 가우스 모델 구현 (깁스 샘플링)), 우리는 이곳에서 더 많은 자료를 발견하고 링크를 클릭하여 보았다
https://qiita.com/Tatsuki-Oike/items/d0931a80676d1bf7c94b
텍스트를 자유롭게 공유하거나 복사할 수 있습니다.하지만 이 문서의 URL은 참조 URL로 남겨 두십시오.
우수한 개발자 콘텐츠 발견에 전념
(Collection and Share based on the CC Protocol.)
데이터는 iris의 데이터 세트를 사용합니다.
X <- iris[,1:4]
D <- ncol(X)
N <- nrow(X)
library(mvtnorm)
library(MCMCpack)
library(cluster)
3. 구현
library(mvtnorm)
library(MCMCpack)
#(1)Kを求める
K <- 3
#(2)muを乱数で初期化
set.seed(100)
mu <- matrix(rep(apply(X,2,mean),each=K)+rnorm(K*D,0,0.1), nrow=K)
#(3)
a0 <- 1
b0 <- 1
alpha0 <- 1
#(4)
pi <- rep(1/K, K)
t <- 1
#(5)(ⅰ)-(ⅳ)を繰り返す
max.iter <- 5
a <- a0 + N*D/2
for(s in 1:max.iter){
#(ⅰ)zのサンプリング
tmp <- t(apply(mu, 1, function(x) dmvnorm(X, x, diag(D)/t)))*as.vector(pi)
z <- apply(tmp, 2, function(x) which.max(rmultinom(1,3,x)))
#(ⅱ)mu(µ)のサンプリング
n <- tapply(z, z, length)
x.k <- apply(X, 2, function(x) tapply(x, z, mean))
mu <- t(apply(cbind(x.k, n), 1,
function(x) rmvnorm(1, x[D+1]*x[1:(D)]/(x[D+1]+1), diag(D)/(t*(x[D+1]+1)))))
#(ⅲ)t(τ)のサンプリング
b <- b0 + sum(unlist(apply(X, 2, function(x) tapply(x, z, function(x) (x-mean(x))^2))))/2 +
sum(as.vector(n/(2*(1+n)))*((x.k%*%t(x.k))*diag(K)))
t <- rgamma(1, a, b)
#(ⅳ)pi(π)のサンプリング
alpha <- alpha0 + n
pi <- rdirichlet(1, alpha)
}
4. 확인
왼쪽이 정답이고 오른쪽이 구현의 결과입니다.
library(cluster)
par(mfrow=c(1,2))
clusplot(X, iris[,5], color=TRUE, shade=FALSE, labels=4, lines=0)
clusplot(X, z, color=TRUE, shade=FALSE, labels=4, lines=0)
Reference
이 문제에 관하여(R에서 혼합 가우스 모델 구현 (깁스 샘플링)), 우리는 이곳에서 더 많은 자료를 발견하고 링크를 클릭하여 보았다
https://qiita.com/Tatsuki-Oike/items/d0931a80676d1bf7c94b
텍스트를 자유롭게 공유하거나 복사할 수 있습니다.하지만 이 문서의 URL은 참조 URL로 남겨 두십시오.
우수한 개발자 콘텐츠 발견에 전념
(Collection and Share based on the CC Protocol.)
library(mvtnorm)
library(MCMCpack)
#(1)Kを求める
K <- 3
#(2)muを乱数で初期化
set.seed(100)
mu <- matrix(rep(apply(X,2,mean),each=K)+rnorm(K*D,0,0.1), nrow=K)
#(3)
a0 <- 1
b0 <- 1
alpha0 <- 1
#(4)
pi <- rep(1/K, K)
t <- 1
#(5)(ⅰ)-(ⅳ)を繰り返す
max.iter <- 5
a <- a0 + N*D/2
for(s in 1:max.iter){
#(ⅰ)zのサンプリング
tmp <- t(apply(mu, 1, function(x) dmvnorm(X, x, diag(D)/t)))*as.vector(pi)
z <- apply(tmp, 2, function(x) which.max(rmultinom(1,3,x)))
#(ⅱ)mu(µ)のサンプリング
n <- tapply(z, z, length)
x.k <- apply(X, 2, function(x) tapply(x, z, mean))
mu <- t(apply(cbind(x.k, n), 1,
function(x) rmvnorm(1, x[D+1]*x[1:(D)]/(x[D+1]+1), diag(D)/(t*(x[D+1]+1)))))
#(ⅲ)t(τ)のサンプリング
b <- b0 + sum(unlist(apply(X, 2, function(x) tapply(x, z, function(x) (x-mean(x))^2))))/2 +
sum(as.vector(n/(2*(1+n)))*((x.k%*%t(x.k))*diag(K)))
t <- rgamma(1, a, b)
#(ⅳ)pi(π)のサンプリング
alpha <- alpha0 + n
pi <- rdirichlet(1, alpha)
}
왼쪽이 정답이고 오른쪽이 구현의 결과입니다.
library(cluster)
par(mfrow=c(1,2))
clusplot(X, iris[,5], color=TRUE, shade=FALSE, labels=4, lines=0)
clusplot(X, z, color=TRUE, shade=FALSE, labels=4, lines=0)
Reference
이 문제에 관하여(R에서 혼합 가우스 모델 구현 (깁스 샘플링)), 우리는 이곳에서 더 많은 자료를 발견하고 링크를 클릭하여 보았다 https://qiita.com/Tatsuki-Oike/items/d0931a80676d1bf7c94b텍스트를 자유롭게 공유하거나 복사할 수 있습니다.하지만 이 문서의 URL은 참조 URL로 남겨 두십시오.
우수한 개발자 콘텐츠 발견에 전념 (Collection and Share based on the CC Protocol.)