\n
![\"Unbenannt-3\"](\"http:\/\/benediktehinger.de\/blog\/science\/upload\/sites\/2\/2016\/03\/Unbenannt-3.png\")
<\/a><\/p>\nWe replace the top 10% and bottom 10% by the remaining most extreme value before calculating the mean (left panel). The right panel shows how the mean is influenced by a single outlier, but the winsorized mean is not (ignore the “yuen”-box”)<\/p>\n
Current Implementation<\/h3>\n
I adapted an implementation from the LIMO toolbox <\/a> based on Original Code from Prof. patrick J Bennett, McMaster University. In this code the dimension is fixed at dim = 3, the third dimension.<\/p>\n They solve it in three steps:<\/p>\n To generalize this to any dimension I have seen two previous solution that feels unsatisfied: Let’s solve it for X = 20 x 10 x 5 x 2 over the third dimension function [x] = winMean(x,dim,percent) % How long is the matrix in our required dimension % replace the left hand side beforeColons = num2cell(repmat(‘:’,dim-1,1)); \n Srep = thus subsref(x,Srep) outputs what x(:,:,2,:) would output. And then we need to repmat it, to fit the number of elements we replace by the winsorizing method.<\/p>\n This is put into subsasgn, where the S here is :<\/p>\n \n Srep = Thus equivalent to x(:,:,[1 2],:). The upper percentile is replaced analogous: x = subsasgn(x,S,repmat(subsref(x,Srep),[ones(1,dim-1) g+1 ones(1,nDim-dim)])); % general case<\/p>\n[\/matlab]\n And in the end we can take the mean, var, nanmean or whatever we need: That finishes the implementation.<\/p>\n But how about speed? I thus generated a random matrix of 200 x 10000 x 5 and measured the timing (n=100 runs) of the original limo implementation and mine:<\/p>\n For the last comparison I permuted the array prior to calculating the winsorized mean, thus the overhead. In my experience, the overhead is greater the larger the arrays are (I’m talking about 5-10GB matrices here).<\/p>\n My generalization seems to work fine. As expected it is slower than the hardcoded version. But it is faster than permuting the whole array.<\/p>\n","protected":false},"excerpt":{"rendered":" This is a function I wrote back in 2014. I think it illustrates an advanced functionality in matlab that I hadn’t found written about before. The problem: Calculate the winsorized mean of a multidimensional matrix over an arbitrary dimension. Winsorized Mean The benefits of the winsorized mean can be seen here: We replace the top 10% and bottom 10% by the remaining most extreme value before calculating the mean (left panel). The right panel shows how the mean is influenced by a single outlier, but the winsorized mean is not (ignore the “yuen”-box”) Current Implementation I adapted an implementation from…<\/p>\n","protected":false},"author":2,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[5],"tags":[],"class_list":["post-69","post","type-post","status-publish","format-standard","hentry","category-blog"],"_links":{"self":[{"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/posts\/69","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/users\/2"}],"replies":[{"embeddable":true,"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/comments?post=69"}],"version-history":[{"count":0,"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/posts\/69\/revisions"}],"wp:attachment":[{"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/media?parent=69"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/categories?post=69"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/tags?post=69"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}\n
\ng=floor((percent\/100)*n);
\nwx(:,:,1:g+1)=repmat(xsort(:,:,g+1),[1 1 g+1]);
\nwx(:,:,n-g:end)=repmat(xsort(:,:,n-g),[1 1 g+1]);
\n[\/matlab]\nGeneralisation<\/h3>\n
\n – Implement it for up to X dimension hardcoded and then use a switch-case to get the solution for the case.
\n – use permute to reorder the array and then go for the first dimension (which can be slow depending on the array)<\/p>\n
\n[matlab]\n
\n% x = matrix of arbitrary dimension
\n% dim = dimension to calculate the winsorized mean over
\n% percent = default 20, how strong to winsorize<\/p>\n
\nn=size(x,dim);
\n% number of items to winsorize and trim
\ng=floor((percent\/100)*n);
\nx=sort(x,dim);<\/p>\n[\/matlab]\nup to here it my and the original version are very similar. The hardest part is to generalize the part, where the entries are overwritten without doing it in a loop.
\nWe are now using the subsasgn command<\/a> and subsref <\/a>
\nWe need to generate a structure that mimics the syntax of
\n[matlab] x(:,:,1:g+1,:) = y [\/matlab]\nfor arbitrary dimensions and we need to construct y<\/em><\/p>\n[matlab]\n% Prepare Structs
\nSrep.type = ‘()’;
\nS.type = ‘()’;<\/p>\n
\nnDim = length(size(x));<\/p>\n
\nafterColons = num2cell(repmat(‘:’,nDim-dim,1));
\nSrep.subs = {beforeColons{:} [g+1] afterColons{:}};
\nS.subs = {beforeColons{:} [1:g+1] afterColons{:}};
\nx = subsasgn(x,S,repmat(subsref(x,Srep),[ones(1,dim-1) g+1 ones(1,nDim-dim)])); % general case
\n[\/matlab]\nThe output of Srep<\/em> is:<\/p>\n
\n type: ‘()’
\n subs: {‘:’ ‘:’ [2] ‘:’ }\n<\/p><\/blockquote>\n
\n type: ‘()’
\n subs: {‘:’ ‘:’ [1 2] ‘:’ }\n<\/p><\/blockquote>\n
\nThe evaluated structure then is:
\n[matlab] x(:,:,[1:2]) = repmat(x[:,:,1],[1 1 2 1]) [\/matlab]\n
\n[matlab]\n% replace the right hand side
\nSrep.subs = {beforeColons{:} [n-g] afterColons{:}};
\nS.subs = {beforeColons{:} [n-g:size(x,dim)] afterColons{:}};<\/p>\n
\n[matlab]\nx = squeeze(nanmean(x,dim));
\n[\/matlab]\nTiming<\/h3>\n
\n\n
\n algorithm<\/strong><\/th>\n timing (95% bootstraped CI of mean)<\/strong><\/th>\n<\/tr>\n<\/thead>\n\n \n limo_winmean<\/td>\n 185 – 188 ms<\/td>\n<\/tr>\n \n my_winmean<\/td>\n 202 – 203ms<\/td>\n<\/tr>\n \n limo_winmean otherDimension than 3 <\/td>\n 218 – 228 ms<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n Conclusion<\/h2>\n