Microsoft PowerPoint - Rによる演習v2.ppt [互換モード]

逆問題スプリングスクール

R言語による演習

2016年3月5日

9：○○～12：○○

逆問題小委員会

：山本真哉（清水建設）

：西村伸一（岡山大学）

土木学会応用力学委員会

講師：大竹

雄（新潟大学）

内 容

◆第一部：プログラミングの基本的理解

(70min)

•R言語とは？

•R言語のインストール

•プログラミングのための基礎

•簡単な例題演習と理解

◆第二部：逆問題の基礎的な例題(80min)

•基礎例題に即したプログラムミング基礎演習

•粒子フィルタープログラミング基礎演習

1

内 容

◆第一部：プログラミングの基本的理解

•R言語とは？

•R言語のインストール

•プログラミングのための基礎

•簡単な例題演習と理解

◆第二部：逆問題の基礎的な例題

•例題に即したプログラムミング演習

•粒子フィルタープログラミング基礎演習

R言語とは？

• R言語(アールげんご)は、オープンソースでフリーソフトウェアの統

計解析向けプログラミング言語、及びその開発実行環境である。

• R言語は、ニュージーランドのオークランド大学のRoss Ihakaと

Robert Gentlemanにより作られた。現在では、R Development

Core Team（S言語開発者であるJohn M. Chambersも参画。 R

Project Contributors）によって、メンテナンスと拡張がなされて

いる。

• なお、R言語仕様を実装した処理系の呼称名はプロジェクトを支援す

るフリーソフトウェア財団によればGNU R（GNU R - Free

Software Directory）だが、他の実装形態が存在しないため当記事

では日本での慣用的呼称にならい仕様・実装をまとめて適宜 “R言

語” や “R” 等と呼ぶ。([Wikipedia] accessed 2014/04/18)

R言語のインストール

• 筑波大学ダウンロードサイト：

http://cran.md.tsukuba.ac.jp/

• おすすめHP：“R-tips”

http://cse.naro.affrc.go.jp/takezawa/r-tips/r.html

• より高度な問題でわからないなら，“R言語

○○○”

と検索

プログラミングのための基礎：データ生成

とにかく打ち込んでみましょう！！ #四則演算 1+3 x <- 1+3 x 18-6 2*5 2/8 # ベクトル・データの生成 x <- 1:5 x x <- c(1:5) x x <- seq(from=-3,to=3,by=0.5) # このようにもできる x <- -3.5 + 0.5*c(1:13) # 上の代替的方法（いろいろあります） 5

プログラミングのための基礎：データ生成

とにかく打ち込んでみましょう！！

# 繰り返し(rep)

rep(1:4,4)

rep(1:4,c(2,2,2,2)) # rep = repeate rep(1:4,rep(2,4)) rep(1:4,1:4) x <- 1:8 dim(x) <- c(2,4) # dimは，行列を生成する関数である x dim(x) <- c(4,2) x # 行列指定によるデータの生成 x <- matrix(1:8,2,4,byrow=F) # matrixも列を生成する関数である x <- matrix(1:8,2,4,byrow=T)

プログラミングのための基礎：データ生成

とにかく打ち込んでみましょう！！ # ベクトルの連結による行列の生成 x <- cbind(c(1,2,3),c(4,5,6)) #cbindは，二つ以上のベクトルを列として行列生成 x y <- rbind(c(1,4),c(2,5),c(3,6)) #rbindは，二つ以上のベクトルを行として行列生成 y

# cbind = colum bind， rbind = row bind

# ベクトルと行列の演算 z <- x + y # 演算子(+ - * / ^2)は，要素に適用される． z x - y x * y x / y x^2 7

プログラミングのための基礎：データ生成

とにかく打ち込んでみましょう！！ # 行列積は「%*%」 x <- matrix(1:4,2,2,byrow=T) z <- c(5,5) x ; z x <- x %*% t(x) x v <- x %*% z v t(x) # 行列の転置 xinv <- solve(x) # 逆行列の計算も簡単 xinv x %*% xinv

プログラミングのための基礎：データ生成

とにかく打ち込んでみましょう！！

# 便利な機能

diag(x) # 対角要素の取り出し

ncol(x) # 列の数 number of columbes nrow(x) # 行の数 number of rows

solve(x) # 逆行列の計算 solve(x,z) # 係数行列x，定数ベクトルzとする連立方程式の解 col(x) # 列番号を要素とする行列の生成 row(x) # 行番号を要素とする行列の生成 svd(x) # 特異値分解 9

