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

Size: px

Start display at page:

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

のぶあきとべ
5 years ago
Views:

1 土木学会応用力学委員会逆問題小委員会逆問題スプリングスクール R 言語による演習 206 年 3 月 5 日 9: ~2: 講師 : 大竹雄 ( 新潟大学 ) : 山本真哉 ( 清水建設 ) : 西村伸一 ( 岡山大学 )

2 内容第一部 : プログラミングの基本的理解 (70min) R 言語とは? R 言語のインストールプログラミングのための基礎簡単な例題演習と理解第二部 : 逆問題の基礎的な例題 (80min) 基礎例題に即したプログラムミング基礎演習粒子フィルタープログラミング基礎演習

3 内容第一部 : プログラミングの基本的理解 R 言語とは? R 言語のインストールプログラミングのための基礎簡単な例題演習と理解第二部 : 逆問題の基礎的な例題例題に即したプログラムミング演習粒子フィルタープログラミング基礎演習 2

4 R 言語とは? R 言語 ( アールげんご ) はオープンソースでフリーソフトウェアの統計解析向けプログラミング言語及びその開発実行環境である R 言語はニュージーランドのオークランド大学の Ross Ihaka と Robert Gentleman により作られた現在では R Development Core eam(s 言語開発者である John M. Chambers も参画 R Project Contributors) によってメンテナンスと拡張がなされているなお R 言語仕様を実装した処理系の呼称名はプロジェクトを支援するフリーソフトウェア財団によれば GNU R(GNU R - Free Software Directory) だが他の実装形態が存在しないため当記事では日本での慣用的呼称にならい仕様実装をまとめて適宜 R 言語や R 等と呼ぶ ([Wikipedia] accessed 204/04/8) 3

5 R 言語のインストール筑波大学ダウンロードサイト : おすすめHP: R-tips より高度な問題でわからないなら, R 言語と検索 4

6 プログラミングのための基礎 : データ生成とにかく打ち込んでみましょう!! # 四則演算 +3 < *5 2/8 # ベクトルデータの生成 <- :5 <- c(:5) <- seq(from=-3,to=3,by=0.5) # このようにもできる < *c(:3) # 上の代替的方法 ( いろいろあります ) 5

7 プログラミングのための基礎 : データ生成とにかく打ち込んでみましょう!! # 繰り返し (rep) rep(:4,4) rep(:4,c(2,2,2,2)) # rep = repeate rep(:4,rep(2,4)) rep(:4,:4) <- :8 dim() <- c(2,4) dim() <- c(4,2) # dim は, 行列を生成する関数である # 行列指定によるデータの生成 <- matri(:8,2,4,byrow=f) # matriも列を生成する関数である <- matri(:8,2,4,byrow=) 6 <- matri(:6,ncol=4)

8 プログラミングのための基礎 : データ生成とにかく打ち込んでみましょう!! # ベクトルの連結による行列の生成 <- cbind(c(,2,3),c(4,5,6)) #cbindは, 二つ以上のベクトルを列として行列生成 y <- rbind(c(,4),c(2,5),c(3,6)) #rbindは, 二つ以上のベクトルを行として行列生成 y # cbind = colum bind, rbind = row bind # ベクトルと行列の演算 z <- + y # 演算子 (+ - * / ^2) は, 要素に適用される. z - y * y / y ^2 7

9 プログラミングのための基礎 : データ生成とにかく打ち込んでみましょう!! # 行列積は %*% <- matri(:4,2,2,byrow=) z <- c(5,5) ; z <- %*% t() v <- %*% z v t() # 行列の転置 inv <- solve() # 逆行列の計算も簡単 inv %*% inv 8

10 プログラミングのための基礎 : データ生成とにかく打ち込んでみましょう!! # 便利な機能 diag() # 対角要素の取り出し ncol() # 列の数 number of columbes nrow() # 行の数 number of rows solve() # 逆行列の計算 solve(,z) # 係数行列, 定数ベクトルzとする連立方程式の解 col() # 列番号を要素とする行列の生成 row() # 行番号を要素とする行列の生成 svd() # 特異値分解 9

