差分

このページの2つのバージョン間の差分を表示します。

この比較画面へのリンク

次のリビジョン		前のリビジョン
次のリビジョン両方とも次のリビジョン
dm:2013 [2013/06/27 09:09] – created wataludm:2013 [2013/06/27 11:28] – watalu
行 1: 行 1:
-<code>+==== レポート課題 (Quiz) #1 (2013.06.27, Due:2013.07.10) ====
  
-# 1. Install R +1対象あたり高々1回しかイベントが発生しない状況での、マルチンゲール中心極限定理の観察。
-# 2. Run R +
-# 3. Copy the lines between "From Here" and "To Here" +
-# 4. Paste those lines to R+
  
-##### From Here ##### +=== 課題 Quizes === 
-# A function to plot the sample path of + 
-# the given event history.+Quizes 
 + 
 +  - Observe how martingale central limit theorems work. 
 +  - Guess the relationships among means, variances, sample sizes (n.target), and simulation times (n.simulation). 
 + 
 +課題 
 + 
 +  - マルチンゲール中心極限定理が作用している様子を観察(確認)しなさい。 
 +  - 背景となる分布についての関心は横においておき、分布ごとに、計数過程の平均、分散、対象数 n.target、シミュレーション回数 n.simulation の関係を調べなさい。 
 + 
 +=== 準備 (Preparation) === 
 + 
 +Instruction 
 +  - Download and install R from [[http://cran.r-project.org/|CRAN]]. ([[http://cran.ism.ac.jp/|The mirror site at ISM (Tachikawa)]] is  highly recommended for your choice.) 
 +  - Run R either by clicking the icon of R or by typing "R" on your console. 
 +  - Copy the following code. 
 +  - Paste those lines to R. 
 + 
 +手順 
 +  - Rというプログラミング言語(インタプリタ)のインストーラを[[http://cran.r-project.org/|CRAN]]からダウンロードしてインストールする (この大学からは、[[http://cran.ism.ac.jp/|統計数理研究所のミラーサイト]]を利用するのが、速くてとてもお勧め) 
 +  - Rのアイコンをクリックするか、コンソールでRと入力して、Rを起動する 
 +  - 下の囲みの中をコードをすべてコピーする 
 +  - 起動したRに貼り付けて実行させておく 
 + 
 +<code> 
 +# A function to plot the sample path of the given event history.
 plot.counting.process <- function(event.history, plot.counting.process <- function(event.history,
                                   origin=c(0,0),                                   origin=c(0,0),
行 108: 行 130:
   return(m.events)   return(m.events)
 } }
 +</code>
  
-##### To Here #####+以上のコードは次の4つの関数を含む。
  
 +  - plot.counting.process: 計数過程のサンプルパスの描画
 +  - plot.counting.process.martingale: 計数過程から作ったマルチンゲールのサンプルパス (Doob-Meyer分解でN=A+Mとあるところを、M=N-AとしてMのプロット)
 +  - count.counting.process: 特定の時点までの累積発生回数の数え上げ
 +  - count.counting.process.martingale: 未完
 +
 +=== シミュレーションデータの生成コード ===
 +
 +シミュレーションデータを生成する前に、対象数をまず設定する。
 +
 +<code>
 # Generate "an event history" # Generate "an event history"
 # The number of targets, each has an event # The number of targets, each has an event
行 117: 行 150:
 n.target <- 100 n.target <- 100
 n.target <- 1000 n.target <- 1000
 +</code>
  
 +今回の課題では「事象発生を支配する確率法則」を3つ提供する。
 +
 +  - ワイブル分布
 +  - 指数分布
 +  - ガンマ分布
 +
 +パラメータは法則ごとに異なる変数で保持しているが、データの生成はgenerator.x()、compensatorの関数はcompensator.x()と、名前を共通にしてある。
 +これはRの関数渡し(関数の引数に関数を持たせることができる)の機能を使っているため。
 +シミュレーションの実施の際には、どの法則に基づいたシミュレーションを行っているか、意識して確認しながら進める必要がある。
 +
 +<code>
 # I offer three underlying models for the counting process. # I offer three underlying models for the counting process.
 # We restrict ourselves to the cases where each target has # We restrict ourselves to the cases where each target has
行 156: 行 201:
   return(n.target * pgamma(x, shape=gamma.shape, scale=gamma.scale))   return(n.target * pgamma(x, shape=gamma.shape, scale=gamma.scale))
 } }
 +</code>
  
 +=== シミュレーションの実施 ===
  
 +<code>
 # Simulations # Simulations
 # The number of repetitions # The number of repetitions
行 186: 行 234:
 } }
 plot(compensator.x, xlim=c(0, n.target), add=TRUE) plot(compensator.x, xlim=c(0, n.target), add=TRUE)
 +</code>
  
 +マルチンゲールをプロットしてみる。ついでに特定の時点でのイベントの累積発生回数も集計させる。
 +たとえば4時点 0.5, 1.0, 1.5, 2 での集計を行うには、
 +<code>
 +time.list <- c(0.5, 1, 1.5, 2)
 +</code>
 +としてから、
 +<code>
 # Plot martingales # Plot martingales
 plot(c(0, max(event.history.data)), c(-10, 10), plot(c(0, max(event.history.data)), c(-10, 10),
行 200: 行 256:
   cp.clt <- rbind(cp.clt,   cp.clt <- rbind(cp.clt,
                   count.counting.process(event.history.data[i,],                   count.counting.process(event.history.data[i,],
-                                         c(1,2,3)))+                                         time.list))
 } }
 lines(c(0, max(event.history.data)), c(0,0)) lines(c(0, max(event.history.data)), c(0,0))