プログラミングのための基礎：データ生成

とにかく打ち込んでみましょう！！ # データの追加:データの後ろに，データを追加する score <- c(45, 46, 28, 64, 45, 64, 88, 32, 54, 65, 59, 76, 83) score newscore1 <- c(score,c(57, 69, 56, 72, 95)) newscore1 newscore1<- append(score,c(57, 69, 56, 72, 95)) newscore1 #データの挿入：11番目と12番目のデータの間に，新しいデータを挿入する

newscore2 <- c(score[1:10],c(57, 69, 56, 72, 95),score[11: length(score)]) newscore2

newscore2 <- append(score,c(57,69,56,72,95),after=10) newscore2

プログラミングのための基礎：データ生成

とにかく打ち込んでみましょう！！ # データの削除：11番目から15番目のデータを削除する newscore3 <- newscore2[-11:-15] newscore3 score # データの置換：1番目のデータを置換する． 11番目から13番目のデータを置換する． newscore3[1] <- 88 newscore3 newscore3[10:13] <- c(88,95,75,100) newscore3 11

プログラミングのための基礎：グラフィック関係

まずは，一気に打ち込んで，コマンドの意味を考えてみてください # 例題： グラフプロット n <- 20 #？？ r1 <- 0.9 # ？？ r2 <- -0.7 # ？？ x1 <- rnorm(n) # ?？ y1 <- r1*x1 + sqrt(1-r1^2)*rnorm(n) # ？？ x2 <- rnorm(n) # ？？ y2 <- r2*x2 + sqrt(1-r2^2)*rnorm(n) # ？？ plot(x1,y1,xlim=c(-4,4),ylim=c(-4,4),xlab="x",ylab="y") points(x2,y2,pch=19) #pchはプロットの種類番号 segments(-5,0,5,0) segments(0,-5,0,5) # ？？

-4 -2 0 2 4 -4 -2 0 2 4 x y

プログラミングのための基礎：グラフィック関係

13 x1-y1 x2-y2 plot points segments segments

プログラミングのための基礎：グラフィック関係

それぞれのコマンドの説明 # 例題： グラフプロット n <- 20 #乱数を生成する数 r1 <- 0.9 #相関係数① r2 <- -0.7 #相関係数② x1 <- rnorm(n) #正規乱数の生成 y1 <- r1*x1 + sqrt(1-r1^2)*rnorm(n) #x1に相関した正規乱数の生成 x2 <- rnorm(n) #正規乱数の生成 y2 <- r2*x2 + sqrt(1-r2^2)*rnorm(n) #x2に相関した正規乱数の生成 plot(x1,y1,xlim=c(-4,4),ylim=c(-4,4),xlab="x",ylab="y") points(x2,y2,pch=19) #pchはプロットの種類番号 segments(-5,0,5,0) segments(0,-5,0,5) # plot：散布図，points：散布図の追加，segments：基準線の追加

プログラミングのための基礎：グラフィック関係

とにかく打ち込んでみましょう！！ # 例題： グラフプロットの種々の例題 lines(c(0.0,-4.0),c(0.0,4.0),lty=2) lines(c(0.0,4.0),c(0.0,4.0),lty=2) #自分でいろいろためしてみてください text(3.5,3.5,"I",cex=2.0) text(-3.5,3.5,"II",cex=2.0) text(-3.5,-3.5,"III",cex=2.0) text(3.5,-3.5,"VI",cex=2.0) legend(-1,-2,c("rho=0.9","rho=-0.7"),pch=c(1,19)) # 困ったときの対応，for help? help("plot") help("hist") 15

-4 -2 0 2 4 -4 -2 0 2 4 y

I

II

III

VI

rho=0.9 rho=-0.7

プログラミングのための基礎：グラフィック関係

簡単な例題演習と課題：台形公式による数値積分

# #****** 問題 数値積分： 台形公式 *************************** # 任意関数の数値積分を考える．ここでは，次の2つの積分を考える． # (1) f(x)= -(x-1)^2 + 1 放物線の(0,2)区間の積分 # (2) f(x)= sqrt(4 - (x-2)^2 ) 半径2，中心(2,0)の半円の(0,4)区間の積分 # # # ※ポイント # FortranのようなDo ループはいりません # この問題を解くことにより， # R言語の特性が理解できます # 17 xs _xe Dx Dx F(xs) F(x)

簡単な例題演習と課題：台形公式による数値積分

