How to Calculate Moving Product?

2 次查看（过去 30 天）
Gaurav Soni
Gaurav Soni 2015-12-18
评论： Gaurav Soni 2020-1-8
Is there any vectorized way of calculating moving products. If A is a vector, how can I calculate its moving product with a specified window size.
For example, for a window size of 3, the moving product should be equal to:
[A(3)*A(2)*A(1), A(4)*A(3)*A(2), A(5)*A(4)*A(3), ..., A(N)*A(N-1)*A(N-2)]
How do I calculate it without using a for loop. I know that the 'filter' function can be used for calculating moving means in a vectorized manner; but I have not seen any function that can calculate a moving product.
Can someone please help. Thanks.
  1 个评论
Guido Marco
Guido Marco 2016-6-27
try somthing like this... function Y=divf(X) % divisive differential
X0=X(1:end-1); X1=X(2:end); Y=X1./X0;

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Jonathan Agg
Jonathan Agg 2020-1-6
The function movprod was added in 17a. There's also a list of similar functions here which are included in MATLAB, like movsum and movmax.

更多回答（3 个）

Guillaume
Guillaume 2015-12-18
One way to achieve what you want:
function mp = movingproduct(v, windowsize)
%add input validation code
vv = repmat([v(:); nan], 1, windowsize);
vv = reshape(vv(1:end-windowsize), [], windowsize);
mp = prod(vv(windowsize:end, :), 2);
if isrow(v), mp = mp.'; end
end
  2 个评论
Gaurav Soni
Gaurav Soni 2015-12-18
Thank you. This is a nice method to avoid a For loop. I was also able to find a way. Here it is:
function mp = movingproduct(v, windowsize)
logv = log(v);
b = (1/windowsize)*ones(1,windowsize);
mp = exp(windowsize*filter(b,1,logv));
mp = mp(windowsize:end);
end
This method is a bit faster than what you devised above, but it suffers from numerical errors that accumulate due to log and exp operations.
Guillaume
Guillaume 2015-12-19
编辑：Guillaume 2016-5-21
Another possible implementation, which has no numerical errors. It does have a loop (through arrayfun) but only over the window size.
function mp = movingproduct(v, windowsize)
%add input validation code
vv = repmat({v(:)}, 1, windowsize);
vv = arrayfun(@(row) circshift(vv{row}, row-1), 1:windowsize, 'UniformOutput', false);
vv = [vv{:}];
mp = prod(vv(windowsize:end, :), 2);
if isrow(v), mp = mp.'; end
end

请先登录，再进行评论。

Image Analyst
Image Analyst 2015-12-19
One possibility is to use nlfilter() where you can define the filter to be used on the window to do absolutely anything you want it to do, such as prod().
output = nlfilter(yourMatrix, [1, 3], @prod);
I attach a full blown demo where I use nlfilter to compute the Otsu threshold on a sliding 3x3 window basis.
  1 个评论
Gaurav Soni
Gaurav Soni 2015-12-19
编辑：Gaurav Soni 2015-12-19
Thanks. I think this would be pretty useful for people who have access to Image Processing toolbox. It looks like nlfilter is only available in the Image Processing Toolbox.

请先登录，再进行评论。

Mark Britten-Jones
Mark Britten-Jones 2016-5-21
编辑：Guillaume 2016-5-21
The problem with the loop approach is that it repeats calculations. For example A(3)*A(2) is calculated twice in a moving product of window size 3. The vectorised versions also suffer from this repetition. It can be avoided using a recursive approach which recognises that e.g. for k even a k moving product can be calculates as a k/2 moving product times the k/2 lag of the k/2 moving product. For large k the increases in speed are dramatic. Here is the code, which relies on a lag function (few lines to code up):
function y = cumprodk(x,k)
y = x;
if k==1; return; end;
y = cumprodk(y, floor(k/2));
y = lag(y, floor(k/2)).*y;
if rem(k,2) == 1 % if k odd need to multiply by lag k-1
y = lag(x,k-1).*y;
end
end
Thats it. No loops and faster for large k

类别

在 Help Center 和 File Exchange 中查找有关 Install Products 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by