文書の過去の版を表示しています。
# 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.