Google Groups no longer supports new Usenet posts or subscriptions. Historical content remains viewable.
Dismiss

Mathematicians-- need help with convolution macro

4,087 views
Skip to first unread message

Curt Dye

unread,
Dec 2, 2001, 9:51:15 PM12/2/01
to
I am trying write a macro to convolute the values in 2 columns. For those
that don't know, given

1 4
2 5
3 6

the convoluted values would be

1*4 = 4
2*4 + 1*5 = 13
3*4 + 2*5 + 1*6 = 28
3*5 + 2*6 + 1*0 = 27
3*6 + 2*0 + 1*0 = 18

If there are a small number of values in the column, I can handle it, but I
would like to be able to handle data sets larger than 10. My main problem
is going from the formula to the values... specifically for the values
greater than the order of the smallest column.

Here is what I have so far... Anyone know of an easier method or how to
accomplish the last part for values? Any help is appreciated :
'----------------------------

Function ConvCol(C1 As Integer, C2 As Integer)
'Puts in convolution formula/value for columns C1 & C2 in ActiveCell.Column
'Dims : omitted

R = ActiveCell.Row
C = ActiveCell.Column
NumC1 = Cells(Rows.Count, C1).End(xlUp).Row - R + 1
NumC2 = Cells(Rows.Count, C2).End(xlUp).Row - R + 1
'Number of values in Column C1 and C2 : needs improvement

For RwIndex = 1 To NumC1 'It is assumed that C1 has fewer values
Frm = "="
V = 0
For J = 1 To RwIndex
Frm = Frm & "+" & ColNm(C1) & "$" & RwIndex + R - J & "*" _
& ColNm(C2) & R + J - 1 'ColNm gives the letter
corresponding to the column #
V1 = ActiveSheet.Cells(RwIndex + R - J, C1).Value
V2 = ActiveSheet.Cells(R + J - 1, C2).Value
V = V + (V1 * V2)
Next J
Cells(RwIndex + R - 1, C).Formula = Frm
Cells(RwIndex + R - 1, C).Offset(0, 1).Value = V ' just to make sure
I'm getting the right values
Next RwIndex

'This part works for the formula, but not for the values...
Cells(NumC1 + R - 1, C).Copy
myRange = Range(Cells(NumC1 + R - 1, C), Cells(NumC1 + NumC2 + R - 2,
C)).Address
ActiveSheet.Paste Destination:=ActiveSheet.Range(myRange)
End Function

' ------------------ Used above
Function ColNm(ColNum As Integer) As String
'Input a given a column number, outputs the corresponding letter
Dim S As String
If ColNum <= 26 Then
S = Chr(64 + ColNum)
Else
S = Chr(64 + Int(ColNum / 26)) & Chr(64 + 1 + (ColNum - 1) Mod 26)
End If
ColNm = S
End Function

Michael Ernstoff

unread,
Dec 3, 2001, 3:44:47 AM12/3/01
to
The Convolution theorems that I've studied have always dealt with continuous
functions rather discrete values, and they've been of the form.