x f( x ) 0.0 0.5 1.0 1.5 2.0 2.5 3.0 -3 -2 -1 0 1 x f(x ) 0 1 2 3 4 5 0 1 2 3 4 5 (1)f(x)= -(x-1)^2 + 1 放物線の(0,2)区間の積分 (2)f(x)= sqrt(4 - (x-2)^2 ) 半径2，中心(2,0)の半円の(0,4)区間の積分

簡単な例題演習と課題：台形公式による数値積分

19

これまでの演習の復習しながら

各自プログラミングしてみてください

まずは(1)の例題，時間があれば(2)へ 必要なコマンド：ベクトルXの足算 sum(X)

# (1)の解答例 xs <- 0 xe <- 2 n <- 100 Dx <- (xe-xs)/n x <- c(1:n) x1 <- xs + (x-1) * Dx x2 <- xs + x * Dx F1 <- -(x1-1)^2 + 1 F2 <- -(x2-1)^2 + 1 x <- (F1+F2)/2*Dx Int <- sum(x) Int

簡単な例題演習と課題：台形公式による数値積分

？？？？ 3行 ※すでにベクトルです ※すでにベクトルです

21 # (2)の解答例 xs <- 0 xe <- 4 n <- 100 Dx <- (xe-xs)/n x <- c(1:n) x1 <- xs + (x-1) * Dx x2 <- xs + x * Dx F1 <- sqrt(4 - (x1-2)^2 ) F2 <- sqrt(4 - (x2-2)^2 ) x <- (F1+F2)/2*Dx Int <- sum(x) Int

簡単な例題演習と課題：台形公式による数値積分

？？？？ 3行 ※すでにベクトルです ※すでにベクトルです

# Function化によりすっきりさせる simpson <- function(xs,xe) { # fun(x)で定義された一変数関数を，(xs,xe)の間で数値積分する． n <- 100 Dx <- (xe-xs)/n x <- c(1:n) x1 <- xs + (x-1) * Dx x2 <- xs + x * Dx x <- (fun(x1)+fun(x2))/2*Dx Int <- sum(x) Int }

簡単な例題演習と課題：台形公式による数値積分

？？？？

23 # Function化によりすっきりさせる fun <- function(x) { -(x-1)^2 + 1 } curve(fun,0,2) simpson(0,2) fun <- function(x) { sqrt(4 - (x-2)^2 ) } curve(fun,0,4) simpson(0,4)

簡単な例題演習と課題：台形公式による数値積分

プログラミングのための基礎：グラフィック関係

# Rでは，確率分布を表す関数はある規則により統一されている # 例えば正規分布を例に取ると， dnorm(x,xmean,xsd) ＃平均 xmean，標準偏差xsdの確率密度関数（PDF） # xに応じた確率密度を出力．dnorm(x)でディフォルトはN(0,1) pnorm(x,xmean,xsd) #平均 xmean，標準偏差xsdの確率分布関数（CDF）xに応じた累積確率を出力． qnorm(p,xmean,xsd) #平均 xmean，標準偏差xsdのquantile関数（CDFの逆関数）pに応じたxを出力． rnorm(n,xmean,xsd) #平均 xmean，標準偏差xsdのn個の乱数を生成する． #すべての関数で，頭文字がd,p,q,rにより，機能が統一されている．

内 容

◆第一部：プログラミングの基本的理解

•R言語とは？

•R言語のインストール

•プログラミングのための基礎

•簡単な例題演習と理解

◆第二部：逆問題の基礎的な例題

•基礎例題に即したプログラムミング基礎演習

•粒子フィルタープログラミング基礎演習

25

例題1：優決定問題（m>n）

Hx

z



x



H

1

z

変位z

₁

, z

₃

, z

₅

を計測し，外力x

₁

とx

₄

を推定する問題を考える．













































































4 1 5 3 1

10

1

6

1

1

1

79

52

10

x

x

z

z

z

























































79

52

10

10

1

6

1

1

1

1 4 1

x

x

1

h

h

₂

h

₃

h

₄

h

₅ 1 2 3 4 5 1 1

, z

x

z

₃

x

₄

z

₅

27

優決定問題の解

0 2 4 6 8 10 02468 1 0 X1 X4 10 4 1 x  x 52 6 ₄ 1 x  x 正解 (7.7,3.4) (8,2) 79 10 ₄ 1 x  x















137

17

17

3

T

H

H



H

H



H

z

x

ˆ



T 1 T

















































79

52

10

10

6

1

1

1

1

137

