Skip to content

Instantly share code, notes, and snippets.

@masayukeeeee
Created October 12, 2016 12:13
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save masayukeeeee/fdb330290067dcf0cb63e3a8e391cb2d to your computer and use it in GitHub Desktop.
Save masayukeeeee/fdb330290067dcf0cb63e3a8e391cb2d to your computer and use it in GitHub Desktop.
# パッケージ読み込み
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