¿Es posible integrar analíticamente $x$ multiplicada por la función de densidad de probabilidad lognormal?

Question

En primer lugar, por integrar analíticamente me refiero a si existe una regla de integración para resolver esto en lugar de análisis numéricos (como las reglas trapezoidal, de Gauss-Legendre o de Simpson).

Tengo una función $\newcommand{\rd}{\mathrm{d}}f(x) = x g(x; \mu, \sigma)$ donde $$ g(x; \mu, \sigma) = \frac{1}{\sigma x \sqrt{2\pi}} e^{-\frac{1}{2\sigma^2}(\log(x) - \mu)^2} $$ es la función de densidad de probabilidad de una distribución lognormal con parámetros $\mu$ y $\sigma$ . A continuación, abreviaré la notación a $g(x)$ y utilizar $G(x)$ para la función de distribución acumulativa.

Necesito calcular la integral $$ \int_{a}^{b} f(x) \,\rd x \>. $$

Actualmente, lo hago con integración numérica utilizando el método de Gauss-Legendre. Debido a que necesito ejecutar esto un gran número de veces, el rendimiento es importante. Antes de mirar en la optimización de los análisis numéricos / otras piezas, me gustaría saber si hay alguna regla de integración para resolver esto.

Intenté aplicar la regla de integración por partes y llegué a esto, donde estoy atascado de nuevo,

1. $\int u \,\mathrm{d}v = u v - \int v \mathrm{d}u$ .

2. $u=x \implies \rd u = \rd x$

3. $\rd v = g(x) \rd x \implies v = G(x)$

4. $u v - \int v \rd x = x G(x) - \int G(x) \rd x$

Estoy atascado, ya que no puedo evaluar la $\int G(x) \rd x$ .

Esto es para un paquete de software que estoy construyendo.

Answer 1

Respuesta corta : No, no es posible, al menos en términos de funciones elementales. Sin embargo, existen algoritmos numéricos muy buenos (¡y razonablemente rápidos!) para calcular tal cantidad y deberían preferirse a cualquier técnica de integración numérica en este caso.

Cantidad de interés en términos de CDF normal

En realidad, la cantidad que te interesa está estrechamente relacionada con la media condicional de una variable aleatoria lognormal. Es decir, si $X$ se distribuye como una lognormal con parámetros $\mu$ y $\sigma$ entonces, usando tu notación, $$ \newcommand{\e}{\mathbb{E}}\renewcommand{\Pr}{\mathbb{P}}\newcommand{\rd}{\mathrm{d}} \int_a^b f(x) \rd x = \int_a^b \frac{1}{\sigma \sqrt{2\pi}} e^{-\frac{1}{2\sigma^2}(\log(x) - \mu)^2} \rd x = \Pr(a \leq X \leq b) \e(X \mid a \leq X \leq b) \>. $$

Para obtener una expresión de esta integral, haz la sustitución $z = (\log(x) - (\mu + \sigma^2))/\sigma$ . Al principio, esto puede parecer un poco inmotivado. Pero, tenga en cuenta que el uso de esta sustitución, $x = e^{\mu + \sigma^2} e^{\sigma z}$ y por un simple cambio de variables, obtenemos $$ \int_a^b f(x) \rd x = e^{\mu + \frac{1}{2}\sigma^2} \int_{\alpha}^{\beta} \frac{1}{\sqrt{2\pi}} e^{-\frac{1}{2} z^2} \rd z \> , $$ donde $\alpha = (\log(a) - (\mu + \sigma^2))/\sigma$ y $\beta = (\log(b) - (\mu + \sigma^2))/\sigma$ .

Por lo tanto, $$ \int_a^b f(x) \rd x = e^{\mu + \frac{1}{2}\sigma^2} \big( \Phi(\beta) - \Phi(\alpha) \big) \>, $$ donde $\Phi(x) = \int_{-\infty}^x \frac{1}{\sqrt{2\pi}} e^{-z^2/2} \rd z$ es la función de distribución acumulativa normal estándar.

Aproximación numérica

A menudo se afirma que no se conoce ninguna expresión de forma cerrada para $\Phi(x)$ existe. Sin embargo, un teorema de Liouville de principios de 1800 afirma algo más fuerte: No existe una expresión cerrada para esta función . (Para la demostración en este caso particular, véase Informe de Brian Conrad .)

Por lo tanto, sólo nos queda utilizar un algoritmo numérico para aproximar la cantidad deseada. Esto puede hacerse dentro de la coma flotante de doble precisión IEEE mediante un algoritmo de W. J. Cody. Se trata de el algoritmo estándar para este problema, y utilizando expresiones racionales de un orden bastante bajo, también es bastante eficiente.

Aquí hay una referencia que discute la aproximación:

W. J. Cody, Aproximaciones racionales de Chebyshev para la función de error , Matemáticas. Comp. , 1969, pp. 631--637.

También es la implementación utilizada tanto en MATLAB como en $R$ entre otros, por si facilitan la obtención de código de ejemplo.

Aquí es una pregunta relacionada, por si le interesa.