2

I'm given a certain seed number, lambda and k number of observations. I'm suppossed to find sum(k), divide it in unit intervals, check how many events occured in each interval, calculate the mean of events and finally calculate the absolute deviation between the experimental value I got of the mean and the theoretical expected value. In exponential distribution, lambda is the expected value. Does this still apply in each interval? Is my theoretical expected value equal to lambda? Here follows my code in R

    set.seed(4731)
    k <- 4979
    lambda <- 9
    data <- rexp(k, lambda)
    sum <- cumsum(data)
    T <- ceiling(sum[k])
    subintervals <- cut(sum, breaks = seq(0, T, 1))
    num_events <- table(subintervals)
    mean <- mean(num_events)
    expected <- lambda
    absolute_deviation <- abs(mean - expected)
    round(absolute_deviation, 4)
Anon
  • 23
  • 4

1 Answers1

2

The answer would be "close but not exactly".

What you have here is a Poisson process with rate $\lambda$, stopping after $k$ events.

If you did not stop, then the number of events in each interval of length $1$ would have have a Poisson distribution with mean $\lambda$. But you do stop, and so the distributions are censored by this, slightly reducing their means: this should be most obvious in the final interval observed but actually affecting them all.

Another issue is that the stopping time itself (you call this the overloaded sum[k] and I will call it $S$) is a random variable which has a gamma distribution with mean $\frac{k}{\lambda}$ but you round it up to $T= \lceil S \rceil$ so $T$ has a expected value which is slightly higher than $\frac{k}{\lambda}$.

$\frac{k}{S}$ has a inverse gamma distribution with mean ${\lambda}\frac{k}{k-1}$, so $\frac{k}{T}$ has a expected value which is slightly lower than ${\lambda}\frac{k}{k-1}$ but not necessarily $\lambda$. This $\frac{k}{T}$ is your average number of events per observed interval.


Your simulation gives $\frac{4979}{549}\approx 9.0692$ for mean(num_events) or $\frac kT$, close to $λ=9$ so an absolute difference of $0.0692$. You would get a different difference with a different seed, sometimes positive and sometimes negative before taking the absolute value (never exactly zero as $4979$ is not a multiple of $9$). Here is a simulation with one million cases suggesting that the expected value may be marginally below $\lambda=9$:

set.seed(2023)
k <- 4979
lambda <- 9
cases <- 10^6
T <- ceiling(rgamma(n=cases, shape=k, rate=lambda)) # number intervals
average_events_per_interval <- k/T
mean(average_events_per_interval)                 # not exactly lambda
# 8.99345
plot.ecdf(average_events_per_interval)

ecdf

Henry
  • 157,058
  • I see. What should I change in my code for it to run the way I want it to then? By the way, calling it soma[k] was just me forgetting to translate it from my native language to english for the post, it was suppossed to be sum[k], now fixed in the post. – Anon May 22 '23 at 16:56
  • 1
    Your code now runs correctly though you have overloaded sum and mean. Guessing that your language may be Portuguese, it might have been better to use soma and medio for the variables, or perhaps cumsum_times and average_events_per_interval, keeping sum and mean for the built-in functions. – Henry May 22 '23 at 21:01
  • Yes, those were the original names lol thank you very much for your help!! – Anon May 23 '23 at 16:31