17

17

3

1





























4 1

4

.

3

7

.

7

ˆ

x

x

x



H

H



H

z

x

ˆ



T 1 T

優決定問題の解

【参考：R言語プログラム】 #Hマトリックスの生成 H <- c(1,1,1,1,6,10) H <- matrix(H,3,2) #Zベクトルの生成 Z <- c(10,52,79) #求解 X <- solve( t(H) %*% H ) %*% t(H) %*% Z #解の表示・確認 X 0 2 4 6 8 10 02468 1 0 X1 X4 10 4 1 x  x 52 6 ₄ 1 x  x 正解 (7.7,3.4) (8,2) 79 10 ₄ 1 x  x



H

H



H

z

x

ˆ



T 1 T

29

優決定問題の解

【参考：R言語プログラム】 >> 前頁の続き F1 <- function(xx){ - xx + 10 } F2 <- function(xx){ - 6*xx + 52 } F3 <- function(xx){ - 10*xx + 79 } xx <- c(0:100)/10 yy <- F1(xx) plot(xx,yy,type="l",xlab="X1",ylab="X4",col=2) yy <- F2(xx) points(xx,yy,type="l",lty=2,col=2) yy <- F3(xx) points(xx,yy,type="l",lty=3,col=2) points(X[2],X[1],cex=2) points(8,2,cex=2,pch=19) 0 2 4 6 8 10 02468 1 0 X1 X4 10 4 1 x  x 52 6 ₄ 1 x  x 正解 (7.7,3.4) (8,2) 79 10 ₄ 1 x  x



H

H



H

z

x

ˆ



T 1 T

0 2 4 6 8 10 02 4 6 8 1 0 X1 X4

優決定問題の解



H

H



H

z

x

ˆ



T 1 T ############################### #コンター図 >> 前頁の続き N1 <- 100 N2 <- 100 J <- matrix(1:N1*N2,N1,N2) X1 <- 0 for( i in 1:N1) { X1 <- X1 + 0.1 X2 <- 0 for( j in 1:N2) { X2 <- X2 + 0.1 XX <- c(X2, X1) J[i,j] <- t(Z - H %*% XX ) %*% ( Z - H %*% XX ) } } x <- seq(from=0.1,to=10,by=0.1) y <- seq(from=0.1,to=10,by=0.1) contour(x,y,J,add=T,levels=seq(0, 10000, by=50),drawlabels = F,col=grey(0.5)) 10 4 1 x  x 52 6 ₄ 1 x  x 正解 (7.7,3.4) (8,2) 79 10 ₄ 1 x  x

0 2 4 6 8 10 02468 1 0 X1 X4 31

計測データの信頼度を考慮：重み付き最小二乗法



H

WH



H

Wz

x

ˆ



T 1 T

変位z

₃

の信頼度が低いとし，その重みを0.2とする．

正解 (7.7,3.4) (8,2) (7.6,2.7)















2

.

108

2

.

12

2

.

12

2

.

2

T

WH

H























1

1

0

0

2

.

0

0

0

0

1

W



H

WH



H

Wz

x

ˆ



T 1 T

















































79

52

10

10

6

1

1

1

1

2

.

108

2

.

12

2

.

12

2

.

2

1





























4 1

7

.

2

6

.

7

ˆ

x

x

x

0 2 4 6 8 10 02468 1 0 X1 X4

計測データの信頼度を考慮：重み付き最小二乗法



H

WH



H

Wz

x

ˆ



T 1 T

変位z

₃

の信頼度が低いとし，その重みを0.2とする．

正解 (7.7,3.4) (8,2) (7.6,2.7) 【参考：R言語プログラム】 >> 前頁から続き #データ数のカウント n <- length(Z) #観測の信頼度（重み）ベクトル W <- matrix(0,n,n) diag(W) <- c(1,0.2,1) #求解 X <- solve( t(H) %*% W %*% H ) %*% t(H) %*% W %*% Z

0 2 4 6 8 10 02468 1 0 X1 X4 33

計測データの信頼度を考慮：重み付き最小二乗法



H

WH



H

Wz

x

ˆ



T 1 T

変位z

₃

の信頼度が低いとし，その重みを0.2とする．

正解 (7.7,3.4) (8,2) (7.6,2.7) 【参考：R言語プログラム】 >> 前頁から続き xx <- c(0:100)/10 yy <- F1(xx) plot(xx,yy,type="l",xlab="X1",ylab="X4",col=2) yy <- F2(xx) points(xx,yy,type="l",lty=2,col=2,lwd=2) yy <- F3(xx) points(xx,yy,type="l",lty=3,col=2) points(X[2],X[1],cex=2) points(8,2,cex=2,pch=19) points(X1.0,X2.0,cex=2,pch=19,col=grey(0.6))

