Inverse Mills ratio for non normal distributions. We have the well known result of the inverse Mills ratio:
$$ \mathbb{E}[\,X\,|_{\ X > k} \,] = 
            \mu + \sigma \frac {\phi\big(\tfrac{k-\mu}{\sigma}\big)}{1-\Phi\big(\tfrac{k-\mu}{\sigma}\big)},$$
where $\phi(.)$ and $\Phi(.)$ are the PDF and CDF of the Gaussian distribution respectively.
The literature seems to confine this ratio to the Gaussian. Are there results for other, more general classes of distributions?
For instance, can one derive it for the general Pearson class class of distributions defined as $$f'(x)=-\frac{\left(a_1 x+a_0\right) }{b_2 x^2+b_1 x+b_0}f(x),$$
and express it as a function of parameters? 
$\textbf{Added:}$
We can derive the inverse Mills ratio as follows. Let $f(.)$ and $F(.)$ be the PDF and CDF, respectively, for a generic distribution with infinite support on the right and $\overline{F}(.)=1-F(.)$. Since $\mathbb{E}(X|_{\ X > k }) = \frac{\int_k^\infty x f(x) \,dx}{\int_k^\infty f(x) \,dx}$, and integrating the numerator by parts :
$$\int_k^\infty x f(x) \,dx= -k\, \overline{F}(k)+\int_k^\infty \overline{F}(x) \,\mathrm{d} x,$$
we get:
$$ \mathbb{E}(X|_{ X > k}) = 
            -k+ \frac{\int_k^\infty \overline{F}(x) \,\mathrm{d} x}{\overline{F}(k)}.$$
Now the idea is to generalize from this equation.
$\textbf{Note:}$ For a Gaussian $\int_k^\infty \overline{F}(x) \,\mathrm{d} x=\frac{1}{2} (\mu -k) \left(1-\text{erf}\left(\frac{k-\mu }{\sqrt{2} \sigma }\right)\right)+\sigma \phi(\frac{k-mu}{\sigma})$.
 A: We start with the Pearson differential equation:
$$f'(x)=-\frac{\left(a_1 x+a_0\right) }{b_2 x^2+b_1 x+b_0}f(x),$$
Define $g(x)=b_2 x^2+b_1 x+b_0$ (using a trick in Diaconis et al.(1991)).
Consider (f g)'(x), also written $(f(x) g(x))'=f'(x) g(x) +f(x) g'(x)$. We have
$$(f g)'(x)=(-a_0 + b_1 - a_1 x + 2 b_2 x) f(x)$$
We assume that the distribution has compact support of the type that $lim_{x\rightarrow +\infty} f(x) g(x)=0$.The idea is to use a probability distribution as a test function in a Schwartz distribution so $\int f' p=-\int f p'$.
Integrating on both sides:
The lhs, by parts, $\int_k^\infty f'(x) g(x) \,dx+\int_k^\infty f(x) g'(x) \, dx = -f(k)g(k)$ 
The rhs: $(b_1-a_0) \overline{F}(k)+ (2 b_2 -a_1)\int_k^\infty x f(x)\,dx$
Therefore: 
$$\int_k^\infty x f(x)\,dx=\frac{ \left(-b_0-b_1 k-b_2 k^2\right)}{2
   b_2-a_1}f(k)-\frac{(b_1-a_0) }{2 b_2-a_1}\overline{F}(k).$$
We note that $\mathbb{E}(X)=-\frac{(b_1-a_0) }{2 b_2-a_1}$, since, integrating by parts, $\int_{-\infty}^\infty (f g)'(x)=0$, and solving for $b_1-a_0+\int_{-\infty}^\infty x\, f(x) (2 b_2-a_1)\, dx=0$.
Hence $$\mathbb{E}(X|_{X>k})=\frac{ \left(-b_0-b_1 k-b_2 k^2\right)}{2
   b_2-a_1}\frac{f(k)}{\overline{F}(k)}+\mathbb{E}(X)$$