f(t) = I [p(x) * q(x - t) dx where I represents the integral sign.

Your computations do not seem to fit this formulation even if the integration
is replaced by a summation. Are you sure the term "convolute" is the most
appropriate term for what you are doing? My unabridged dictionary did not have
a mathermatical definition for convolute... Do you mean convolve or
convolution?

What is the application for the process? Certainly others have had the same
problem in the past, and viewing their solutions might lead to greater clarity.

Finally, I suspect that some sort of matrix algebra operations may the best way
to expedite the computations. Consider that the product of a row and column
matrix is a single value equal to the sum of the products. The determinant of
a 2 x 2 matrix is the sum of the products of the cross members. Excel has
canned functions for handling matrices and determinants.

buczacz in Los Angeles

Curt Dye wrote:

--
Michael
Los Angeles, CA

Curt Dye

unread,
Dec 3, 2001, 10:28:58 PM12/3/01
to
I attempting to use this in the calculation of earnings. I'm not sure if
'convolution' is the term I want, but that's what I've always heard it
referred to.

Let's try another simplified example: say we take in a reserve of $1,000
every year for the next four years. Of this $1,000 we expect losses of 70%
in the 1st year of the contracts and 10% in each of the next 3 years. Our
losses per year would thus be:

Year c1 c2 (c1 * c2) = c^(2) = "c1 convoluted
with c2"
1 1000 .7 700
2 1000 .1 Losses of contracts in 1st year +
Losses of contracts in 2nd year = 1000 * .7 + 1000 * .1 = 800
3 1000 .1 1000 * .7 + 1000 * .1 + 1000 * .1 =
900
4 1000 .1 1000 * .7 + 1000 * .1 + 1000 * .1 +
1000 *.1 = 1000
5 1000 * .1 + 1000 * .1 + 1000 *.1
= 300
6 1000 * .1 + 1000 * .1 = 200
7 1000 * .1 = 100
8 0

Notice, the sum of (c1 * c2) = 4000 = sum of c1, since the sum of c2 = 1.

When dealing with probability densities for c1 and c2, I believe you get the
form you show. For example

f^(2)
F1 * f2
f1 f2 F1 F2 f1*f2 F^(2) F^(2)
0.6 0.7 0.6000 0.7000 0.4200 0.4200 0.4200
0.2 0.1 0.8000 0.8000 0.2000 0.6200 0.6200
0.1 0.1 0.9000 0.9000 0.1500 0.7700 0.7700
0.1 0.1 1.0000 1.0000 0.1600 0.9300 0.9300

I don't think the ends taper off the same though (?)


"Michael Ernstoff" <buc...@yahoo.com> wrote in message
news:3C0B39EC...@yahoo.com...

Nico Sterk

unread,
Dec 4, 2001, 5:24:32 AM12/4/01
to
Curt,

You are not really clear about the meaning of the table. Further more a
notation like f^(2) is ambiguous. Please tell me the meaning of f1, f2, F1,
F2, f1*f2 and F^(2). Where can I read probability densities for c1 and c2???

Regards, Nico Sterk.

S...@internet.org

unread,
Dec 4, 2001, 10:13:49 AM12/4/01
to
Curt,

m years of reserves and
k years to payout
results in
(m + k - 1) payout years

A is the payout ratio matrix of dimension m x (m + k - 1)
B is the reserve matrix (column vector) of dimension m x 1

If the result you desire is the total dollars paid out by payout year,
then the matrix operation
Transpose(B) * A
results in a 1 x (m + k - 1) total dollars by payout year row vector.


--S...@internet.org

"Curt Dye" <cd...@tampabay.rr.com> wrote in message news:DMBO7.181848$zK1.49...@typhoon.tampabay.rr.com...

Curt Dye

unread,
Dec 4, 2001, 8:22:42 PM12/4/01
to
Sorry, the formatting didn't send as I viewed it.
In any case, given 2 independent random variables X1 and X2 whose
distributions defined by columns f1(x) and f2(x),
I believe one can derive the p.f. and d.f. of S = X1 + X2 with convolution.
[In terms of probability, think of f1(x) as the discrete probability
function Pr(X1=x)]

f1(0) = .6
f1(1) = .2
f1(2) = .1
f1(3) = .1

f2(0) = .7
f2(1) = .1
f2(2) = .1
f2(3) = .1

F1 is the density function of f1
F2 is the density funtion of f2

F1(0) = f1(0) = .6
F1(1) = f1(0) + f1(1) = .8
F1(2) = .9
F1(3) = 1

Similarly, the values for F2(x) are .7, .8, .9, and 1

f1*f2 is the convolution of f1 and f2. I denoted it f^(2) [as I have seen
elsewhere].
Following this notation, (f1*f2)*f3 is f^(3). In this example, we have

f1*f2(0) = f1(0) * f2(0) = .6 * .7 = .42
f1*f2(1) = f1(0) * f2(1) + f2(0) * f1(1) = .2
f1*f2(2) = .15
f1*f2(3) = .16

Finally, then, F^(2) is the density function of S = X1 + X2 ( or more
generally, F^(n) is the d.f. of X1 + X2 + .... + Xn)
and this can be calculated as:

F^(2) (0) = f1*f2(0) = .42
F^(2) (1) = f1*f2(0) + f1*f2(1) =.42 +.2 = .62
F^(2) (2) = f1*f2(0) + f1*f2(1) + f1*f2(2) = ..42+ .2 + .15 = .77
F^(2) (3) = .93

or by the relation F^(2) = F1 * f2
F^(2) (0) = F1*f2(0) = F1(0)* f2(0) = .6 * .7 = .42
F^(2) (1) = F1*f2(1) = F1(0) * f2(1) + F1(1) * f2(0) = .6*.1 + .8*.7 = .06 +
.56 = .62
... etc

Hope that helps....


"Nico Sterk" <n.s...@chello.nl> wrote in message
news:Av1P7.29841$Iw1.1...@nlnews00.chello.com...

Dana DeLouis

unread,
Dec 6, 2001, 11:34:07 AM12/6/01
to
I have been trying to understand this subject for quite some time. I am
currently reading a few books on the subject, but it is still over my head.
From what I understand so far, your description seems to agree with what I
have read.
Are you trying to improve your VBA routine to do Convolution?

This is a little off subject, but I thought I would try to follow along your
first example using a program called Mathematica.
Here are the two lists of numbers from your first example.

v1 = {3, 2, 1};
v2 = {4, 5, 6};

I reversed your first list from 1, 2, 3 to 3, 2, 1. to match a
"shifting" type of algorithm.
The following command does the Convolution, but with the options to pad both
ends with 0's.

ListCorrelate[v1, v2, {-1, 1}, 0]

{4,13,28,27,18}
The results agree with your answers.


I am guessing that your abbreviation of d.f. is the Discrete Fourier
Transform.
I am currently reading about this, but as I have said, it is still over my
head.
I thought I would use your first example again to examine this concept.

Since this is for my own benefit as well, let me divert for a moment and
study the difference in Fourier in Excel and Mathematica.
In Excel, let's fill in A1:A8 with 1,2,3,0,0,0,0,0
And B1:B8 with 4,5,6,0,0,0,0,0

(Note: Excel can only work with 2^n numbers. However, Mathematica can work
with any size list.)

In Excel, make sure you have the analysis took pack loaded. Then, do
<Tools><Data Analysis> and take the Fourier of A1:A8 and put the results in
C1:C8.
The results I get are:

6
2.41421356237309-4.4142135623731i
-2-2i
-0.414213562373095+1.58578643762691i
2
-0.414213562373096-1.5857864376269i
-2+2i
2.4142135623731+4.41421356237309i


Let's do the same thing using Mathematica.


[Fourier[{1, 2, 3, 0, 0, 0, 0, 0}]]

{2.121320343559643, 0.8535533905932737 +
1.5606601717798214*I, -0.7071067811865476 +
0.7071067811865476*I, -0.14644660940672627 -
0.5606601717798213*I, 0.7071067811865476,
-0.14644660940672627 + 0.5606601717798213*I,
-0.7071067811865476 - 0.7071067811865476*I,
0.8535533905932737 - 1.5606601717798214*I}


Hmmm. The answers are different.
Here is a very short copy from Help." In different scientific and technical
fields different conventions are often used for defining discrete Fourier
transforms. The option FourierParameters in Mathematica allows you to choose
any of these conventions you want."

The default form of the equation this program uses are coefficients of. 1/
Sqrt(n).
"Help" goes on to talk about the field of signal processing, where the
coefficients are 1/n. (plus some other stuff.)

Let's set a variable to hold this option.

signal = FourierParameters -> {1, -1}

Fourier[{1,2,3,0,0,0,0,0},signal]

{6.,
2.41421-4.41421 I,
-2.-2.I,
-0.414214+1.58579 I,
2.,
-0.414214-1.58579 I,
-2.+2. I,
2.41421+4.41421 I}


Now the answers match. Although I have looked everywhere before, I could
not find any information on what convention Excel uses.
Therefore, it looks like Excel's Fourier function uses the definition as
used in the field of signal processing.

Let me test your theory of Convolution and the discrete Fourier transform.
Since the list should be 3+3-1 = 5 (5 numbers, which agrees with the list of
answers you gave) I changed the two lists to this.
I padded the ends with zeros.

v1 = {1, 2, 3, 0, 0};
v2 = {4, 5, 6, 0, 0};

I then multiplied the two Fourier's together, and took the inverse Fourier.
Now this is the interesting part. Taking the Inverse Fourier of the above
data returns the same Convolution that was done at the beginning.
InverseFourier[Fourier[v1,signal] * Fourier[v2,signal],signal]

{4.,13.,28.,27.,18.}

Wow! The answers agree.

Let's do this in Excel.
We have already done the transform of A1:A8 and placed the results in C1:C8.
Take the Transform of B1:B8 and place the results in D1:D8.
Now, in E1, enter the following engineering function:
=IMPRODUCT(C1,D1)
and copy this down E1:E8.

Now, do the Transform of E1:E8, except set the option to take the inverse
Transform and place the results F1:F8.

The results I get are:
4
13
28
27
18
0
0
0

Again, I may be wrong on my understanding. This was an interesting example
to try as I attempt to understand this concept.
I hope that someone with some real knowledge of this can jump in and add to
this discussion.


--
Dana DeLouis Windows Me & Office XP

<snip>


Curt Dye

unread,
Dec 9, 2001, 11:20:42 AM12/9/01
to
I wish I had more to say about your response, but I don't really know
anything about Fourier transforms. In any case, I'll reply to what I can:

>>>Are you trying to improve your VBA routine to do Convolution?

Yes.

>>>I am guessing that your abbreviation of d.f. is the Discrete Fourier
>>>Transform.

Actually, I meant d.f. as an abbreviation of density function.

I tried to follow your example, but I wasn't getting the same results. How
do we get the transform of the values in column A? Was that from
mathematica or Excel? I used Mathematica some in college, but I no longer
have a copy of it. (I have Maple around here somewhere but not at work....
anyway)
When I go to <Tools><Data Analysis><Fourier Analysis> I choose Input: A1:A8
=>> Ouput C1 and Input B1:B8 =>> Output D1, I get:

In Column C (as you showed)


6
2.41421356237309-4.4142135623731i
-2-2i
-0.414213562373095+1.58578643762691i
2
-0.414213562373096-1.5857864376269i
-2+2i
2.4142135623731+4.41421356237309i

and in Column D
15
7.53553390593273-9.53553390593274i
-2-5i
0.464466094067262+2.46446609406727i
5
0.46446609406726-2.46446609406726i
-2+5i
7.53553390593274+9.53553390593273i


The IMPRODUCT, however, yeilds
90
-23.8994949366118-56.2842712474619i
-6+14i
-4.10050506338836-0.284271247461902i
10
-4.10050506338832+0.284271247461908i
-6-14i
-23.8994949366115+56.2842712474619i

????

This is with A1:A8 = {1, 2, 3, 0, 0, 0, 0, 0}
Did your result come from the transforms of the 5 element array {1, 2, 3, 0,
0}? or what did I do wrong?

If this works, I could easily write the macro around this process. My other
options have all involved making matrices-- which is what I was trying to
aviod with the macro.


Dana DeLouis

unread,
Dec 9, 2001, 12:00:43 PM12/9/01
to
Yes, I have the same results as you do with column D, and (say column E) for
IMPRODUCT.
I think you are missing the "INVERSE" of the IMPRODUCT().
In Column F, I took the INVERSE of the Fourier Transform of Column E.
It is the same Excel procedure as Transform, except there is an option for
"Inverse."

I put the inverse in Column G and got these results.

4
13
28
27
18
0
0

0

> Did your result come from the transforms of the 5 element array {1, 2, 3,
0,
> 0}?

No, they were from the array of size 8 (1,2,3,0,0,0,0,0)
I suppose I could have used a size of 4.
I believe you have to have a size equal to a power of 2 for Excel to work.
(ie. 2,4,8,16,32...)

You can Take the transform of 3 different arrays and use IMPRODUCT on these
3 ranges, then do an inverse to do the same thing.
(I believe one of your examples was to do this 3 times.)
Write back, as this is an area I am trying to study, so this helps me also.

What was your definition of d.f.? I thought it was Discrete Fourier
Transform.


--
Dana DeLouis Windows Me & Office XP

"Curt Dye" <cd...@tampabay.rr.com> wrote in message
news:ubMQ7.28992$oj3.6...@typhoon.tampabay.rr.com...

Curt Dye

unread,
Dec 9, 2001, 9:28:29 PM12/9/01
to
You were correct-- all I was missing was the Inverse Transform of E.
Now, in certain cases, this works really well, but in others it doesn't work
for what I'm doing.

Before getting to that, however, by "d.f.", I really meant the density
function. In terms of probability, where X is a continous random variable,
the function f is called the density function of X if 1) f(x) >= 0 and 2) I
f(x) dx = 1 where I is the integral over all real values. Accordingly, I
believe I should have actually called F1 and F2 'cummulative density
functions' instead of density functions in my earlier example.