0 2 4 6 8 10 0 2 468 1 0 X1 X4

計測データの信頼度を考慮：重み付き最小二乗法



H

WH



H

Wz

x

ˆ



T 1 T 正解 (7.7,3.4) (8,2) (7.6,2.7) ############################### #コンター図 >> 前頁の続き N1 <- 100 N2 <- 100 J <- matrix(1:N1*N2,N1,N2) X1 <- 0 for( i in 1:N1) { X1 <- X1 + 0.1 X2 <- 0 for( j in 1:N2) { X2 <- X2 + 0.1 XX <- c(X2, X1) J[i,j] <- t(Z - H %*% XX ) %*% W %*% ( Z - H %*% XX ) } } x <- seq(from=0.1,to=10,by=0.1) y <- seq(from=0.1,to=10,by=0.1) contour(x,y,J,add=T,levels=seq(0, 10000, by=50),drawlabels = F,col=grey(0.5))

35

例題2：粒子フィルターのプログラミング基礎

劣決定問題（m<n）

変位z

₁

, z

₅

を計測し，外力x

₁

, x

₃

, x

₅

を推定する問題を考える．

Hx

z



x



H

1

z

1 2 3 4 5 1 1

, z

x

x

₃

x

₅

, z

₅ 1

h

h

₂

h

₃

h

₄

h

₅





























































5 3 1 5 1

15

6

1

1

1

1

250

14

x

x

x

z

z

















































110

14

15

6

1

1

1

1

1 5 3 1

x

x

x

劣決定問題の解：ノルム最小解



HH



z

H

x

ˆ



T T 1 14 5 3 1x x  x 110 15 6 _{3 5} 1 x  x  x (2,6,10) 正解 (4.2,4.6,5.2)















262

22

22

3

T

HH



HH



z

H

x

ˆ



T T 1

















































250

14

262

22

22

3

15

1

6

1

1

1

₁













































₃ 1

2

.

5

6

.

4

2

.

4

ˆ

x

x

x

x

37

例題2：粒子フィルターのプログラミング基礎







_N i i t t i t t i

x

y

p

x

y

p

1

)

(

)

(











_

_

_





))

(

(

))

