Density of Non-Homogeneous Poisson Process

Question

Given $\lbrace Y_i\rbrace$ a non-homogenous Poisson process with mean density $\theta y^{-1}e^{-y}$ where $y>0$ $(\theta>0)$. I.e., the number of points of $\lbrace Y_i\rbrace$ in $(a,b)$ with $0<a<b\le \infty$ has a Poisson distribution with mean $\int_a ^b \theta y^{-1}e^{-y}dy$.

How can the density of $Y_j$ be expressed as $$\theta y^{-1}\frac{(\theta E_1(y))^{j-1}}{(j-1)!}e^{-\theta E_1(y)}, y>0 \tag{$*$}$$

where $E_1(y)$ is the exponential integral $$E_1(y)=\int_y ^\infty \frac{e^{-w}}{w}dw\,\,?$$ I found this in page 3 of http://www.stats.ox.ac.uk/~griff/pd.pdf

I know how to find the number of $Y_j$ in a given interval: $$\mathbb{P}(N(b)-N(a)=n)=\frac{(\int_a ^b \theta y^{-1}e^{-y}dy)^n \exp(-\int_a ^b\theta y^{-1}e^{-y}dy )}{n!}$$

Although this formula has some similarities to the desired formula, I am not sure how to derive $(*)$.

Answer 1

the desired density $\rho_j(y)d\epsilon$ is the probability that there is an event in the interval $(y,y+d\epsilon)$, with $d\epsilon$ infinitesimal, multiplied by the independent probability that there are $j-1$ events in the interval $(y,\infty)$:

$$\rho_j(y)=\left(\lim_{\epsilon\rightarrow 0}\frac{1}{\epsilon}\mathbb{P}(N(y+\epsilon)-N(y)=1)\right)\left(\lim_{b\rightarrow\infty}\mathbb{P}(N(b)-N(y)=j-1)\right)$$ $$=\theta y^{-1}e^{-y}\frac{1}{(j-1)!}\bigl(\theta E_1(y)\bigr)^{j-1} \exp\bigl(-\theta E_1(y) \bigr)$$

the formula $(*)$ in the OP is missing a factor $e^{-y}$.