Efficient Weighted Moving average

88 views (last 30 days)
Matlab User
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 Comments
Matlab User
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.

Sign in to comment.

Accepted Answer

John BG
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?
if you find this answer useful would you please be so kind to mark my answer as Accepted Answer?
To any other reader, please if you find this answer
please click on the thumbs-up vote link
thanks in advance
John BG
  9 Comments
dpb
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...
"Addendum. An addendum, in general, is an addition required to be made to a document by its author subsequent to its printing or publication. It comes from the Latin verbal phrase addendum est, being the gerundive form of the verb addo, addere, addidi, additum, "to give to, add to", meaning "(that which) must be added"."

Sign in to comment.

More Answers (2)

dpb
dpb on 13 Mar 2017
Edited: dpb on 17 Mar 2017
N=5;
yavN=conv(v(:,2).*v(:,3),ones(1,N)/N,'valid');
discarding the zero-padded end values. See
doc conv % for other options, details...
ADDENDUM
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');
  1 Comment
Jan
Jan on 20 Mar 2017
+1. This is much faster than the methods using loops.

Sign in to comment.


Jan
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 Comments
Jan
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.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!