有人询问函数的卷积怎么计算,这是一个问题。
两个函数卷积,把积分离散求和,硬算也能算出来。
但是N个呢?复杂度太高了!所以用Fourier transformation可以大大简化计算。
对应统计的问题是,已知
<bblatex>X_1,X_2,...iid \sim X \sim F,</bblatex>计算<bblatex>S_n = X_1+...+X_n</bblatex>的分布函数----即F的n重卷积?对于特殊的分布,可能能够得到显示解。但是,一般的分布只能依赖数值解。
用特征函数的知识,先将问题简化---请注意特征函数本质是Fourier transformation。
<bblatex>\phi(S_n,t) = \phi(X,t)^n</bblatex>,
再由反转公式:<bblatex>f(S_n,y) = \frac{1}{2\pi} \int exp\{-iyt\}\phi(S_n,t)dt</bblatex>即可得到<bblatex>S_n</bblatex>之密度(Durrett书上有,这本书刚有人贴出了pdf版)。
当然上面的等式成立需要一些正则化条件,实际问题我们先不考虑数学的严谨性。
下面以正态分布为例,因为正态时候我们有<bblatex>S_n</bblatex>的显示表达,因此可以验证算法是否正确。
<br />
##convolution<br />
##实部需要的函数<br />
psi1 = function(x,t)<br />
{<br />
dnorm(x)*cos(t*x)<br />
}<br />
##虚部需要的函数<br />
psi2 = function(x,t)<br />
{<br />
dnorm(x)*sin(t*x)<br />
}<br />
Phi1 = 1:201<br />
Phi2 = 1:201<br />
T = -10 + (Phi1 - 1)*0.1<br />
##characteristic function<br />
for(i in 1:201)<br />
{<br />
Phi1[i] = integrate(psi1, t = -10+(i-1)*0.1, lower = -Inf, upper = Inf)$value<br />
Phi2[i] = integrate(psi2, t = -10+(i-1)*0.1, lower = -Inf, upper = Inf)$value<br />
}<br />
##plot<br />
pdf("convolution based on Chf 100.pdf")<br />
par(mfrow=c(2,2))<br />
plot(x = seq(-10,10,0.1), y= Phi1,type="l")<br />
title("Real part of characteristic")<br />
plot(x = seq(-10,10,0.1), y= Phi2,type="l")<br />
title("Imaginary part of characteristic")<br />
##从图可以看出,这步很容易算准</p>
<p>##power n = 100: 100次卷积<br />
n = 100<br />
Phi = complex(real = Phi1, imaginary = Phi2)<br />
nPhi = Phi^n<br />
nPhi1 = Re(nPhi)<br />
nPhi2 = Im(nPhi)</p>
<p>##density:用反转公式,数值积分<br />
Y = T<br />
M = matrix(1:(201*201),ncol=201)<br />
for(i in 1:201)<br />
{<br />
M[i,] = cos(T*Y[i])*nPhi1 + sin(T*Y[i])*nPhi2<br />
}</p>
<p>den = apply(M,1,sum)*0.1/(2*pi)<br />
plot(x = T, y = den, type="l")<br />
lines(x = T, y = dnorm(T,0,sqrt(n))+0.0005, col = "red")<br />
##画高0.005是为了看出有两条线,否则重合<br />
title("True(red) vs Estimate(black)")</p>
<p>plot(dnorm(T,0,sqrt(n))~den, xlim=c(min(den),max(den)),ylim=c(min(den),max(den)))<br />
title("True vs Estimated")<br />
dev.off()</p>
<p>
</p>
上面其实每步都很粗糙,有进一步改进的余地,这里只是提供一种方法而已。