【问题标题】:Calculating a density from the characteristic function using fft in R使用 R 中的 fft 从特征函数计算密度
【发布时间】:2014-12-23 04:13:26
【问题描述】:

我想计算特征函数已知的分布的密度函数。以正态分布为例。

norm.char<-function(t,mu,sigma) exp((0+1i)*t*mu-0.5*sigma^2*t^2)

然后我想使用 R 的 fft 函数。但我没有得到正确的乘法常数,我必须重新排序结果(取第二半,然后取值的前半)。我尝试了类似

 xmax = 5
 xmin = -5
 deltat = 2*pi/(xmax-xmin)
 N=2^8
 deltax = (xmax-xmin)/(N-1)
 x = xmin + deltax*seq(0,N-1)
 t = deltat*seq(0,N-1)
 density = Re(fft(norm.char(t*2*pi,mu,sigma)))
 density = c(density[(N/2+1):N],density[1:(N/2)])

但这仍然不正确。在密度计算的背景下,有人知道 R 中 fft 的一个很好的参考吗?显然,问题在于连续 FFT 和离散 FFT 的混合。有人可以推荐一个程序吗? 谢谢

【问题讨论】:

  • density 函数帮助页面说它使用 FFT。为什么不审查代码?
  • 究竟有什么不正确的?如果您真正的问题只是“在离散傅里叶变换期间应用了哪些常数?”然后查看fft 的帮助页面,我相信它给出了方程式。

标签: r fft


【解决方案1】:

这很麻烦:拿笔和纸, 写出你要计算的积分 (特征函数的傅里叶变换), 离散化它,并重写这些术语,使它们看起来像 离散傅立叶变换(FFT 假设区间开始 为零)。

请注意,fft 是非标准化变换:没有 1/N 因子。

characteristic_function_to_density <- function(
  phi, # characteristic function; should be vectorized
  n,   # Number of points, ideally a power of 2
  a, b # Evaluate the density on [a,b[
) {
  i <- 0:(n-1)            # Indices
  dx <- (b-a)/n           # Step size, for the density
  x <- a + i * dx         # Grid, for the density
  dt <- 2*pi / ( n * dx ) # Step size, frequency space
  c <- -n/2 * dt          # Evaluate the characteristic function on [c,d]
  d <-  n/2 * dt          # (center the interval on zero)
  t <- c + i * dt         # Grid, frequency space
  phi_t <- phi(t)
  X <- exp( -(0+1i) * i * dt * a ) * phi_t
  Y <- fft(X)
  density <- dt / (2*pi) * exp( - (0+1i) * c * x ) * Y
  data.frame(
    i = i,
    t = t,
    characteristic_function = phi_t,
    x = x,
    density = Re(density)
  )
}

d <- characteristic_function_to_density(
  function(t,mu=1,sigma=.5) 
    exp( (0+1i)*t*mu - sigma^2/2*t^2 ),
  2^8,
  -3, 3
)
plot(d$x, d$density, las=1)
curve(dnorm(x,1,.5), add=TRUE)

【讨论】:

    猜你喜欢
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2017-09-20
    • 1970-01-01
    • 2018-03-17
    • 2015-06-28
    • 2020-05-17
    • 1970-01-01
    相关资源
    最近更新 更多