>Is there any function or library for simulating
>the Non-Homogenous Poisson Process?
You can generate observations from a homogeneous Poisson
process, then transform the times.
If the intensity is given by f(t), with F(t) = \int_0^t f(u) du,
and (u1, u2, ...) are generated from a Poisson process
with constant rate 1, then solve
F(tk) = uk
for k = 1, 2, ...
The 'solveNonlinear' function in library(resample) would be
useful for this. If you can't express F in in closed form,
then a combination of 'approx' and 'CumIntegrate' could be used.
CumIntegrate <- function(f, a = 0, b = 1, nintervals = 5, ...){
# Cumulative integral. f is a function, a and b are endpoints
# return list(x,y) of two vectors,
# where x = seq(a, b, length = nintervals)
# and y = \int_a^{x} f(t) dt
#
# 10-point Gaussian on each of nintervals subintervals
# Assume that f can take a vector argument
ends <- seq(a, b, len = nintervals + 1)
h <- (b - a)/(2 * nintervals) #half-width of a single subinterval
mids <- (ends[-1] + ends[1:nintervals])/2
x <- c(0.14887433898163099, 0.43339539412924699, 0.67940956829902399,
0.86506336668898498, 0.97390652851717197)
wt <- c(0.29552422471475298, 0.26926671930999602, 0.21908636251598201,
0.149451349150581, 0.066671344308687999)
xvalues <- outer(h * c( - x, x), mids, "+")
list( x = ends,
y = c(0,h*cumsum(colSums( matrix( wt*f(c(xvalues), ...), 10)))))
}
========================================================
| Tim Hesterberg Senior Research Scientist |
| timh@insightful.com Insightful Corp. |
| (206)802-2319 1700 Westlake Ave. N, Suite 500 |
| (206)283-8691 (fax) Seattle, WA 98109-3044, U.S.A. |
| www.insightful.com/Hesterberg |
========================================================
Download S+Resample from www.insightful.com/downloads/libraries
I'll teach short courses:
Advanced Programming in S-PLUS: San Antonio TX, March 26-27, 2008.
Bootstrap Methods and Permutation Tests: San Antonio, March 28, 2008.
Links for more info are at www.insightful.com/Hesterberg
|