Your reports should be submitted to the Support O ce
Your report must be typed - not handwritten - and be no more than 5 sides of A4 in length
excluding gures. You will lose marks if you do not follow these requirements. Please make
sure you include your student ID code and lab group number on the front sheet, and NOT
your name!
1. Poisson Arrival Process
In this exercise we study random arrivals on a time line. This example is inspired
by real world phenomena. The arrivals can be something like customers entering a
store, phone calls coming into a phone exchange, radioactive particles arriving at a
counter, earthquakes occurring in a certain region, or typos you make in your R code
for this exercise sheet.
Typical simplifying assumptions are: there are no exactly identical arrival times (so
there’s never two events that happen simultaneously), their intensity is homogeneous
over time (you don’t have \slow" periods and \fast" periods), and their numbers in
disjoint time intervals are independent (knowing how many events have happened
already tells you nothing about how many will happen in future, as long as you know
the intensity). This leads to a model called Poisson arrival process with intensity
(for some positive constant ):
(i) The number N(I) of points in an interval I is a Poisson random variable with
mean length(I).
(ii) For disjoint intervals I1;:::;Ij, the numbers of arrivals N(I1);:::;N(Ij), are
mutually independent.
An equivalent description of this process involves the exponential distribution. It
is obtained by shifting the perspective from arrivals to waiting times. Instead of
counting arrivals in a given interval, we now measure inter-event times, that is, the
time between consecutive arrivals:
The distribution of waiting time W1 until the rst arrival is Exp( ). W1 and the sub-
sequent waiting times W2;W3;::: between each arrival and the next are independent,
all with the same exponential distribution.
(a) Consider the following piece of code and explain in words what each line of code
is doing.
[1] v <- rexp(100,1)
[2] vcum <- cumsum(v)
[3] count <- length(which(vcum <= 10))
[4 marks]
1
(b) Looking at the code in part (a), we can interpret the values in v to be the rst
100 inter-event times of a Poisson process with intensity = 1 and the values
in count to be the number of events that occur before time t = 10.
(i) Adjust the code so that count returns the number of events that occur
before time t = 0:0001. Include the adjusted code in your report. What
result do you expect to get for count in this case? Brie y justify your
answer. Run the code a few times to verify your answer but do not include
the R output for these runs in your report.
(ii) Assume that we want to count the number of events before time t = 1000
and we use the following code:
[1] v <- rexp(100,1)
[2] vcum <- cumsum(v)
[3] count <- length(which(vcum <= 1000))
Why would this approach not work? How would you adjust the code?
[4 marks]
(c) Write a function countevents() that, given a vector v of inter-event times
and a speci c time t, returns the number of events that occur before time t.
Note that you need to pay attention to the case where the total of the inter-
event times may be shorter than the interval. Make sure that your code returns
a warning message when such case occurs. Also, test your function on three
numerical examples; (v1;t1) = ((1;4;2;5;3);10), (v2;t2) = ((1;2;1);10) and
(v3;t3) = ((15;6;4;1);10), and include the output for these examples in your
report. [4 marks]
(d) For this part, you are asked to simulate a Poisson arrival process by using the
waiting time description given above. This involves simulating inter-event times
using rexp() and using countevents() to determine the number of counts.
To do this, write a function poissonsimulation() which calculates the number
of arrivals until time t in n such trials. The input parameters should be the
exponential parameter , the time t, and the number of trials n. The function
should return a vector of length n, where the ith element of this vector should
correspond to the number of events before time t for the ith trial. For this
section, you may assume 0:1 1, and 0 < t 1, otherwise the number
of inter-event times you need to calculate will get too high. Even with this
limitation though you will have to think carefully about how many inter-event
times you need to calculate. 1 [3 marks]
(e) Once complete, poissonsimulation() should return a vector with n elements,
each of which is a count. Run your function for = 0:5, t = 1 and n = 50000
to get a sample of n = 50000 counts and nd the proportion of each count value
1For a random variable X, the value x such that P(X x) = 0:01 is called the rst percentile of the
distribution. For 0:1 1, we can show that the rst percentile of Exp( ) is above 0.01. This implies
that with high probability inter-event times following Exp( ) (0:1 1) will be larger than 0.01. Which
in turn implies, that we can be fairly certain that a time period of length t will contain no more than 100t
events.
2
in your sample. In your report, only include the corresponding proportions and
the code you used to get the proportions. How do these proportions compare
with the probability mass function of a Poisson(0:5) distribution? [4 marks]
(f) Explore the distribution of the process you have generated with your simulation
function. Since all the values are discrete, bar plots are better than histograms.
Note you may not only invoke barplot() but also apply table() to the result
vector to get the data into the right shape. What happens if you change the
number of trials, the time or the intensity? Although you are free to experiment
as much as you like, in your report only include the bar plots for the proportions
for the following combinations of n, and t:
n t
500 0.1 0.1
50000 0.1 0.1
50000 0.1 0.5
50000 1 0.5
Brie y comment on what you observe. [3 marks]
2. Poisson Scatter
In this exercise we study randomly scattered points on a unit square. This spatial
example is inspired by real world phenomena such as positions of stars on a photo-
graphic plate, defects on a sheet of metal, positions of cells on a microscope slide or
raindrops hitting your head.
In many situation it is acceptable to make a few simplifying assumptions: there are
no points occurring more than once, the intensity of points is homogeneous across the
entire square, and that whatever happens in one area does not a ect what happens
elsewhere. This leads to a model called Poisson scatter with intensity (for some
positive constant ):
(i) The number N(B) of points in a subset B of the unit square is a Poisson random
variable with mean area(B).
(ii) For disjoint subsets B1;:::;Bj, the numbers of points N(B1);:::;N(Bj), are
mutually independent.
We now turn our attention to the computational implementation and exploration of
the random scatters described by this model.
(a) Generate 100 independent uniformly distributed points on a unit square and
visualise this in a plot using plot(). You may use the build-in R function
runif() that generates uniformly distributed points on an interval. Repeat this
a few times to get an idea of what such patterns look like. In you report, only
include the code and one of your plots (NOT the generated values). [2 marks]
(b) Write a function runifrect() that generates n independent uniformly dis-
tributed points on a rectangle [x1;x2] [y1;y2] within the unity square, that is,
with 0 x1 3
should be n and the coordinates of the rectangle. The function should return a
matrix of dimensions n 2. You may build in checks that activate warnings to
users who plug in coordinates that create non-admissible or degenerated rectan-
gles; ideally by writing a separate function so that this can be used routinely in
the functions you will use in the other parts of this exercise. [5 marks]
(c) Write a function plotrectpoints() that visualises your result within the unit
square by tracing the rectangle in red dotted lines and plotting the simulated
points in black. The input arguments of this function should be the coordinates
of the rectangle and a matrix containing the points as its rows. You may nd
the function rect() useful; note the arguments border and lty. [4 marks]
(d) Write a function poissonrect() that generates the realisations of a Poisson
random scatter with positive intensity parameter for rectangles. The input
arguments of this function should be the coordinates of the rectangle and .
The output value of this function should be a matrix that you can plug into
plotrectpoints(). Enjoy your hard work by exploring what such random
pattern looks like. Vary the rectangle sizes and locations and run through a
range of . Although you are free to experiment as much as you like, in your
report only include the corresponding plots for the following combinations:
x1 x2 y1 y2
100 0.1 0.5 0.3 0.8
100 0 0.5 0.3 0.4
1000 0 0.5 0.3 0.4
1000 0 0.5 0.3 1
Note that corresponds to the expected number of points for the whole unit
square, so you need to expect smaller numbers of points for smaller rectangles,
and for small you can easily end up with no points at all. [7 marks]
TOTAL: 40 MARKS