-hist(cp.clt[,1]) +colnames(cp.clt) <- time.list 
-hist(cp.clt[,2]) +</code> 
-hist(cp.clt[,3])+を実行する。これで表示されるグラフは縦軸が(-10, 10)の範囲で固定されているので、「c(-10, 10)」という箇所を変更して調整してほしい。 
 +これを例えば「c(-5,5)」とすると、(-5, 5)の範囲の描画になる。 
 + 
 +<code> 
 +par(mfrow=c(ceiling(sqrt(length(time.list))), 
 +            ceiling(sqrt(length(time.list)))) ) 
 +par(cex=0.5) 
 +for( i in c(1: length(time.list)) ) { 
 +  hist(cp.clt[,i]) 
 +} 
 +par(mfrow=c(1,1), cex=1.0) 
 +</code> 
 + 
 +<code>
 pairs(cp.clt) pairs(cp.clt)
 +</code>
  
 +その他、
  
 +<code>
 +summary(cp.clt)
 +apply(cp.clt,2,mean)
 +sqrt(apply(cp.clt,2,var))
 +</code>
 +
 +なども参考になるかもしれない。
 +
 +=== データの出力 ===
 +
 +データの処理に、RではなくExcelもしくは他の言語が便利だと思う人は、下記のようにCSVファイルに出力して、そちらで集計や作図を行ってもよい。 (If you are familier with spread sheet software like Excel or any other data analysis environments, you can export the simulation data as CSV files.)
 +
 +<code>
 # If you want to export # If you want to export
 write.csv(event.history.data, "event-history-data.csv") write.csv(event.history.data, "event-history-data.csv")
行 214: 行 298:
 </code> </code>
  
-Quizes +=== レポート提出 Report Submission ===
- +
-  - Observe how martingale central limit theorems work. +
-  - Guess the relationships among means, variances, sample sizes, and simulation times.+
  
 +電子メールにて、[[mailto:watalu_at_inf|山本]]まで送付すること。メールアドレスには適切な修正が必要だが、たぶん推測できるはず。
  
 +(Please send your report file to [[mailto:watalu_at_at_inf|me]]. You need to correct my email address.  But I believe you can guess it.)