11 プログラミングのための基礎 : データ生成とにかく打ち込んでみましょう!! # データの追加 : データの後ろに, データを追加する score <- c(45, 46, 28, 64, 45, 64, 88, 32, 54, 65, 59, 76, 83) score newscore <- c(score,c(57, 69, 56, 72, 95)) newscore newscore<- append(score,c(57, 69, 56, 72, 95)) newscore # データの挿入 : 番目と2 番目のデータの間に, 新しいデータを挿入する newscore2 <- c(score[:0],c(57, 69, 56, 72, 95),score[: length(score)]) newscore2 newscore2 <- append(score,c(57,69,56,72,95),after=0) newscore2 0

12 プログラミングのための基礎 : データ生成とにかく打ち込んでみましょう!! # データの削除 : 番目から5 番目のデータを削除する newscore3 <- newscore2[-:-5] newscore3 score # データの置換 : 番目のデータを置換する. 番目から3 番目のデータを置換する. newscore3[] <- 88 newscore3 newscore3[0:3] <- c(88,95,75,00) newscore3

13 プログラミングのための基礎 : グラフィック関係まずは, 一気に打ち込んで, コマンドの意味を考えてみてください # 例題 : グラフプロット n <- 20 #?? r <- 0.9 #?? r2 < #?? <- rnorm(n) #?? y <- r* + sqrt(-r^2)*rnorm(n) #?? 2 <- rnorm(n) #?? y2 <- r2*2 + sqrt(-r2^2)*rnorm(n) #?? plot(,y,lim=c(-4,4),ylim=c(-4,4),lab="",ylab="y") points(2,y2,pch=9) #pchはプロットの種類番号 segments(-5,0,5,0) segments(0,-5,0,5) #?? 2

14 プログラミングのための基礎 : グラフィック関係 segments y y2 points segments plot -y

15 プログラミングのための基礎 : グラフィック関係それぞれのコマンドの説明 # 例題 : グラフプロット n <- 20 # 乱数を生成する数 r <- 0.9 # 相関係数 r2 < # 相関係数 2 <- rnorm(n) # 正規乱数の生成 y <- r* + sqrt(-r^2)*rnorm(n) #に相関した正規乱数の生成 2 <- rnorm(n) # 正規乱数の生成 y2 <- r2*2 + sqrt(-r2^2)*rnorm(n) #2に相関した正規乱数の生成 plot(,y,lim=c(-4,4),ylim=c(-4,4),lab="",ylab="y") points(2,y2,pch=9) #pchはプロットの種類番号 segments(-5,0,5,0) segments(0,-5,0,5) # plot: 散布図,points: 散布図の追加,segments: 基準線の追加 4

16 プログラミングのための基礎 : グラフィック関係とにかく打ち込んでみましょう!! # 例題 : グラフプロットの種々の例題 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) # 自分でいろいろためしてみてください tet(3.5,3.5,"i",ce=2.0) tet(-3.5,3.5,"ii",ce=2.0) tet(-3.5,-3.5,"iii",ce=2.0) tet(3.5,-3.5,"vi",ce=2.0) legend(-,-2,c("rho=0.9","rho=-0.7"),pch=c(,9)) # 困ったときの対応,for help? help("plot") help("hist") 5

17 プログラミングのための基礎 : グラフィック関係 y II III rho=0.9 rho=-0.7 I VI

18 簡単な例題演習と課題 : 台形公式による数値積分 # #****** 問題数値積分 : 台形公式 *************************** # 任意関数の数値積分を考える. ここでは, 次の2つの積分を考える. # () f()= -(-)^2 + 放物線の (0,2) 区間の積分 # (2) f()= sqrt(4 - (-2)^2 ) 半径 2, 中心 (2,0) の半円の (0,4) 区間の積分 # # # ポイント # FortranのようなDo ループはいりません # この問題を解くことにより, # R 言語の特性が理解できます # F() F(s) D D s e 7

19 簡単な例題演習と課題 : 台形公式による数値積分 ()f()= -(-)^2 + 放物線の (0,2) 区間の積分 (2)f()= sqrt(4 - (-2)^2 ) 半径 2, 中心 (2,0) の半円の (0,4) 区間の積分 5 4 f() f()

20 簡単な例題演習と課題 : 台形公式による数値積分これまでの演習の復習しながら各自プログラミングしてみてくださいまずは () の例題, 時間があれば (2) へ必要なコマンド : ベクトルXの足算 sum(x) 9