In any case, I have been messing around with this process and came up with a
macro for this process for a 32 element array: (As soon as I get happy with
it, I will make it more general)

Sub ConvolTester()
Range("C1:G32").ClearContents
Application.Run "ATPVBAEN.XLA!Fourier",
ActiveSheet.Range("$A$1:$A$32"), _
ActiveSheet.Range("$C$1"), False, False
Application.Run "ATPVBAEN.XLA!Fourier",
ActiveSheet.Range("$B$1:$B$32"), _
ActiveSheet.Range("$D$1"), False, False
Range("E1:E32").FormulaR1C1 = "=IMPRODUCT(RC[-2],RC[-1])"
Application.Run "ATPVBAEN.XLA!Fourier",
ActiveSheet.Range("$E$1:$E$32"), _
ActiveSheet.Range("$F$1"), True, False
For I = 1 To 32
Cells(I, 7).Formula = "=IMREAL(IMSUM($F$1:F" & I & "))"
Next I
End Sub

Following this, column G is now the "cummulative density function."

Using this, we get correct results unless A or B is a constant or near
constant. Using my example earlier of 1000 premium each year for 4 years,
and losses of .7, .2, and .1 .... if we fill in 0's for the tail, everything
works out. Now, however, what happens if we get 1000 premium for 32 years?
For some reason, this process tries to tell me that I lose 1000 per year,
which is obviously incorrect.

