On Thu, 15 Sep 2011 22:12:54 +1200, Rolf Turner
<rolf.tur...@xtra.co.nz> wrote:
On 15/09/11 19:24, Torbjørn Ergon wrote:
On Thu, 15 Sep 2011 09:47:35 +1200, Rolf Turner
<rolf.tur...@xtra.co.nz> wrote:
On 15/09/11 07:21, Torbjørn Ergon wrote:
Dear list,
I'm looking for a function to generate (simulate) a random Weibull
point process. Can anyone help?
Cheers,
Torbjørn Ergon, University of Oslo
Do you mean a renewal process with the inter-event times having
a Weibull distribution? Should be trivial to code up, using
rweibull.
cheers,
Rolf Turner
Thanks for your reply!
Yes, in a way - but the parameters of the Weibull must change after
each event. I'm trying to simulate events that can happen multiple
times thru the life of an individual and where the hazard rate is age
dependent according to a Weibull hazard. Hence, after the first event
has occurred at time t1 the hazard rate should be h(x; t1,
shape,scale) = (shape/scale)*((x-t1)/scale)^(shape-1). This doesn't
seem trivial to do using 'rweibull' - but perhaps I'm missing
something (trivial).
Unless I am terribly confused your hazard function is just the hazard
function
of a Weibull distribution ``starting off'' from t1 rather than
starting off from 0.
So a random variable having that hazard function would have the same
distribution
as t1 + X where X is Weibull with the given shape and scale.
So you can just create the points of your process as the cumulative
sum of
a number of independent Weibull variates.
Perhaps I'm not seeing things correctly. If so, would some wiser
person please chip in and
set me straight?
cheers,
Rolf Turner
Thanks again for your reply and sorry I was so bad at explaining this!
I'm trying to simulate a process where an event can happen multiple
times and where the hazard rate declines with age of the individual (not
just with time since last event) [I wasn't quite correct in saying that
"the parameters of the Weibull must change after each event", and it
should be x+t1 in the formula I gave].
After some thinking, I have come to that I (think I) can simulate this
by discretizing and drawing events from a Bernoulli with probablity
1-exp(-integrated hazard) for each time bin - see script below. An
alternative way of simulating this would be to draw X number of values
from a Weibull distribution where X is stochastic, but I don't know what
distribution X should have.
Does this make sense?
rweibull.point.proces = function(end.time, shape, scale, delta =
end.time/1000){
event.times = NULL
for(t in seq(0,end.time, by=delta)){
p = 1 - exp(-(((t+delta)/scale)^shape - (t/scale)^shape)) #
probability of an event occuring betweem time t and t+delta
z = rbinom(1,1,p)
if(z==1) event.times = c(event.times, t+delta/2)
}
return(event.times)
}
shape = 0.5
scale = 4
end = 20
weibull.hazard = function(x,k,L) (k/L)*(x/L)^(k-1)
par(mfrow=c(1,2))
x = seq(0,20,0.1)
plot(x,weibull.hazard(x,k=shape, L=scale), type="l")
y = rweibull.point.proces(end, shape, scale)
plot(y,y ,xlim = c(0,end), ylim= c(0,end))
Cheers,
Torbjørn Ergon
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
