文書の過去の版を表示しています。
# 1. Install R # 2. Run R # 3. Copy the lines between "From Here" and "To Here" # 4. Paste those lines to R ##### From Here ##### # A function to plot the sample path of # the given event history. plot.counting.process <- function(event.history, origin=c(0,0), printing=FALSE, add=FALSE) { n.events <- length(event.history) if(printing==TRUE) print(length(event.history)) time.init <- origin[1] time.history <- append(time.init, time.init) time.history <- append(time.history, sort(append(event.history, event.history))) if(printing==TRUE) print(length(time.history)) count.init <- origin[2] count.history <- append(count.init, sort(append(c(1:n.events), c(1:n.events) ) )-1+count.init) count.history <- append(count.history, n.events+count.init) if(printing==TRUE) print(length(count.history)) if(add==TRUE) { lines(time.history, count.history) } else { plot(time.history, count.history, type="l") } } # A function to plot the sample path of # the martingale calculated from the event history # and the compensator. plot.counting.process.martingale <- function(event.history, compensator, origin=c(0,0), printing=FALSE, add=FALSE) { n.events <- length(event.history) if(printing==TRUE) print(length(event.history)) time.init <- origin[1] time.history <- append(time.init, time.init) time.history <- append(time.history, sort(append(event.history, event.history))) if(printing==TRUE) print(length(time.history)) count.init <- origin[2] count.history <- append(count.init, sort(append(c(1:n.events), c(1:n.events) ) )-1+count.init) count.history <- append(count.history, n.events+count.init) martingale <- count.history - compensator(time.history) if(printing==TRUE) print(length(count.history)) if(add==TRUE) { lines(time.history, martingale, col="grey") } else { plot(time.history, martingale, type="l") } } # A function to count the sample path of # the martingale calculated from the event history # and the compensator. count.counting.process <- function(event.history, time.to.count) { n.time <- length(time.to.count) if( n.time > 1 ) { c.events <- NULL for( i in c(1:n.time) ) { c.events <- append(c.events, sum(event.history<=time.to.count[i])) } } else { c.events <- sum(event.history<=time.to.count) } return(c.events) } count.martingale <- function(event.history, compensator, time.to.count) { n.time <- length(time.to.count) if( n.time > 1 ) { c.events <- NULL for( i in c(1:n.time) ) { c.events <- append(c.events, sum(event.history<=time.to.count[i])) } } else { c.events <- sum(event.history<=time.to.count) } m.events <- c.events - compensator(time.to.count) return(m.events) } ##### To Here ##### # Generate "an event history" # The number of targets, each has an event # during its lifetime. n.target <- 10 n.target <- 100 n.target <- 1000 # I offer three underlying models for the counting process. # We restrict ourselves to the cases where each target has # an event once at most. # Weibull distribution: # with shape parameter 2 and scale parameter 1/gamma(1+1/2) # so as to we have E[X] = 1 wb.shape <- 2 wb.scale <- 1/gamma(1+1/2) generator.x <- function(n) { return(rweibull(n=n.target, shape=wb.shape, scale=wb.scale)) } compensator.x <- function(x) { return(n.target * pweibull(x, shape=wb.shape, scale=wb.scale)) } # Exponential distribution: # with rate 1 # so as to we have E[X] = 1 exp.rate <- 1 generator.x <- function(n) { return(rexp(n=n.target, rate=exp.rate)) } compensator.x <- function(x) { return(n.target * pexp(x, rate=exp.rate)) } # Gamma distribution: # with shape parameter 2 and scale parameter 1/2 # so as to we have E[X] = 1 gamma.shape <- 2 gamma.scale <- 1/2 generator.x <- function(n) { return(rgamma(n=n.target, shape=gamma.shape, scale=gamma.scale)) } compensator.x <- function(x) { return(n.target * pgamma(x, shape=gamma.shape, scale=gamma.scale)) } # Simulations # The number of repetitions n.simulation <- 10 n.simulation <- 100 n.simulation <- 1000 # Generate event histories with "generator.x" function event.history.data <- NULL for( i in c(1:n.simulation) ) { event.history.data <- rbind(event.history.data, generator.x(n.target)) } # You can check the sample mean and the sample standard deviation # by running the following two lines. mean(event.history.data) sqrt(var(c(event.history.data))) # Plot event histories and the compensator plot(c(0, max(event.history.data)), c(0, n.target), xlab="event time", ylab="event counts", type="n") for( i in c(1:dim(event.history.data)[1]) ) { plot.counting.process(event.history.data[i,], printing=FALSE, add=TRUE) } plot(compensator.x, xlim=c(0, n.target), add=TRUE) # Plot martingales plot(c(0, max(event.history.data)), c(-10, 10), xlab="event time", ylab="martingales", type="n") cp.clt <- NULL for( i in c(1:dim(event.history.data)[1]) ) { plot.counting.process.martingale(event.history.data[i,], compensator.x, printing=FALSE, add=TRUE) cp.clt <- rbind(cp.clt, count.counting.process(event.history.data[i,], c(1,2,3))) } lines(c(0, max(event.history.data)), c(0,0)) hist(cp.clt[,1]) hist(cp.clt[,2]) hist(cp.clt[,3]) pairs(cp.clt) # If you want to export write.csv(event.history.data, "event-history-data.csv") write.csv(cp.clt, "counts-at-1-2-3.csv")
Quizes
- Observe how martingale central limit theorems work.
- Guess the relationships among means, variances, sample sizes, and simulation times.