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 $(\*)$.



Density of Non-Homogeneous Poisson Process
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}$.


