在R中使用fft从特征函数计算密度

问题描述

我想计算其特征函数已知的分布的密度函数.作为一个简单的例子，采用正态分布.

I would like to calculate a density function of a distribution whose characteristics function is known. As a simple example take the normal distribution.

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

，然后我想使用R的fft函数.但是我没有正确地获得乘法常数，因此我不得不重新排列结果的顺序(取值的第二个一半，然后是第一个一半).我尝试过

and then I would like to use R's fft function. but I don't get the multiplicative constants right and I have to reorder the result (take the 2nd half and then the first half of the values). I tried something like

 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的混合.有人可以推荐一个程序吗?谢谢

But this is still not correct. Does anybody know a good reference on the fft in R in the context of density calculations? Obviously the problem is the mixture of the continuous FFT and the discrete one. Can anybody recommend a procedure?Thanks

推荐答案

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

It is just cumbersome: take a pen and paper,write the integral you want to compute(the Fourier transform of the characteristic function),discretize it, and rewrite the terms so that they look likea discrete Fourier transform (the FFT assumes that the interval startsat zero).

请注意，fft是未归一化的转换:没有1/N因子.

Note that fft is an unnormalized transform: there is no 1/N factor.

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)

这篇关于在R中使用fft从特征函数计算密度的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！