La distribución estable positiva en R

Las distribuciones estables positivas se describen mediante cuatro parámetros: el parámetro de asimetría , el parámetro de escala , el parámetro de ubicación y llamado parámetro de índice . Cuando es cero, la distribución es simétrica alrededor de , cuando es positiva (resp. negativa) la distribución está sesgada a la derecha (resp. a la izquierda) Distribuciones estables permiten colas gruesas cuando disminuye.σ > 0 μ ∈ ( - ∞ , ∞ ) α ∈ ( 0 , 2 ] β μ $\beta\in[-1,1]$$\sigma>0$$\mu\in(-\infty,\infty)$$\alpha\in(0,2]$$\beta$$\mu$$\alpha$

Cuando es estrictamente menor que uno y el soporte de la distribución se limita a .β = 1 ( μ , ∞ $\alpha$$\beta=1$$(\mu,\infty)$

La función de densidad solo tiene una expresión de forma cerrada para algunas combinaciones particulares de valores para los parámetros. Cuando , , y es (vea la fórmula (4.4) aquí ):α < 1 β = 1 σ = $\mu=0$$\alpha<1$$\beta=1$$\sigma=\alpha$

$f(y) = -\frac{1}{\pi y} \sum_{k=1}^{\infty} \frac{\Gamma(k\alpha+1)}{k!} (-y^{-\alpha})^k \sin(\alpha k \pi)$

Tiene media y varianza infinitas.

Pregunta

Me gustaría usar esa densidad en R. Yo uso

> alpha <- ...
> dstable(y, alpha=alpha, beta=1, gamma=alpha, delta=0, pm=1)

donde la función dstable viene con el paquete fBasics.

¿Puedes confirmar que esta es la forma correcta de calcular esa densidad en R?

¡Gracias de antemano!

EDITAR

Una razón por la que sospecho es que, en la salida, el valor de delta es diferente al de la entrada. Ejemplo:

> library(fBasics)
> alpha <- 0.4
> dstable(4, alpha=alpha, beta=1, gamma=alpha, delta=0, pm=1)
[1] 0.02700602
attr(,"control")
   dist alpha beta gamma    delta pm
stable   0.4    1   0.4 0.290617  1

La respuesta corta es que su está bien, pero su está mal. Para obtener la distribución estable positiva dada por su fórmula en R, debe establecer γ γ = | 1 - i tan ( π α / 2 ) | - 1 / α$\delta$$\gamma$

$$\gamma = |1 - i \tan \left(\pi \alpha / 2\right)|^{-1/\alpha}.$$

El primer ejemplo que pude encontrar de la fórmula que dio fue en (Feller, 1971), pero solo encontré ese libro en forma física. Sin embargo (Hougaard, 1986) proporciona la misma fórmula, junto con la transformación de Laplace Del manual ( se usa en ), la parametrización es de (Samorodnitsky y Taqqu, 1994), otro recurso cuya reproducción en línea me ha eludido. Sin embargo (Weron, 2001) da la función característica en Samorodnitsky y la parametrización de Taqqu para que sea α ≠ 1 φ ( t ) = E [ exp ( i t X ) ] = exp [ i δ t - γ α | t | α ( 1 - i β s i

$$\mathrm{L}(s) = \mathrm{E}\left[\exp(-sX)\right] = \exp\left(-s^\alpha\right).$$ stablediststabledistfBasicspm=1$\alpha \neq 1$μδσγβ=1δ=0φ(t)=exp[-γα| t| α(1-isign(t)tanπ $$\varphi(t) = \mathrm{E}\left[\exp(i t X) \right] = \exp\left[i \delta t - \gamma^\alpha |t|^\alpha \left(1 - i \beta \mathrm{sign}(t) \tan{\frac{\pi \alpha}{2}} \right) \right].$$ Cambié el nombre de algunos parámetros del papel de Weron a coinide con la notación que estamos usando. Utiliza para y para . En cualquier caso, conectando y , obtenemos $\mu$$\delta$$\sigma$$\gamma$$\beta=1$$\delta=0$ $$\varphi(t) = \exp\left[-\gamma^\alpha |t|^\alpha \left(1 - i \mathrm{sign}(t) \tan \frac{\pi \alpha}{2} \right) \right].$$

Tenga en cuenta que para y que . Formalmente, , entonces estableciendo en obtenemos Un punto interesante a tener en cuenta es que la que corresponde a también es , por lo que si tuviera que probar o$(1 - i \tan (\pi \alpha / 2)) / |1 - i \tan(\pi \alpha / 2)| = \exp(-i \pi \alpha / 2)$$\alpha \in (0, 1)$$i^\alpha = \exp(i \pi \alpha / 2)$γ = | 1 - i tan ( π α / 2 ) | - 1 / α φ ( t ) φ ( i s ) = exp ( - s α ) = L ( s ) . γ alpha = 1 / 2 1 / 2 γ = alpha γ = 1 - alpha $\mathrm{L}(s)=\varphi(is)$$\gamma = |1 - i \tan \left(\pi \alpha / 2\right)|^{-1/\alpha}$$\varphi(t)$

$$\varphi(is) = \exp\left(-s^\alpha\right) = \mathrm{L}(s).$$ $\gamma$$\alpha=1/2$$1/2$$\gamma=\alpha$$\gamma=1-\alpha$, que en realidad no es una mala aproximación, terminas exactamente correcto para .$\alpha=1/2$

Aquí hay un ejemplo en R para verificar la corrección:

library(stabledist)

# Series representation of the density
PSf <- function(x, alpha, K) {
  k <- 1:K
  return(
    -1 / (pi * x) * sum(
      gamma(k * alpha + 1) / factorial(k) * 
        (-x ^ (-alpha)) ^ k * sin(alpha * k * pi)
    )
  )
}

# Derived expression for gamma
g <- function(a) {
  iu <- complex(real=0, imaginary=1)
  return(abs(1 - iu * tan(pi * a / 2)) ^ (-1 / a))
}

x=(1:100)/100
plot(0, xlim=c(0, 1), ylim=c(0, 2), pch='', 
     xlab='x', ylab='f(x)', main="Density Comparison")
legend('topright', legend=c('Series', 'gamma=g(alpha)'),
       lty=c(1, 2), col=c('gray', 'black'),
       lwd=c(5, 2))
text(x=c(0.1, 0.25, 0.7), y=c(1.4, 1.1, 0.7), 
     labels=c(expression(paste(alpha, " = 0.4")),
              expression(paste(alpha, " = 0.5")),
              expression(paste(alpha, " = 0.6"))))

for(a in seq(0.4, 0.6, by=0.1)) {
  y <- vapply(x, PSf, FUN.VALUE=1, alpha=a, K=100)
  lines(x, y, col="gray", lwd=5, lty=1)
  lines(x, dstable(x, alpha=a, beta=1, gamma=g(a), delta=0, pm=1), 
        col="black", lwd=2, lty=2)
}

$\hskip1in$

1. Feller, W. (1971). Una introducción a la teoría de la probabilidad y sus aplicaciones , 2 , 2ª ed. Nueva York: Wiley.
2. Hougaard, P. (1986). Modelos de supervivencia para poblaciones heterogéneas derivadas de distribuciones estables , Biometrika 73 , 387-396.
3. Samorodnitsky, G., Taqqu, MS (1994). Procesos aleatorios no gaussianos estables , Chapman & Hall, Nueva York, 1994.
4. Weron, R. (2001). Distribuciones estables a Levy revisitadas: índice de cola> 2 no excluye el régimen estable a Levy , International Journal of Modern Physics C, 2001, 12 (2), 209-223.

Lo que creo que está sucediendo es que en la salida deltase puede informar un valor de ubicación interna, mientras que en la entrada deltase describe el cambio. [Parece haber un problema similar con gammacuándo pm=2.] Entonces, si intentas aumentar el cambio a 2

> dstable(4, alpha=0.4, beta=1, gamma=0.4, delta=2, pm=1)
[1] 0.06569375
attr(,"control")
   dist alpha beta gamma    delta pm
 stable   0.4    1   0.4 2.290617  1

luego agrega 2 al valor de ubicación.

Con beta=1y pm=1tiene una variable aleatoria positiva con una distribución de límite inferior en 0.

> min(rstable(100000, alpha=0.4, beta=1, gamma=0.4, delta=0, pm=1))
[1] 0.002666507

Cambia en 2 y el límite inferior aumenta en la misma cantidad

> min(rstable(100000, alpha=0.4, beta=1, gamma=0.4, delta=2, pm=1))
[1] 2.003286

Pero si desea que la deltaentrada sea el valor de ubicación interna en lugar del desplazamiento o límite inferior, debe usar una especificación diferente para los parámetros. Por ejemplo, si intenta lo siguiente (con pm=3y probando delta=0y lo delta=0.290617que encontró anteriormente), parece que obtiene lo mismo deltadentro y fuera. Con pm=3y delta=0.290617obtienes la misma densidad de 0.02700602 que encontraste anteriormente y un límite inferior en 0. Con pm=3y delta=0obtienes un límite inferior negativo (de hecho -0.290617).

> dstable(4, alpha=0.4, beta=1, gamma=0.4, delta=0, pm=3)
[1] 0.02464434
attr(,"control")
   dist alpha beta gamma delta pm
 stable   0.4    1   0.4     0  3
> dstable(4, alpha=0.4, beta=1, gamma=0.4, delta=0.290617, pm=3)
[1] 0.02700602
attr(,"control")
   dist alpha beta gamma    delta pm
 stable   0.4    1   0.4 0.290617  3
> min(rstable(100000, alpha=0.4, beta=1, gamma=0.4, delta=0, pm=3))
[1] -0.2876658
> min(rstable(100000, alpha=0.4, beta=1, gamma=0.4, delta=0.290617, pm=3))
[1] 0.004303485

Puede que le resulte más fácil simplemente ignorar deltaen la salida, y siempre y cuando continúe beta=1usando pm=1medios deltaen la entrada es el límite inferior de distribución, que parece que quiere ser 0.

También es de destacar: Martin Maechler simplemente refactorizó el código para la distribución estable y agregó algunas mejoras.

Su nuevo paquete stabledist también será utilizado por fBasics, por lo que es posible que desee echarle un vistazo también.