Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

To resolve issues starting MATLAB on Mac OS X 10.10 (Yosemite) visit: http://www.mathworks.com/matlabcentral/answers/159016

Calculate Sn (like median absolute deviation) using matlab.

Asked by Witold on 18 Apr 2013

Hi, I'm trying to compute Sn value based on : "Alternatives to the Median Absolute Deviation Peter J. Rousseeuw and Christophe Croux" paperwork.

According to this article:

Sn = c*median_i{median_j|xi - xj|}

For each i we compute the median of {|xi-xj|; j = 1,...,n}. This yields n numbers, the median of which gives our final estimate Sn. C is constant value.

I have Y matrix of size 10000 x 2000 (i = 2000, j = 10000) Because of the size, the simplest loop method:

for i = 1 : size(Y,2)
    for j = 1 : size(Y,1)
        a(j) = median( abs( Y(:,i)-Y(j,i) ) );
    end
    Sn(i) = c * median(a);
end

is time consuming, so is no good at all. I'm matlab newbie and I don't know how to use _repmat _ function which will be - as I guess - very helpful. Can I ask you for help, how to write this so that the computation time will be comparable to _mad _ function?

0 Comments

Witold

Products

No products are associated with this question.

2 Answers

Answer by Andrei Bobrov on 18 Apr 2013
Edited by Andrei Bobrov on 18 Apr 2013
Accepted answer
s = [3 1130;
     4 1527;
     3  907;
     2  878;
     4  995];
k = permute(abs( bsxfun(@minus,s,permute(s,[3 2 1])) ),[1 3 2]);
out = squeeze(median(median(k)));

ADD

out = zeros(1,size(s,2));
for jj = 1:size(s,2)
    out(jj) = median(median(abs(bsxfun(@minus,s(:,jj),s(:,jj)'))));
end

ADD2 use function mad from Statistics Toolbox

out = mad(s,1);

2 Comments

Witold on 18 Apr 2013

Thank you very much, that's what I was thinking about, unfortunately I got the error message:

Error using bsxfun
Out of memory. Type HELP MEMORY for your options. 

memory available: win 8 x64 / 8GB/ r1012b

Maximum possible array:      6163 MB (6.462e+09 bytes) *
Memory available for all arrays:      6163 MB (6.462e+09 bytes) *
Memory used by MATLAB:      1009 MB (1.058e+09 bytes)
Physical Memory (RAM):      8157 MB (8.553e+09 bytes)
*  Limited by System Memory (physical + swap file) available.

should I give up? or mayby try to split columns of the matrix

Andrei Bobrov on 18 Apr 2013

see ADD

Andrei Bobrov
Answer by Image Analyst on 18 Apr 2013

I think you're totally misinterpreting the formula. You don't look at the difference of all possible pairs of points - that would take forever and isn't what the formula says. Look at Wikipedia: http://en.wikipedia.org/wiki/Median_absolute_deviation For each element, you'd have a million differences. And there are a million elements so you'd have a billion values. No wonder it takes forever! The i and j DO NOT REFER TO ROW AND COLUMN. They refer to element, or "linear index" in MATLAB lingo.

I think what they mean is to take the mean over all elements (pixels) and then take the median of the differences, which is totally different. So

% Get the median value of all elements (pixels).
medianOfWholeMatrix = median(Y(:));
% Find the difference between each element and the overall median.
differenceImage = double(Y) - medianOfWholeMatrix;
% Now take the median of that:
madValue = median(differenceImage (:));

This should be substantially faster, not to mention correct.

3 Comments

Witold on 18 Apr 2013

the output (Sn) should be a vector of size 1 x n where n is num. of columns

Image Analyst on 18 Apr 2013

Their adaption of MAD might do something different. You could do the same thing on a column-by-column basis. If you want to do that, you could just extract a column at a time and put it in a loop

for col = 1 : size(Y, 2)
  % Get the median value of all elements (pixels) in one column.
  medianOfWholeColumn = median(Y(:, col));
  % Find the difference between each element and the overall median.
  differenceImage = double(Y(:, col)) - medianOfWholeColumn ;
  % Now take the median of that:
  madValue(col) = median(differenceImage (:));
end
Sn = c * madValue;
Witold on 18 Apr 2013

I'm afraid my interpretation is correct, I founded an example of Sn algorithm: http://www.maplesoft.com/support/help/AddOns/view.aspx?path=Statistics/RousseeuwCrouxSn

consider M 5x2 matrix:

M = [3 1130;
     4 1527;
     3  907;
     2  878;
     4  995]

for each column we calculate x(:,i) - x(i,j) and its median, for 1st column it is:

abs(x(:,1) - x(1,1)) = [0;1;0;1;1] -> med(..) = [1]
abs(x(:,1) - x(2,1)) = [1;0;1;2;0] -> med(..) = [1]
abs(x(:,1) - x(3,1)) = [0;1;0;1;1] -> med(..) = [1]
abs(x(:,1) - x(4,1)) = [1;2;1;0;2] -> med(..) = [1]
abs(x(:,1) - x(5,1)) = [1;0;1;2;0] -> med(..) = [1]

now we calculate median of those medians which is 1, and we move to second column - using my simplest loops we get 117 which is correct answer.

I agree that their adaption of MAD is something different, I was just wondering if there is any method to calculate it faster... it takes up to 5 minutes to compute it for 1000 x ~2000 matrix, which is ok if worse comes to worst, but much over hours(?) for size 10k /: Anyway, thank you for help!

Image Analyst

Contact us