正如评论中所建议的，您可以使用fitdistr, 和 中的密度函数fBasics。

# Sample data
x <- rt(100,df=4)

# Density (I reparametrize it to remove the constraints 
# on the parameters)
library(fBasics)
library(stabledist)
f <- function(u,a,b,c,d) {
  cat(a,b,c,d,"\n")  # Some logging (it is very slow)
  dstable(u, 2*exp(a)/(1+exp(a)), 2*exp(b)/(1+exp(b))-1, exp(c), d)
}

# Fit the distribution
library(MASS)
r <- fitdistr(x, f, list(a=1, b=0, c=log(mad(x)), d=median(x)))
r

# Graphical check
plot(
  qstable(ppoints(100),
    2*exp(r$estimate[1])/(1+exp(r$estimate[1])),
    2*exp(r$estimate[2])/(1+exp(r$estimate[2]))-1,
    exp(r$estimate[3]),
    r$estimate[4]
  ),
  sort(x)
)
abline(0,1)

拟合参数的一种方法是通过 Okoneshnikov 描述的 Nagaev变换。另一种方法是Mantegna 和 Stanley的“返回概率”方法，该方法要容易得多。$\alpha$

编辑：另一种“经典”方法是 Kogon & Williams 的（SM Kogon，Douglas B. Williams，“基于特征函数的稳定分布参数估计技术”），另见K&W 的 matlab 实现

@Vincent 的回答听起来不错，但这是另一种方法：既然您知道分布的傅立叶变换，请对数据进行适当的傅立叶变换，并找到最适合傅立叶空间的参数。

我认为这种方法在理论上应该同样有效，在实践中会避免大量的数值积分来获得稳定分布的形式。我现在不编写测试代码，抱歉。有人对此有任何见解吗？