La réponse courte est que votre va bien, mais votre est faux. Pour obtenir la distribution stable positive donnée par votre formule dans R, vous devez définir
γ γ = | 1 - i tan ( π α / 2 ) | - 1 / α .δγ
γ= | 1 - je bronze( πα / 2 ) |- 1 / α.
Le premier exemple que j'ai pu trouver de la formule que vous avez donnée était (Feller, 1971), mais je n'ai trouvé ce livre que sous forme physique. Cependant (Hougaard, 1986) donne la même formule, avec la transformée de Laplace
À partir du manuel ( utilisé dans ), le paramétrage provient de (Samorodnitsky et Taqqu, 1994), une autre ressource dont la reproduction en ligne m'a échappé. Cependant (Weron, 2001) donne la fonction caractéristique dans le paramétrage de Samorodnitsky et Taqqu pour que soit
L (s)= E [ exp( - s X) ] = exp( - sα) .
stabledist
stabledist
fBasics
pm=1
α ≠ 1φ ( t ) = E [ exp( i t X) ] = exp[ i δt - γα| t |α( 1 - i βs i g n (t)tanπα2) ] .
J'ai renommé certains paramètres du papier de Weron en coin avec la notation que nous utilisons. Il utilise pour et pour . Dans tous les cas, en branchant et , nous obtenons
μδσγβ= 1δ= 0φ ( t ) = exp[ - γα| t |α( 1 - i s i g n ( t ) tanπα2) ] .
Notez que pour et que . Formellement, , donc en définissant dans nous obtenons
Un point intéressant à noter est que le qui correspond à est également , donc si vous essayiez ou( 1 - je bronzais( πα / 2 ) ) / | 1 - je bronze( πα / 2 ) | = exp( - i πα / 2 )α ∈ ( 0 , 1 )jeα= exp( i πα / 2 )γ = | 1 - i tan ( π α / 2 ) | - 1 / α φ ( t ) φ ( i s ) = exp ( - s α ) = L ( s ) . γ alpha = 1 / 2 1 / 2 γ = alpha γ = 1 - alpha alphaL (s)=φ(is)γ= | 1 - je bronze( πα / 2 ) |- 1 / αφ ( t )
φ ( i s ) = exp( - sα) = L ( s ) .
γα = une / deux1 / 2γ=αγ=1−α, ce qui n'est en fait pas une mauvaise approximation, vous vous retrouvez exactement correct pour .
α=1/2
Voici un exemple dans R pour vérifier l'exactitude:
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)
}
- Feller, W. (1971). An Introduction to Probability Theory and Its Applications , 2 , 2e éd. New York: Wiley.
- Hougaard, P. (1986). Modèles de survie pour les populations hétérogènes dérivées de distributions stables , Biometrika 73 , 387-396.
- Samorodnitsky, G., Taqqu, MS (1994). Stable Non-Gaussian Random Processes , Chapman & Hall, New York, 1994.
- Weron, R. (2001). Revues stables de Levy revisitées: indice de queue> 2 n'exclut pas le régime stable de Levy , International Journal of Modern Physics C, 2001, 12 (2), 209-223.
Je pense que ce qui se passe est que dans la sortie, il
delta
peut y avoir une valeur d'emplacement interne, tandis que dans l'entréedelta
décrit le changement. [Il semble y avoir un problème similaire avecgamma
quandpm=2
.] Donc, si vous essayez d'augmenter le décalage à 2puis vous ajoutez 2 à la valeur d'emplacement.
Avec
beta=1
etpm=1
vous avez une variable aléatoire positive avec une distribution borne inférieure à 0.Décaler de 2 et la borne inférieure augmente du même montant
Mais si vous souhaitez que l'
delta
entrée soit la valeur d'emplacement interne plutôt que le décalage ou la limite inférieure, vous devez utiliser une spécification différente pour les paramètres. Par exemple, si vous essayez ce qui suit (avecpm=3
et en essayantdelta=0
et quedelta=0.290617
vous avez trouvé plus tôt), vous semblez obtenir la même chosedelta
dans et hors. Avecpm=3
etdelta=0.290617
vous obtenez la même densité de 0,02700602 que vous avez trouvée plus tôt et une borne inférieure à 0. Avecpm=3
etdelta=0
vous obtenez une borne inférieure négative (en fait -0,290617).Vous pouvez trouver plus facile d'ignorer simplement
delta
dans la sortie, et tant que vous continuez àbeta=1
utiliser despm=1
moyensdelta
dans l'entrée, la distribution est la borne inférieure, qui semble vouloir être 0.la source
A noter également: Martin Maechler vient de refactoriser le code de l'écurie distribuée et a ajouté quelques améliorations.
Son nouveau package stabledist sera également utilisé par fBasics, vous pouvez donc également y jeter un œil.
la source