Hello, May be somebody can help me... I am trying to find a solution of a convolution equation using fft (and unfortunately I do not have a good background for this). So I am just trying to figure out how it can be implemented in R. I have two multidimensional independent variables X and Z and I know their densities fx and fz, which are multidimensional arrays. So I have to find the density of Y, such that X+Y=Z. So, first I tried on a simple example, where the variables are one-dimensional, say X is N(0,1) and Z is N(7,3). So I want to find the density of Y (which should be N(7, sqrt(8)). I did: x <- seq(-6, 20, length=500) fx <- dnorm(x) z <- seq(-6, 20, length=500) fz <- dnorm(z, mean=7, sd=3) ffty <- fft(fz)/fft(fx) fy <- fft(ffty, inverse=T)/length(ffty) plot(Re(fy)) But the plot is far from being normal. May be the problem is with the lengths of fx and fz? should they be of different lengths and fx padded with zeros? May be somebody could point out that?? Thanks! Anna
You need to pad both fx and fz with zeros to at least a length of length(fx) + length(fz) - 1, because you're really computing a circular convolution. The same holds in higher dimensions (for each dimension). Reid Huntsinger -----Original Message----- From: r-help-bounces at stat.math.ethz.ch [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Anna Oganyan Sent: Thursday, September 29, 2005 3:23 PM To: r-help at stat.math.ethz.ch Subject: [R] solution of convolution equation Hello, May be somebody can help me... I am trying to find a solution of a convolution equation using fft (and unfortunately I do not have a good background for this). So I am just trying to figure out how it can be implemented in R. I have two multidimensional independent variables X and Z and I know their densities fx and fz, which are multidimensional arrays. So I have to find the density of Y, such that X+Y=Z. So, first I tried on a simple example, where the variables are one-dimensional, say X is N(0,1) and Z is N(7,3). So I want to find the density of Y (which should be N(7, sqrt(8)). I did: x <- seq(-6, 20, length=500) fx <- dnorm(x) z <- seq(-6, 20, length=500) fz <- dnorm(z, mean=7, sd=3) ffty <- fft(fz)/fft(fx) fy <- fft(ffty, inverse=T)/length(ffty) plot(Re(fy)) But the plot is far from being normal. May be the problem is with the lengths of fx and fz? should they be of different lengths and fx padded with zeros? May be somebody could point out that...? Thanks! Anna ______________________________________________ R-help at stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Anna Oganyan wrote:>Hello, >May be somebody can help me... >I am trying to find a solution of a convolution equation using fft (and >unfortunately I do not have a good background for this). >So I am just trying to figure out how it can be implemented in R. I have >two multidimensional independent variables X and Z >and I know their densities fx and fz, which are multidimensional arrays. >So I have to find the density of Y, such that X+Y=Z. >So, first I tried on a simple example, where the variables are >one-dimensional, say X is N(0,1) and Z is N(7,3). >So I want to find the density of Y (which should be N(7, sqrt(8)). >I did: >x <- seq(-6, 20, length=500) >fx <- dnorm(x) >z <- seq(-6, 20, length=500) >fz <- dnorm(z, mean=7, sd=3) >ffty <- fft(fz)/fft(fx) >fy <- fft(ffty, inverse=T)/length(ffty) >plot(Re(fy)) > >But the plot is far from being normal. May be the problem is with the >lengths of fx and fz? should they be of different lengths and fx padded >with zeros? May be somebody could point out that?? Thanks! >Anna > >______________________________________________ >R-help at stat.math.ethz.ch mailing list >https://stat.ethz.ch/mailman/listinfo/r-help >PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html > >Hello Anna, in our R package "distr" (on CRAN) we have implemented a convolution algorithm via fft. See also: http://www.uni-bayreuth.de/departments/math/org/mathe7/KOHL/pubs/comp.pdf respectively library(distr) getMethod("+", signature=c("AbscontDistribution","AbscontDistribution")) (or see the sources) Unfortunatelly, we haven't implemented our algorithms for multidimensional distributions so far. hope that helps, Matthias