I also tried an example where A1 = 1 and each subsequent cell is 5% more
(i.e., A2 = A1*1.05, A3=A2*1.05, etc.). I left B1 = .7, B2 = .2, B3=.1
and B4:B32=0.

Strangely, I got:
C1 =75.2988293720722
D1=1
E1= 75.2988293720722
F1=2.03980213629671

Clearly F1 should be 1*.7 = .7 , but the other values are pretty good:

2.04
1.39
1.08
1.14
1.19
1.25
1.31
1.38
1.45
1.52
1.60
1.68
1.76
1.85
1.94
2.04
2.14
2.25
2.36
2.48
2.60
2.73
2.87
3.01
3.16
3.32
3.49
3.66
3.85
4.04
4.24
4.45

Using my previous code, I get

0.70
0.94
1.08
and then the same values as above.

I'm not sure what happens in this case.... Should the arrays be of order
2^n, where 2^n-1 < # of Non-Zero Values in A + # of Non-Zero Values in B - 1
< 2^n? It seems like the more zeros I put in, the better it does.

Dana DeLouis

unread,
Dec 12, 2001, 2:51:30 AM12/12/01
to
Hi Curt. I think we need a real Mathematician for anything more
complicated. However, I believe we're not doing to bad yet. :>)

I am not sure, but I think I know what is happening to your model.
In Column B, you have the numbers .7, .2, .1
You have also filled the Range A1:A32 with 1000.

