prml 10章 実装

# はじめに
prml(パターン認識と機械学習)の10章で紹介されている図をRｄ実装しましたので、こちらにその結果をまとめます。順次更新予定です。

# 図10.2の左 正規分布の平均場近似
```
library( mvnfast )

# 近似される正規分布の平均と精度
mu <- c( 1, 3 )
Sigma <- matrix( c( 4, 0.9*3*2, 0.9*3*2, 9 ), ncol = 2 )
Lambda <- solve( Sigma )

# サンプリングによる表示
# サンプリング数の指定とプロット
n <- 5000
result <- mvnfast::rmvn( n, mu, Sigma )
plot( result[,1], result[,2] )

####################################

# 初期値の設定
mu1 <- 0
mu2 <- 0
# mu1とmu2の変化を記録する用のベクトル
mu1_store <- mu1
mu2_store <- mu2

# 反復による推定
m <- 1 # カウンター
crite <- 1e-4 # 収束の判断基準
dif1 <- dif2 <- 1 # dif1,dif2は前の推定値との差分の絶対値を記録。初期値として1を代入
while ( ! ( dif1 < crite & dif2 < crite ) ) {
  # 値の更新
  mu1 <- mu[ 1 ] - Lambda[ 1, 1 ]^( -1 ) * Lambda[ 1, 2 ] * ( mu2 - mu[ 2 ] )
  mu2 <- mu[ 2 ] - Lambda[ 2, 2 ]^( -1 ) * Lambda[ 2, 1 ] * ( mu1 - mu[ 1 ] )

  # 更新値の記録
  mu1_store <- c( mu1_store, mu1 )
  mu2_store <- c( mu2_store, mu2 )
  
  # 差分の絶対値の計算
  dif1 <- abs( mu1_store[ m + 1 ] - mu1_store[ m ] )
  dif2 <- abs( mu2_store[ m + 1 ] - mu2_store[ m ] )
  
  # カウンターの更新
  m <- m + 1
}

# mu1, mu2の変化のプロット
plot( mu1_store, type="l" )
plot( mu2_store, type="l" )

# 近似された分布からのサンプリング
result1 <- mvnfast::rmvn( n, mu1[ length( mu1 ) ], Lambda[ 1, 1 ]^( -1 ) )
result2 <- mvnfast::rmvn( n, mu2[ length( mu2 ) ], Lambda[ 2, 2 ]^( -1 ) )

# 近似される分布のサンプリング
plot( result[ ,1 ], result[ , 2 ], col = 1, xlab = "", ylab = "" )
# 平均場近似によるサンプリング。赤で示している。
points( result1, result2, col = 2 )
```
mu1の変化過程
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/657720/448968eb-2fc7-4e37-8c44-fe05937b2fc5.png)

mu2の変化過程
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/657720/39c85d26-e79a-4205-a6c8-161f40b17364.png)

近似前と近似後のサンプリング結果の違い(黒：近似前、赤：近似後)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/657720/f0e56afa-98b6-41a2-8032-f2fd36689e79.png)

# 図10.2の右 逆のKLダイバージェンスによる近似
```
library( mvnfast )

# 近似される正規分布の平均と精度
mu <- c( 1, 3 )
Sigma <- matrix( c( 4, 0.9*3*2, 0.9*3*2, 9 ), ncol = 2 )
Lambda <- solve( Sigma )

# サンプリングによる表示
# サンプリング数の指定とプロット
n <- 5000
result <- mvnfast::rmvn( n, mu, Sigma )
plot( result[ , 1 ], result[ , 2 ] )

#######################################

# 近似された分布からのサンプリング
result1 <- mvnfast::rmvn( n, mu[ 1 ], Sigma[ 1, 1 ] )
result2 <- mvnfast::rmvn( n, mu[ 2 ], Sigma[ 2, 2 ] )

# 近似される分布のサンプリング
plot( result[ , 1 ], result[ , 2 ], col = 1, xlab = "", ylab = "" )
# 平均場近似によるサンプリング。赤で示している。
points( result1, result2, col = 2 )
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/657720/9a7bd950-b936-4695-8edd-0cdd74aee7ea.png)

# 図10.4の(a)
特にパラメータは書かれていなかったので適当にきんじしｔ
```
n_samp <- 10000 # number of sampling

#####################################################
# True distribution

# definition of parameters
shape_t = 4
scale_t = 1/5 
mu0_t <- 0
lambda0_t <- 2

# sampling
a_t <-rgamma( n_samp, shape = shape_t, scale = scale_t ) 
b_t <- numeric( length( a_t ) )
for( i in 1:length( a_t ) ){
  tau_t <- a_t[ i ]
  var_t <- ( lambda0_t * tau_t )^( -1 )
  b_t[ i ] <- rnorm( 1, mu0_t, var_t )
}

#####################################################
# The approximated distribution with initial values

# definition of parameters
shape_init = 150
scale_init = 1/100
mu_init <- 1
var_init <- 0.2

# sampling independently
a_init <- rgamma( n_samp, shape = shape_init, scale = scale_init ) 
b_init <- rnorm( n_samp, mu_init, var_init )

#####################################################
# True distribution
plot( b_t, a_t, xlim = c( -2, 2 ), ylim = c( 0, 2 ), xlab = "mu", ylab = "tau" )
# The approximated distribution with initial values
points( b_init, a_init, col = 2 )
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/657720/141ac980-5960-4f18-bb46-be489feeb1fb.png)