/
limo_trimmed_mean.m
118 lines (94 loc) · 3.25 KB
/
limo_trimmed_mean.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
function TM = limo_trimmed_mean(varargin)
% Compute a trimmed mean and CI on the 3rd dimension of a 3D matrix
%
% INPUT
% TM = limo_trimmed_mean(data,percent, alpha)
% data is a 3D matrix
% options are percent [0-100] the percentage of trimming to do
% e.g. 20 will trim 20% on each side, thus a total of 40% of the data.
% alpha [0-100] is the error rate of the confidence interval
% if alpha is not used, only the trimmed mean is computed
%
% OUTPUT
% TM is a 2D or 3D matrix with the lower CI, the trimmed mean and the high CI
% ------------------------------
% Copyright (C) LIMO Team 2019
% Original R code by Rand Wilcox - See also Wilcox p.71, 139
% GA Rousselet, University of Glasgow, Dec 2007
% C. Pernet rewrote for 3D matrices
% Version 1 June 2010
%% checkings
data = varargin{1};
percent = 20/100;
if nargin > 1
percent = varargin{2};
if percent > 0
percent = percent / 100;
end
end
if nargin == 3
option = 1;
alpha = varargin{3};
else
option = [];
end
reduced_dim = 0;
if length(size(data)) ==2
tmp = zeros(2,size(data,1),size(data,2));
tmp(1,:,:) = data; tmp(2,:,:) = data;
clear data; data = tmp; clear tmp
reduced_dim = 1;
end
if length(size(data)) ~=3
error('data in must be a 2D or 3D matrix')
end
%% do the trimming
n = size(data,3);
g = floor(percent*n); % g trimmed elements
datasort = sort(data,3);
if option == 1
TM = NaN(size(data,1),size(data,2),3);
TM(:,:,2) = nanmean(datasort(:,:,(g+1):(n-g)),3);
else
TM = NaN(size(data,1),size(data,2));
TM = nanmean(datasort(:,:,(g+1):(n-g)),3);
end
%% compute confidence intervals
if option == 1
for i=1:size(data,1)
tmp_data = squeeze(datasort(i,:,:)); tmp = tmp_data(~isnan(tmp_data));
[tv,g] = tvar(reshape(tmp,size(tmp_data,1),length(tmp)/size(tmp_data,1)),percent); % trimmed squared standard error + g trimmed elements
se = sqrt(tv); % trimmed standard error
df = n - 2.*g - 1; % n-2g = number of observations left after trimming
TM(i,:,1) = TM(i,:,2)+tinv(alpha./2,df).*se; %
TM(i,:,3) = TM(i,:,2)-tinv(alpha./2,df).*se; %
end
end
if reduced_dim == 1
TM = TM(1,:,:);
end
% t=(tmdata-nullvalue)./se;
% p=2.*(1-tcdf(abs(t),df)); % 2-tailed probability
% tcrit=tinv(1-alphav./2,df); % 1-alpha/2 quantile of Student's distribution with df degrees of freedom
end
function [tv,g]=tvar(x,percent)
% function [tv]=tvar(x,percent)
% The trimmed variance is calculated from the winsorized variance, vw,
% according to the formula vw/(k*length(x)), where k=(1-2*percent/100)^2.
% Original code provided by Prof. Patrick J. Bennett, McMaster University
% See Rand R. Wilcox (2001), Fundamentals of Modern Statisical Methods, page 164.
% See also Rand R. Wilcox (2005), p.61-63
% Edit input checks: GA Rousselet - University of Glasgow - Nov 2008
% Merge function tvar, winvar, winsample - make it work 3D - C Pernet June 2010
g=floor((percent/100)*size(x,2));
xsort=sort(x,2);
loval=xsort(:,g+1);
hival=xsort(:,size(x,2)-g);
for i=1:size(x,1)
x(i,find(x(i,:)<=loval(i))) = loval(i);
x(i,find(x(i,:)>=hival(i))) = hival(i);
wv(i)=var(x(i,:));
end
k=(1-2*percent/100)^2;
tv=wv/(k*length(x));
end