I too get a column of 1,000's for results.
I believe the reason for this is that you have exceeded the size limit for
an array with only 32 elements.
When we tested 3 elements in column A, and 3 elements in column B, we ended
up with 3+3-1 =5 elements. However, Excel's Fourier requires an array size
rounded up to 8.

In your new example, I believe the Fourier size should now be 32+3-1=34, but
rounded up to 64!
(Trailing numbers zero.)

If you put 1000 in A1:A30, then the size of the array should be 30+3-1 = 32,
which should work.
And it appears to work correctly for me.
(ie 700, 900, 1000, 1000 ...etc)

What you saw earlier was probably some sort of cut-off.
I think if you up the array size to 64, then this should work for your
example. (keep in mind the formula)

Here is my attempt at the same macro idea you were using.
Anyway, I have not even looked at making your original macro more efficient.
I got side tracked by this interesting study.
Do you still want to study your original macro?
I am not sure, but this way may be able to handle 3 or more columns of data
using this same technique.

' = = = = = = = = = = = = = = = = = = = = = =
Sub ConvolTester_Version2()
Dim cell As Range
Dim NRows As Long

Const Fx_Fourier As String = "ATPVBAEN.XLA!Fourier"

Const Fourier As Boolean = False
Const InverseFourier As Boolean = True
Const Labels As Boolean = True

