5

特性関数が既知の分布の密度関数を計算したいと考えています。簡単な例として、正規分布を取り上げます。

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の混合です。誰でも手順を推奨できますか？ありがとう

4

1 に答える 1

11

ただ面倒です: ペンと紙を取り、計算したい積分 (特性関数のフーリエ変換) を書き、それを離散化し、離散フーリエ変換のように見えるように項を書き直します (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)
于 2012-04-06T01:27:41.317 に答える