[R] convolve bug?

[Martin Maechler]

>>>>> On Fri, 19 Nov 1999 15:09:17 +0000, yudi at ucd.ie (Yudi Pawitan) said:

    Yudi> I'm not sure if the group gets this.
    >> Date: Fri, 19 Nov 1999 10:09:17 +0000 (GMT) From: Bill Simpson
    >> <wsi at gcal.ac.uk> BTW I don't see why you say it is OK to define
    >> convolution with no reversal of the function!  If you don't reverse,
    >> the operation is CORRELATION. Or do statisticians make no
    >> distinction between convolution and correlation? Anyway I like the
    >> definition of convolution given on the help page. (Even though it is
    >> not what is happening)

    >> This is really sad that the engineers and stats types have such
    >> totally different conventions.

    >> To me, the linear system output y(t) is the convolution of the input
    >> x(t) and the impulse response h(t): y(t) = x(t) * h(t) where * is
    >> convolution.  But this means something completely different to a
    >> statistician, who in my language would interpret this as "linear
    >> system output y(t) is the correlation of the input x(t) and the
    >> impulse response h(t)", which is wrong.
    >> 
    >> The convolution and correlation theorems are also screwed up (the
    >> complex conjugation).
    >> 
    >> I'm glad I know about this sad state of affairs now! The help page
    >> refs Brillinger, a stats type (though I know he has published papers
    >> in IEEE journals, so he must constantly have to be flipping things
    >> in his head!).  I will look at his book.

    Yudi> Bill, I can assure you that statisticians use 'convolution' in
    Yudi> the same way as engineers do. It is actually the defaults in R
    Yudi> that are rather nonstandard (and the help is wrong).

yes.  I have changed the help already and was trying to add more useful
 explanation (and I was the one who added the wrong help part for "open",...).

    Yudi> I agree with Bill that convolve() should mean convolution
    Yudi> operation. So the default 'conj' should be set to F (Fourier
    Yudi> transform of a convolution is a product of the transfroms), and
    Yudi> the zero padding is on the right.

I tend to agree....  see below

    Yudi> E.g., target:  (2,2,3,3,4)*(1,1,2)  =  (2  4  9 10 13 10  8 )

    Yudi> convolve(c(2,2,3,3,4),c(2,1,1),type='o')
    Yudi> #[1]  2  4  9 10 13 10  8  
    Yudi> convolve(c(2,2,3,3,4),c(1,1,2),conj=F,type='o')
    Yudi> #[1] 10  8  2  4  9 10 13      # almost 
    Yudi> x _ c(2,2,3,3,4,0,0)   # zero padding
    Yudi> y _ c(1,1,2,0,0,0,0)
    Yudi> aa_ fft(fft(x) * fft(y), inv = TRUE)
    Yudi> Re(aa)/length(aa)
    Yudi> #[1]  2  4  9 10 13 10  8   # wanted
(there are indeed several ways to get this result,  using
 combinations of  Conj(), rev() and/or  0- padding...      

    Yudi> I would suggest that the 'conj' option not be provided; when
    Yudi> conj=F, the type='f' option gives meaningless answer:
?? Contradiction:
Above you say that  conj=F should be the *default* ...

    Yudi> convolve(c(2,2,3,3,4),c(1,1,2),conj=F,type='f') #[1] 2 4 9

I agree that at least    convolve(x,y) 
should really give what everyone expects..

One way to achieve this would be a new type, 
e.g. "regular" , "default" or "polynomial".

And we would change the *default* 
   -   from type = "circular" to  type = "regular" (or whatever we call it).
   -   from conj = TRUE  to   conj = type != "circular"
	    (conj defaulting to FALSE when type="open", ..)

or do you think conj should default to FALSE even for circular?
(I must admit not having seen a definition of a circular convolution in a
 text book..)

Here is the longer story:
  "convolve()" was part of R even before there were version numbers
  (starting with 0.1)..., only with the circular type.
  It is has been used inside  density() since "ever", and
  its use inside density *is* the one with the current defaults.

  Then, "open"/"filter" were added in June, where I was following a
  proposition from Martyn Plummer.

--
Yes, I'm waiting for further propositions, and I agree the defaults should
be changed.

Unfortunately, it's too late for the imminent 0.90.0 release (Monday'ish).

Martin Maechler <maechler at stat.math.ethz.ch>	http://stat.ethz.ch/~maechler/
Seminar fuer Statistik, ETH-Zentrum  LEO D10	Leonhardstr. 27
ETH (Federal Inst. Technology)	8092 Zurich	SWITZERLAND
phone: x-41-1-632-3408		fax: ...-1228			<><
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._