On Error Resume Next
With [A1:B64]
.SpecialCells(xlBlanks) = 0
NRows = .Rows.Count
End With
On Error GoTo 0

Columns("C:K").Clear

With Application
.Run Fx_Fourier, [A1].Resize(NRows, 1), [C1], Fourier, Not Labels
.Run Fx_Fourier, [B1].Resize(NRows, 1), [D1], Fourier, Not Labels
[E1].Resize(NRows, 1).FormulaR1C1 = "=IMPRODUCT(RC[-2],RC[-1])"
.Run Fx_Fourier, [E1].Resize(NRows, 1), [F1], InverseFourier, Not
Labels

'// Pause here at the inverse Fourier to round the results...
'// This also converts these strings back to numbers.
For Each cell In [F1].Resize(NRows, 1).Cells
cell = Round(cell, 5)
Next

'// Add a running total
[G1].Resize(NRows, 1).FormulaR1C1 = "=SUM(R1C[-1]:RC[-1])"
End With

Columns("C:G").Columns.AutoFit
End Sub

' = = = = = = = = = = = = = = = = = = = = = =

--
Dana DeLouis Windows Me & Office XP

"Curt Dye" <cd...@tampabay.rr.com> wrote in message

news:h5VQ7.43357$Ga5.7...@typhoon.tampabay.rr.com...

Christopher Eberhardt

unread,
Dec 12, 2001, 2:51:06 PM12/12/01
to
I tried a different approach than the one that you listed. Perhaps
this will help:

Here is my simplistic convolution function. Note that it assumes the
data is in columns and there is minimal error checking for clarity.
The first and second arguments are the data ranges that you want to
convolve, and the third argument is which result you want.

Function Conv(f1 As Range, f2 As Range, iDex As Integer) As Double
Dim dblResult As Double
Dim jf1 As Integer
Dim jf2 As Integer

jf2 = iDex
For jf1 = 1 To f1.Rows.Count
If (jf2 <= f2.Rows.Count) And (jf2 > 0) Then
dblResult = dblResult + f1.Cells(jf1, 1).Value *
f2.Cells(jf2, 1).Value
End If
jf2 = jf2 - 1
Next jf1
cmeConv = dblResult
End Function
A B C
1 Func1 Func2 Conv(f1,f2)
2 1 4 =Conv(Func1, Func2, 1)
3 2 5 =Conv(Func1, Func2, 2)
4 3 6 =Conv(Func1, Func2, 3)
5 =Conv(Func1, Func2, 4)
6 =Conv(Func1, Func2, 5)

Where Func1 is the named range $A$1:$A$4 and Func2 = $B$1:$B$4

Hope this helps

My example worksheet looks like

"Curt Dye" <cd...@tampabay.rr.com> wrote in message news:<DMBO7.181848$zK1.49...@typhoon.tampabay.rr.com>...

0 new messages