21 簡単な例題演習と課題 : 台形公式による数値積分 # () の解答例 s <- 0 e <- 2 n <- 00 D <- (e-s)/n <- c(:n) <- s + (-) * D 2 <- s + * D すでにベクトルですすでにベクトルです F <- -(-)^2 +???? F2 <- -(2-)^2 + 3 行 <- (F+F2)/2*D Int <- sum() Int 20

22 簡単な例題演習と課題 : 台形公式による数値積分 # (2) の解答例 s <- 0 e <- 4 n <- 00 D <- (e-s)/n <- c(:n) <- s + (-) * D 2 <- s + * D すでにベクトルですすでにベクトルです F <- sqrt(4 - (-2)^2 )???? F2 <- sqrt(4 - (2-2)^2 ) 3 行 <- (F+F2)/2*D Int <- sum() Int 2

23 簡単な例題演習と課題 : 台形公式による数値積分 # Function 化によりすっきりさせる simpson <- function(s,e) { # fun() で定義された一変数関数を,(s,e) の間で数値積分する. n <- 00 D <- (e-s)/n <- c(:n) <- s + (-) * D???? 2 <- s + * D <- (fun()+fun(2))/2*d Int <- sum() Int } 22

24 簡単な例題演習と課題 : 台形公式による数値積分 # Function 化によりすっきりさせる fun <- function() { -(-)^2 + } curve(fun,0,2) simpson(0,2) fun <- function() { sqrt(4 - (-2)^2 ) } curve(fun,0,4) simpson(0,4) 23

25 プログラミングのための基礎 : グラフィック関係 # Rでは, 確率分布を表す関数はある規則により統一されている # 例えば正規分布を例に取ると, dnorm(,mean,sd) # 平均 mean, 標準偏差 sdの確率密度関数 (PDF) # に応じた確率密度を出力.dnorm() でディフォルトはN(0,) pnorm(,mean,sd) # 平均 mean, 標準偏差 sdの確率分布関数 (CDF)に応じた累積確率を出力. qnorm(p,mean,sd) # 平均 mean, 標準偏差 sdのquantile 関数 (CDFの逆関数)pに応じたを出力. rnorm(n,mean,sd) # 平均 mean, 標準偏差 sdのn 個の乱数を生成する. # すべての関数で, 頭文字がd,p,q,rにより, 機能が統一されている. 24

26 内容第一部 : プログラミングの基本的理解 R 言語とは? R 言語のインストールプログラミングのための基礎簡単な例題演習と理解第二部 : 逆問題の基礎的な例題基礎例題に即したプログラムミング基礎演習粒子フィルタープログラミング基礎演習 25

27 26 例題 : 優決定問題 (m>n) H z z H 変位 z, z 3, z 5 を計測し, 外力と 4 を推定する問題を考える z z z h 2 h 3 h 4 h 5 h , z 3 z 4 5 z

28 優決定問題の解 ˆ H H H z X (7.7,3.4) 正解 (8,2) 52 3 H H 7 ˆ 7 37 H H H z 3 7 ˆ X 27

29 優決定問題の解 ˆ H H H z 参考 :R 言語プログラム X (7.7,3.4) 正解 (8,2) 52 #H マトリックスの生成 H <- c(,,,,6,0) H <- matri(h,3,2) #Z ベクトルの生成 Z <- c(0,52,79) # 求解 X <- solve( t(h) %*% H ) %*% t(h) %*% Z # 解の表示確認 X X 28

30 優決定問題の解 ˆ H H H z 参考 :R 言語プログラム >> 前頁の続き X (7.7,3.4) 正解 (8,2) 52 F <- function(){ } F2 <- function(){ - 6* + 52 } F3 <- function(){ - 0* + 79 } <- c(0:00)/0 yy <- F() plot(,yy,type="l",lab="x",ylab="x4",col=2) yy <- F2() points(,yy,type="l",lty=2,col=2) yy <- F3() points(,yy,type="l",lty=3,col=2) points(x[2],x[],ce=2) points(8,2,ce=2,pch=9) X 29

