# Efficient Weighted Moving average

88 views (last 30 days)
Matlab User on 13 Mar 2017
Edited: Jan on 21 Mar 2017
I have a matrix with x, y, and an uncertainty value (25x3 double). I want to compute a moving average of these points weighted by the third value, (V(:,3)), with a time window of say 5 points. What is the most efficient way to achieve this?
##### 2 CommentsShowHide 1 older comment
Matlab User on 13 Mar 2017
I have attached the input as V=x,y,uncertainty. The uncertainty is what I wish to weight the average by - the lower the value the higher the weight.

John BG on 16 Mar 2017
Edited: John BG on 16 Mar 2017
Hi Matlab User
1.
simulating data
clear all,clc
V=zeros(25,3);x=V(:,1);y=V(:,2);
w=randi([1 9],25,1);
w=w/sum(w);V(:,3)=w
V(:,1)=randi([-10 10],25,1) % x
V(:,2)=randi([-10 10],25,1) % y
V(:,3)=w;
.
the data
format short
V =
-3.0000 3.0000 0.0857
-5.0000 8.0000 0.0190
-10.0000 3.0000 0.0476
-2.0000 -1.0000 0.0476
4.0000 1.0000 0.0571
0 3.0000 0.0190
1.0000 -7.0000 0.0190
10.0000 -1.0000 0.0857
10.0000 -1.0000 0.0762
-6.0000 -3.0000 0.0762
-5.0000 -9.0000 0.0381
0 4.0000 0.0095
-7.0000 1.0000 0.0190
-8.0000 -10.0000 0.0286
3.0000 8.0000 0.0381
-4.0000 -4.0000 0.0095
2.0000 6.0000 0.0190
-9.0000 0 0.0095
5.0000 6.0000 0.0286
-4.0000 -2.0000 0.0286
-10.0000 6.0000 0.0667
-9.0000 -2.0000 0.0381
-9.0000 1.0000 0.0095
7.0000 4.0000 0.0476
9.0000 -6.0000 0.0762
2.
averaging with the shifting window and applying the weight V(:,:,3)
N=25
R=zeros(N,1)
for k=3:1:N-2
s=[k-2 k-1 k k+1 k+2]
R(k)= mean(.5*(V(s,2)+V(s,1)).*V(s,3))
end
R =
0
0
-0.0133
-0.0076
-0.0248
0.0857
0.1686
0.0714
0.0124
0.0276
-0.0610
-0.1810
-0.0705
-0.0248
-0.0133
-0.0105
0.0724
0.0133
-0.0057
-0.0629
-0.0619
-0.0410
-0.0010
0
0
This does not apply the shifting window on the first 2 and last 2 sames, there aren't 5 samples to apply the window on head and tail anyway, aren't they?
John BG
dpb on 20 Mar 2017
However is fine... :) I've noticed quite a few have typed 'dpd' rather than the trailing 'b'; I'm not sure why. I kinda' like the geometry of lower case string however as being mirror image/rotations of a single basic form...

dpb on 13 Mar 2017
Edited: dpb on 17 Mar 2017
N=5;
yavN=conv(v(:,2).*v(:,3),ones(1,N)/N,'valid');
doc conv % for other options, details...
If not clear how to use mean in the above, simply compute it "on the fly" as well instead of using just a single column from V...note the use of the optional DIM second argument to the MEAN function...
C=conv(mean(V(:,1:2),2).*V(:,3),ones(1,N)/N,'valid');
Jan on 20 Mar 2017
+1. This is much faster than the methods using loops.

Jan on 20 Mar 2017
Edited: Jan on 20 Mar 2017
Some further speed measurements:
function speedtest
N = 10000;
tic; for k = 1:N; R = test1(V); end, toc
tic; for k = 1:N; R = test2(V); end, toc
tic; for k = 1:N; R = test3(V); end, toc
tic; for k = 1:N; R = test4(V); end, toc
function R = test1(V) % John BG
N=size(V, 1);
R=zeros(N,1);
for k=3:1:N-2
s=[k-2 k-1 k k+1 k+2];
R(k)= mean(.5*(V(s,2)+V(s,1)).*V(s,3));
end
function R = test2(V) % Jan
% 1. repeated calcultations moved outside the loop
% 2. Indexing by a:b is much faster than using a pre-calculated vector
N = size(V, 1);
R = zeros(N,1);
X = 0.5 * sum(V(:, 1:2), 2) .* V(:, 3);
for k = 3:1:N-2
R(k) = sum(X(k-2:k+2)) * 0.2;
end
function R = test3(V) % dpb
N = 5;
R = conv(mean(V(:,1:2), 2) .* V(:,3), ones(1,N) / N,'valid');
function R = test4(V) % dpb without MEAN()
N = 5;
R = conv(0.5 * (V(:,1) + V(:, 2)) .* V(:,3), ones(1,N) / N, 'valid');
Results (R2016b/64, Core2Duo):
Elapsed time is 3.138929 seconds. John BG
Elapsed time is 0.444694 seconds. Improved loop
Elapsed time is 0.193468 seconds. dpb
Elapsed time is 0.119400 seconds. dpb without MEAN()
What happens, if the input is not tiny:
V = rand(25000, 3);
N = 10;
Elapsed time is 3.667920 seconds.
Elapsed time is 0.418300 seconds.
Elapsed time is 0.018349 seconds.
Elapsed time is 0.014222 seconds.
Conclusion: Loops profit from smart indexing and avoiding repeated calculations, but conv rules. mean has a remarkable overhead.
##### 2 CommentsShowHide 1 older comment
Jan on 21 Mar 2017
x = rand(1000, 2);
tic; for k = 1:100000; m = x(:, 1) + x(:, 2); end; toc
tic; for k = 1:100000; m = sum(x, 2); end; toc
Elapsed time is 0.120969 seconds.
Elapsed time is 0.156458 seconds.
The difference can be measured. SUM is parallelized for more than 88999 elements. Then the number of cores matter also.
x = rand(100000, 2);
tic; for k = 1:1000; m = x(:, 1) + x(:, 2); end; toc
tic; for k = 1:1000; m = sum(x, 2); end; toc
Elapsed time is 0.103634 seconds.
Elapsed time is 0.148599 seconds.
But rmember, we are talking about 0.00004 seconds. "sum(x,2)" can save minute of hourse during the debugging. Therefore I prefer sum and apply the dirty tricks like v(:,1)+v(:,2) only for bottlenecks and with an exhaustive documentation and a unit-test function, which specifically cares about these lines.