在R中使用fft从特征函数计算密度 [英] Calculating a density from the characteristic function using fft in R

问题描述

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

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 like a discrete Fourier transform (the FFT assumes that the interval starts at 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从特征函数计算密度的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！