Created
October 12, 2016 12:13
-
-
Save masayukeeeee/fdb330290067dcf0cb63e3a8e391cb2d to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# パッケージ読み込み | |
library(MASS) | |
# 2次元正規分布からサンプルを500個抽出 | |
x <- mvrnorm(n=500, mu=c(10, (pi/2)), Sigma=matrix(c(50,1,1,0.025), 2, 2) ) | |
# シミュレーション用の状態方程式 | |
system1 <- function(x){ y <- c(x[1]*cos(x[2]), x[1]*sin(x[2])); return(t(data.frame(y=y))) } | |
# 観測値ベクトルの行列を作成 | |
y <- matrix(0, 500, 2) | |
for(i in 1:500){ | |
y[i,] <- t(system1(x[i,])) | |
} | |
# シグマポイントを計算する関数を呼び出し | |
app.U <- matrix(0, 500, 2) # U変換による近似の値を保存する行列を作成 | |
# 先ほど作ったU変換の関数を利用 | |
for(i in 1:500){ | |
tmp1 <- UT(f=system1, x=x[i,], dim=2, P=var(x), k=1, method="chol") # シグマポイントの計算結果をtmp1に保存 | |
app.U[i,] <- tmp1$ym # 近似結果をapp.Uに保存 | |
} | |
# デバイスを分割 | |
par(mfrow=c(1,2)) | |
# それぞれの共分散行列の標準偏差に対応する楕円をプロット | |
tmp <- ellipse(cov(x), centre=colMeans(x), level=0.95) | |
plot(x[,1],x[,2], xlim=c(min(c(tmp[,1], x[,1])), max(c(tmp[,1], x[,1]))), ylim=c( min(c(tmp[,2], x[,2])), max(c(tmp[,2], x[,2])) ), ylab="", xlab="", main="multivaliate normal distribution") | |
points(colMeans(x)[1], colMeans(x)[2], pch=19, col=4, cex=2) | |
polygon(ellipse(cov(x), centre=colMeans(x), level=0.95), lty=1, lwd=1, border=4) | |
tmp1 <- ellipse(cov(app.U), centre=colMeans(app.U), level=0.95) | |
tmp2 <- ellipse(cov(y), centre=colMeans(y), level=0.95) | |
xlim=c( min(c(tmp1[,1], tmp2[,1], y[,1])), max(c(tmp1[,1], tmp2[,1], y[,1]) ) ) | |
ylim=c( min(c(tmp1[,2], tmp2[,2], y[,2])), max(c(tmp1[,2], tmp2[,2], y[,2]) ) ) | |
plot(y[,1],y[,2], xlim=xlim, ylim=ylim, ylab="", xlab="", main="nonlinear function") | |
points(colMeans(y)[1], colMeans(y)[2], pch=19, col=4, cex=2) | |
polygon(ellipse(cov(y), centre=colMeans(y), level=0.95), lty=1, lwd=1, border=4) | |
points(colMeans(app.U)[1], colMeans(app.U)[2], pch=19, col=2, cex=1) | |
polygon(ellipse(cov(app.U), centre=colMeans(app.U), level=0.95), lty=2, lwd=1, border=2) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment