[R] Weibull point process

[Torbjørn Ergon]

 On Thu, 15 Sep 2011 22:12:54 +1200, Rolf Turner 
 <rolf.turner at 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.turner at 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