/
quantile_new.m
72 lines (71 loc) · 2.18 KB
/
quantile_new.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
function [y,b]=quantile(X,u,dim,mode)
% quantile - Quantile of a vector
% [y] = quantile(X,u)
% With 0 <= u <= 1 : Find the value(s) 'y' so that:
% u = [ The proportion of values in X <= y ]
% If X is a matrix, the quantile will be taken along the 1st dimension
% If u is a matrix, outputs numel(u) values
%
% [y] = quantile(X,u,dim) along dimension dim
% [y,b]=quantile(X,u) also outputs a logical array of those values below
% the quantile.
%
% Examples
% >> quantile(X,0.5) outputs the median
% >> quantile(X,0.95) outputs the level of the highest 5%
%
% Author: K. N'Diaye (kndiaye01<at>yahoo.fr)
% Copyright (C) 2006
% This program is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the
% Free Software Foundation; either version 2 of the License, or (at your
% option) any later version: http://www.gnu.org/copyleft/gpl.html
%
% ----------------------------- Script History ---------------------------------
% KND 2006-01-05 Creation
% KND 2008-10-13 NxMxPx... matrix outputs a 1xMxPx... matrix
% KND 2009-03-05 Accepts list vector in u
% ----------------------------- Script History ---------------------------------
u=u(:);
if nargin<2
error('Threshold needed!')
end
sX=size(X);
if nargin<3
if prod(sX)==max(sX)
X=X(:);
else
error('Choose a dimension on which to compute quantiles');
end
else
pX=[dim setdiff(1:ndims(X),dim)];
X=permute(X,pX);
end
[X,i]=sort(X);
q1 = min(floor(u*sX(1)+1),sX(1));
q2 = max(ceil(u*sX(1)),1);
y1 = X(q1,:);
y2 = X(q1,:);
y =(y1+y2)/2;
if prod(sX)~=max(sX) % X wasn't a vector
nu = numel(u);
%if nu>1
% y is made into a P-row matrix (1-by-N-by-...)
sX(dim)=1;
y=reshape(y, [nu sX 1]);
%else
% y=reshape(y, [sX(2:end) 1]);
%end
end
if nargout > 1
b = i < round((q1+q2)/2);
if prod(sX)~=max(sX) % X wasn't a vector
nu = numel(u);
if nu>1
b=reshape(b, [nu sX(2:end) 1]);
else
b=reshape(b, [sX(2:end) 1]);
end
end
b=permute(b,[2:ndims(b) 1]);
end