(

(

2

1

exp

)

2

(

1

)

(

_{t t}i T 1 _{t t}i l i t t

y

H

x

R

y

H

x

R

x

y

p



土木学会応用力学委員会逆問題小委 員会ホームページ逆問題副読本 山本真哉より引用 尤度

例題2：粒子フィルターのプログラミング基礎

################################################################# #Hマトリックスの生成 H <- c(1,1,1,6,1,15) H <- matrix(H,2,3) #Zベクトルの生成 Z <- c(14,110) ################################################################# #ここからスタート #粒子の数 N <- 10000 #事前分布（一様乱数，1～10） x1 <- runif(N,1,10) x2 <- runif(N,1,10) x3 <- runif(N,1,10) XX <- cbind(x1,x2,x3) #変数を束ねる #尤度の計算 W <- matrix(1:N,N) R <- matrix(rep(0:0),2,2) #観測誤差 diag(R) <- 1

for(i in 1:N) W[i] <- exp( -1/2 * t( Z - H %*% XX[i,] ) %*% solve(R) %*% ( Z - H %*% XX[i,] ) ) W <- W/sum(W)

39

例題2：粒子フィルターのプログラミング基礎

################################################################# >> 前頁の続き #リサンプリング・サブルーチン resampling <- function(W) { W.cum <- cumsum(W) #累積和

sapply(runif(length(W)), function(x) which(x < W.cum)[1]) } No.samp <- resampling(W) No <- c(1:N) TEST <- data.frame(No,XX,W) TEST.d2 <- TEST[No.samp, ] #################################################################

例題2：粒子フィルターのプログラミング基礎

################################################################# #描画 windows(height=10,width=10) frame() plot.new() par(mfrow=c(2,2)) # グラフの外の余白 par(oma = c(5, 5, 2, 2)) # 上段と下段のグラフ間の余白(ピッタリくっつけるときは"0") par(mar = c(5, 5, 0, 0)) br <- c(0:20)/2 hist(TEST.d2$x1,xlim=c(0,10),breaks=br,main="",xlab="x1",probability=T) lines(c(4.2,4.2),c(0,1),col=2) hist(TEST.d2$x2,xlim=c(0,10),breaks=br,main="",xlab="x3",probability=T) lines(c(4.6,4.6),c(0,1),col=2) hist(TEST.d2$x3,xlim=c(0,10),breaks=br,main="",xlab="x5",probability=T) lines(c(5.2,5.2),c(0,1),col=2)

41

例題2：粒子フィルターのプログラミング基礎

x1 De n s it y 0 2 4 6 8 10 12 0. 0 0 .1 0. 2 0 .3 0. 4 0 .5 x2 De n s it y 0 2 4 6 8 10 12 0. 0 0 .1 0. 2 0 .3 0. 4 0 .5 x3 D ens it y 0 2 4 6 8 10 12 0 .0 0 .1 0 .2 0 .3 0 .4 0 .5

例題2：粒子フィルターのプログラミング基礎

0 2 4 6 8 10 02 468 1 0 X1 X3 0 2 4 6 8 10 02 468 1 0 X1 X5 0 2 4 6 8 10 0 2468 1 0 X5

これまでの演習の復習しながら

各自プログラミングしてみてください

43

例題3：粒子フィルターのプログラミング基礎

非線形問題

































































2 1 2 1 2 1

2

2

1

2

1

1

7

5

2

x

x

x

x

x

x





























0.7142

0.5555

ˆ

2 1

x

x

x

例題3：粒子フィルターのプログラミング基礎

非線形問題

#Zベクトルの生成 Z <- c(2,5,7) #粒子の数 N <- 10000 #事前分布（一様乱数，0～5） x1 <- runif(N,0,5) x2 <- runif(N,0,5) XX <- cbind(x1,x2) #変数を束ねる #尤度の計算 W <- matrix(1:N,N) R <- matrix(rep(0:0),3,3) #観測誤差 diag(R) <- 1 RES <- Z[1] - (1/x1+1/x2) RES <- cbind(RES,Z[2] - (2/x1+1/x2)) RES <- cbind(RES,Z[3] - (2/x1+2/x2))

for(i in 1:N) W[i] <- exp( -1/2 * t( RES[i,] ) %*% solve(R) %*% ( RES[i,] ) ) W <- W/sum(W) ################################################################# No.samp <- resampling(W) No <- c(1:N) TEST <- data.frame(No,XX,W) TEST.d2 <- TEST[No.samp,] #################################################################

45

例題3：粒子フィルターのプログラミング基礎

非線形問題

ans.x1 <- 0.5555 ans.x2 <- 0.7142 windows(height=10,width=10) frame() plot.new() par(mfrow=c(2,2)) # グラフの外の余白 par(oma = c(5, 5, 2, 2)) # 上段と下段のグラフ間の余白(ピッタリくっつけるときは"0") par(mar = c(5, 5, 0, 0)) br <- c(0:20)/4 hist(x1,breaks=br,col=grey(0.5),border=F,probability=T,xlim=c(0,6),ylim=c(0,2),main="") hist(TEST.d2$x1,breaks=br,xlab="x1",probability=T,add=T) lines(c(ans.x1,ans.x1),c(0,2),col=2) hist(x2,breaks=br,col=grey(0.5),border=F,probability=T,xlim=c(0,6),ylim=c(0,2),main="") hist(TEST.d2$x2,xlim=c(0,10),breaks=br,xlab="x3",probability=T,add=T) lines(c(ans.x2,ans.x2),c(0,2),col=2) plot(TEST.d2$x1, TEST.d2$x2,xlim=c(0,2),ylim=c(0,2),xlab="X1",ylab="X2",col=grey(0.6),pch=19) points(ans.x1,ans.x2,col=2,pch=19)

例題3：粒子フィルターのプログラミング基礎

非線形問題

































































2 1 2 1 2 1

2

2

1

2

1

1

7

5

2

x

x

x

x

x

x

x1 De n s it y 0 1 2 3 4 5 6 0 .0 0 .5 1 .0 1 .5 2 .0 x2 De n s it y 0 1 2 3 4 5 6 0 .0 0 .5 1 .0 1 .5 2 .0 0.0 0.5 1.0 1.5 2.0 0 .00 .5 1 .01 .5 2 .0 X2 0.0 0.5 1.0 1.5 2.0 0 .00 .5 1 .01 .5 2 .0 X2





























0.7142

0.5555

2

ˆ

1

x

x

x