31 優決定問題の解 ˆ H H H z ############################### # コンター図 >> 前頁の続き N <- 00 N2 <- 00 X (7.7,3.4) 正解 (8,2) 52 J <- matri(:n*n2,n,n2) X <- 0 for( i in :N) { X <- X + 0. X2 <- 0 for( j in :N2) { X2 <- X XX <- c(x2, X) J[i,j] <- t(z - H %*% XX ) %*% ( Z - H %*% XX ) } } <- seq(from=0.,to=0,by=0.) y <- seq(from=0.,to=0,by=0.) X contour(,y,j,add=,levels=seq(0, 0000, by=50),drawlabels = F,col=grey(0.5)) 30

32 計測データの信頼度を考慮 : 重み付き最小二乗法 X4 変位 z 3 の信頼度が低いとし, その重みを0.2とする. H WH ˆ H Wz 0 0 W X 正解 (8,2) (7.7,3.4) (7.6,2.7) 2.2 H WH 2.2 ˆ H WH H Wz ˆ

33 計測データの信頼度を考慮 : 重み付き最小二乗法変位 z 3 の信頼度が低いとし, その重みを 0.2 とする. ˆ H WH H Wz 参考 :R 言語プログラム X 正解 (8,2) (7.7,3.4) (7.6,2.7) >> 前頁から続き # データ数のカウント n <- length(z) # 観測の信頼度 ( 重み ) ベクトル W <- matri(0,n,n) diag(w) <- c(,0.2,) # 求解 X <- solve( t(h) %*% W %*% H ) %*% t(h) %*% W %*% Z X 32

34 計測データの信頼度を考慮 : 重み付き最小二乗法変位 z 3 の信頼度が低いとし, その重みを 0.2 とする. ˆ H WH H Wz 参考 :R 言語プログラム X 正解 (8,2) (7.7,3.4) (7.6,2.7) >> 前頁から続き <- c(0:00)/0 yy <- F() plot(,yy,type="l",lab="x",ylab="x4",col=2) yy <- F2() points(,yy,type="l",lty=2,col=2,lwd=2) yy <- F3() points(,yy,type="l",lty=3,col=2) points(x[2],x[],ce=2) points(8,2,ce=2,pch=9) points(x.0,x2.0,ce=2,pch=9,col=grey(0.6)) X 33

35 計測データの信頼度を考慮 : 重み付き最小二乗法 ############################### # コンター図 >> 前頁の続き ˆ H WH H Wz N <- 00 N2 <- 00 X 正解 (8,2) (7.7,3.4) (7.6,2.7) J <- matri(:n*n2,n,n2) X <- 0 for( i in :N) { X <- X + 0. X2 <- 0 for( j in :N2) { X2 <- X XX <- c(x2, X) J[i,j] <- t(z - H %*% XX ) %*% W %*% ( Z - H %*% XX ) } } <- seq(from=0.,to=0,by=0.) y <- seq(from=0.,to=0,by=0.) X contour(,y,j,add=,levels=seq(0, 0000, by=50),drawlabels = F,col=grey(0.5)) 34

36 例題 2: 粒子フィルターのプログラミング基礎劣決定問題 (m<n) 変位 z, z 5 を計測し, 外力, 3, 5 を推定する問題を考える., z 3 5, z h h2 h3 h4 h5 z z z H H 5 z

