# Simulating a (discretized) Cox process via binomial sampling

Let X be a Cox process (doubly-stochastic Poisson process) driven by a Poisson process with fixed intensity(rate) $$\lambda=50$$ , and choose some small time interval $$dt=0.01$$ . Is the proper way to simulate this, by letting Y be a binomial distribution with the number of trials equal to $$n=1$$ and the probability equal to $$p=dt*\lambda$$ and then drawing from $$Y$$ to populate each element of the time-discretized grid $x=[0..dt,dt..dt*2,dt*2..dt*3,..] ? There is an additional stipulation here, that the variable $$dt$$ must be chosen small enough such that $$p<1$$ http://people.math.aau.dk/~rw/Papers/cox2.RR.ps seems to be a good reference on Cox process inference. • If this is a fixed-intensity Cox process, isn't it just a Poisson process? A Cox process produces a counter$C(t)$that increments over time as a Poisson process with intensity$\lambda(t)\$, and the funny bit is just that the intensity itself is drawn randomly. Does this match the way you use the word "intensity"? Sep 24, 2019 at 12:42
• Not exactly. The author of the paper I'm reading claims to start using a Cox process and then simulates a Poisson process instead . When I simulate a cox process by drawing from an exponential whose parameter is also drawn randomly, I get a process whose mean is the same but variance is twice as large. What I mean is that if you simulate a Poisson process usually I draw the waiting times from an exponential with rate 1/lambda. To simulate a cox process, you actually draw the waiting times from an exponential with a different rate each time drawn from another exponential with rate 1/lambda
– crow
Sep 24, 2019 at 14:09
• The paper I'm trying to understand the usage of is at tel.archives-ouvertes.fr/tel-00778458v2/document if your curious
– crow
Sep 24, 2019 at 14:12
• To simulate a cox process, you actually draw the waiting times from an exponential with a different rate each time drawn from another exponential with rate 1/lambda. Why do you say this? Sep 24, 2019 at 17:14
• And if you already know how you want to simulate a Cox process, then what more do you need to know to answer your original question? Sep 24, 2019 at 17:16

Cox Processes normally do not have fixed intensity. The intensity is usually random. A Cox Process with fixed intensity is just a Poisson process. So, your question as currently phrased is asking how to sample from a Poisson process.

The method you propose is called a Bernoulli Process. It is not correct, but it is a reasonable approximation for a Poisson process provided that $$dt$$ is very small. To see that it is not exact, note that a Poisson random variable has no upper bound, but the Bernoulli Process approximation can only produce a maximum count of $$T/dt$$. In the limit as $$dt$$ goes to 0, your method becomes exact. This is nicely explained here (for count distributions with no temporal component) or here in section 2.2.5 (includes temporal component).

Edit: To sample exactly from a homogeneous Poisson process, one can sample the wait times between events as exponential with the rate parameter of the poisson process used as the rate parameter of the exponential draw. This is already discussed in many places and I don't think I can add much. One interesting detail is that OP wants to draw the sample one time-step at a time. To do this, one can "leapfrog" the time-steps with the exponential draws: continue drawing wait times until the next event falls past the end of the current interval. Just make sure to start from the most recent event, not from the beginning of the interval. This leapfrogging trick can be used in conjunction with the new answer below.

Comments clarified that by a "Cox process with fixed intensity", OP means to sample from the following model, where HPP means a homogeneous Poisson Process and IPP means an inhomogeneous Poisson process.

$$\lambda(t) \sim HPP(50)$$

$$N(t) \sim IPP(\lambda(t))$$

Sampling $$N(t)$$ can be done by first sampling $$\lambda$$ and then plugging that in to the conditional distribution of $$N(t)$$. The old answer above discusses how (and how not) to sample from a homogeneous Poisson process, so the missing piece is how to sample from an inhomogeneous Poisson process. There are nice discussions of that here and here, and once again I don't have much to add, except a couple of tricks to apply what they say.

The first method discussed by both sources is time-warping: 1) take out your rolling pin and stretch the time-axis to render the intensity constant; 2) sample from the resulting homogeneous PP; 3) transform the samples back onto the original scale. This process requires the integral of the intensity (and its inverse). For OP, the intensity will be piece-wise constant and nonnegative, so its integral will be piece-wise linear and non-decreasing. This is simple to invert, except for the pesky piece at the start where $$\lambda(t)=0$$. For that piece, the simulation method is simple: since the intensity is 0, skip it.

You can use this tactic along with the "leapfrogging" outlined above, but it's kind of tricky.

First, you have to remember the portions of both $$N(t)$$ and $$\lambda(t)$$ that you have already sampled. If you erase either and sample it over again, you'll be sampling from the wrong distribution, and you could even accidentally generate a $$\lambda(t)$$ or $$N(t)$$ that sometimes decreases. (Aside from being wrong, it will also probably render your algorithm quadratic because you'll have to backtrack every time.)

Second, you have to figure out the details of this time-warping procedure when $$\Lambda(t) \equiv \int_{t_0}^t \lambda(z)dz$$ is partially unknown. Suppose you're working on the interval from $$t_0$$ to $$t_1$$. First, sample $$\lambda(t)$$ until the most recent event is after $$t_1$$. Then calculate $$\Lambda$$ and its inverse between $$t_0$$ and $$t_1$$ (resp. $$y_0$$ and $$y_1$$, where $$y_i=\Lambda(t_i)$$). Then sample new wait times from your warped homogeneous process but start at $$y_0$$, not at 0. Transform the new events back to $$t$$-space. Once you reach an event in $$N(t)$$ that is past $$t_1$$ in $$t$$-space or equivalently past $$y_1$$ in $$y$$-space, stop and discard it, just as you would discard the final sample when sampling a homogeneous Poisson Process on a fixed interval. (This is not the only possible course of action, but it is simpler than alternatives where you would use this last event.)

Since you discarded the last sample but gave it the opportunity to fall short of $$t_1$$, you start sampling $$N(t)$$ not at its most recent event, but at $$t_1$$ (which acts as $$t_0$$ during the next iteration). Sorry for the confusion on this point; comments, edits, and formal proofs are welcome!

• I answered your question as stated, but I suspect this won't be satisfactory once I understand the question better, so I can edit or post a new answer at that point. Sep 25, 2019 at 19:51
• I think you misunderstood. The cox process is a poisson process whose intensity is another poisson process with fixed intensity. The driving process with fixed intensity drives the parent poisson process with random intensity
– crow
Sep 25, 2019 at 21:04
• That is not always true, as can be seen from section 3 examples 1 and 2 in the reference you provided. (The intensities in them are Gamma and Gaussian respectively.) But, your reply helps me understand your question, and I will soon edit my answer accordingly. Reference: people.math.aau.dk/~rw/Papers/cox2.RR.ps Sep 25, 2019 at 21:10
• Ahh yes you are right. Thank you erik, I've enjoyed the enlightening banter . Yes I was referring to the parts of that paper where he uses the types of processes that i described but in general the driving process can be almost any type of stochastic processes
– crow
Sep 25, 2019 at 21:33
• This answer now lays out one method, but there is a lot more to say on this topic, and there might be a better method for your purposes. Sep 25, 2019 at 22:07