Computationally, all we're doing is convolution. We are computing A*B*C*...*J where * here means convolution not multiply. We can do this with the R convolution routine, and it's much faster because convolution is performed internally with compiled algorithms that do fast Fourier transforms. FFTs have complexity of n*log(n) instead of n^2. R defines convolution a little screwy, so we need the first line to define it the way we want it. So we want to basically do this:
Code:
conv = function(f,g) convolve(f, rev(c(0, g)), type = "o")
a = conv(A,B)
b = conv(a,C)
c = conv(b,D)
d = conv(c,E)
e = conv(d,F)
f = conv(e,G)
g = conv(f,H)
h = conv(g,I)
i = conv(h,J)
i
I made this better by putting all of your probabilities in a single matrix and looping through it instead of having separate letters. Then you can do different numbers of integers by just changing the matrix instead of the code itself. You can edit that in the code, or separately if you type edit(M) which brings up a spreadsheet type thing. So now the code is just a few lines aside from the matrix, and it responds instantaneously.
Code:
M = matrix(c(
.0405, .0374, .1468, .1063, .2368, .1265, .1630, .0445, .0628, .0121, .0192, .0020, .0000, .0020,
.0427, .0305, .1442, .0998, .2329, .1109, .1803, .0555, .0577, .0144, .0200, .0089, .0022, .0000,
.0427, .0305, .1442, .0998, .2329, .1109, .1803, .0555, .0577, .0144, .0200, .0089, .0022, .0000,
.0351, .0229, .1340, .1024, .2378, .1003, .1827, .0673, .0566, .0172, .0265, .0129, .0043, .0000,
.0351, .0229, .1340, .1024, .2378, .1003, .1827, .0673, .0566, .0172, .0265, .0129, .0043, .0000,
.0297, .0222, .1342, .0941, .2350, .1097, .1853, .0675, .0586, .0148, .0319, .0126, .0044, .0000,
.0297, .0222, .1342, .0941, .2350, .1097, .1853, .0675, .0586, .0148, .0319, .0126, .0044, .0000,
.0253, .0214, .1404, .0673, .2495, .0897, .1979, .0750, .0780, .0205, .0283, .0049, .0019, .0000,
.0253, .0214, .1404, .0673, .2495, .0897, .1979, .0750, .0780, .0205, .0283, .0049, .0019, .0000,
.0284, .0262, .1223, .0306, .2314, .1092, .2162, .0742, .1092, .0197, .0240, .0087, .0000, .0000),
ncol=14, byrow=TRUE)
conv = function(f,g) convolve(f, rev(c(0, g)), type = "o")
c = conv(M[1,],M[2,])
if (nrow(M) > 2) for (i in 3:nrow(M)) c = conv(c,M[i,])
c[1:(nrow(M)-1)] = 0
plot(c,xlab='Sum',ylab='Probability')
c
One thing about this is that numbers that are 0 or very very tiny won't be exact due to roundoff error, and some will be very small and negative. For example, for 140 you get -7.137148e-18. But that only affects numbers that are around 1e-17, so you probably don't care. We could just round those to 0 too. I added a line to force 1-9 to be 0 in this case.
I added a plot to the code. Here it is:
BTW, if you want to check that the probabilities add to 1, you can do for example, sum(A) with the old code, or with the new code sum(M[1,]), sum(M[2,]), etc. sum(M) should sum the whole matrix and give 10 right now, except it gives 10.0002.
Last edited by BruceZ; 04-10-2014 at 03:26 AM.
Reason: Added plot to code