36 劣決定問題の解 : ノルム最小解 z HH H ˆ 4 5 3 0 5 6 5 3 (2,6,0) 正解 (4.2,4.6,5.

37 36 劣決定問題の解 : ノルム最小解 z HH H ˆ (2,6,0) 正解 (4.2,4.6,5.2) HH z HH H ˆ ˆ

38 37 例題 2: 粒子フィルターのプログラミング基礎 N i i t t i t t i y p y p ) ( ) ( )) ( ( )) ( ( 2 ep ) (2 ) ( i t t i t t l i t t H y R H y R y p 土木学会応用力学委員会逆問題小委員会ホームページ逆問題副読本山本真哉より引用尤度

39 例題 2: 粒子フィルターのプログラミング基礎 ################################################################# #H マトリックスの生成 H <- c(,,,6,,5) H <- matri(h,2,3) #Z ベクトルの生成 Z <- c(4,0) ################################################################# # ここからスタート # 粒子の数 N < # 事前分布 ( 一様乱数,~0) <- runif(n,,0) 2 <- runif(n,,0) 3 <- runif(n,,0) XX <- cbind(,2,3) # 変数を束ねる # 尤度の計算 W <- matri(:n,n) R <- matri(rep(0:0),2,2) # 観測誤差 diag(r) <- for(i in :N) W[i] <- ep( -/2 * t( Z - H %*% XX[i,] ) %*% solve(r) %*% ( Z - H %*% XX[i,] ) ) W <- W/sum(W) 38

40 例題 2: 粒子フィルターのプログラミング基礎 ################################################################# >> 前頁の続き # リサンプリングサブルーチン resampling <- function(w) { W.cum <- cumsum(w) # 累積和 sapply(runif(length(w)), function() which( < W.cum)[]) } No.samp <- resampling(w) No <- c(:n) ES <- data.frame(no,xx,w) ES.d2 <- ES[No.samp, ] ################################################################# 39

41 例題 2: 粒子フィルターのプログラミング基礎 ################################################################# # 描画 windows(height=0,width=0) 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(es.d2$,lim=c(0,0),breaks=br,main="",lab="",probability=) lines(c(4.2,4.2),c(0,),col=2) hist(es.d2$2,lim=c(0,0),breaks=br,main="",lab="3",probability=) lines(c(4.6,4.6),c(0,),col=2) hist(es.d2$3,lim=c(0,0),breaks=br,main="",lab="5",probability=) lines(c(5.2,5.2),c(0,),col=2) XX <- cbind(""=es.d2$, "2"=ES.d2$2, "3"=ES.d2$3) 40

42 例題 2: 粒子フィルターのプログラミング基礎 4 Density Density Density

43 例題 2: 粒子フィルターのプログラミング基礎 X X X X X X3

44 これまでの演習の復習しながら各自プログラミングしてみてください 43 例題 3: 粒子フィルターのプログラミング基礎非線形問題 ˆ 2

45 例題 3: 粒子フィルターのプログラミング基礎非線形問題 #Z ベクトルの生成 Z <- c(2,5,7) # 粒子の数 N < # 事前分布 ( 一様乱数,0~5) <- runif(n,0,5) 2 <- runif(n,0,5) XX <- cbind(,2) # 変数を束ねる # 尤度の計算 W <- matri(:n,n) R <- matri(rep(0:0),3,3) # 観測誤差 diag(r) <- RES <- Z[] - (/+/2) RES <- cbind(res,z[2] - (2/+/2)) RES <- cbind(res,z[3] - (2/+2/2)) for(i in :N) W[i] <- ep( -/2 * t( RES[i,] ) %*% solve(r) %*% ( RES[i,] ) ) W <- W/sum(W) ################################################################# No.samp <- resampling(w) No <- c(:n) ES <- data.frame(no,xx,w) ES.d2 <- ES[No.samp,] ################################################################# 44

46 例題 3: 粒子フィルターのプログラミング基礎非線形問題 ans. < ans.2 < windows(height=0,width=0) 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(,breaks=br,col=grey(0.5),border=f,probability=,lim=c(0,6),ylim=c(0,2),main="") hist(es.d2$,breaks=br,lab="",probability=,add=) lines(c(ans.,ans.),c(0,2),col=2) hist(2,breaks=br,col=grey(0.5),border=f,probability=,lim=c(0,6),ylim=c(0,2),main="") hist(es.d2$2,lim=c(0,0),breaks=br,lab="3",probability=,add=) lines(c(ans.2,ans.2),c(0,2),col=2) plot(es.d2$, ES.d2$2,lim=c(0,2),ylim=c(0,2),lab="X",ylab="X2",col=grey(0.6),pch=9) points(ans.,ans.2,col=2,pch=9) 45

47 例題 3: 粒子フィルターのプログラミング基礎非線形問題 Density Density ˆ X X X X 46

「統計数学３」

「統計数学３」関数の使い方 1 関数と引数関数の構造関数名 ( 引数 1, 引数 2, 引数 3, ) 例 : マハラノビス距離を求める関数 mahalanobis(data,m,v) 引数名を指定して記述する場合 mahalanobis(x=data, center=m, cov=v) 2 関数についてのヘルプ基本的な関数のヘルプの呼び出し? 関数名例 :?mean 例 :?mahalanobis 指定できる引数を確認する関数