sumskipnan.m 5.78 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
function [o,count,SSQ] = sumskipnan(x, DIM, W)
% SUMSKIPNAN adds all non-NaN values.
%
% All NaN's are skipped; NaN's are considered as missing values.
% SUMSKIPNAN of NaN's only  gives O; and the number of valid elements is return.
% SUMSKIPNAN is also the elementary function for calculating
% various statistics (e.g. MEAN, STD, VAR, RMS, MEANSQ, SKEWNESS,
% KURTOSIS, MOMENT, STATISTIC etc.) from data with missing values.
% SUMSKIPNAN implements the DIMENSION-argument for data with missing values.
% Also the second output argument return the number of valid elements (not NaNs)
%
% Y = sumskipnan(x [,DIM])
% [Y,N,SSQ] = sumskipnan(x [,DIM])
% [...] = sumskipnan(x, DIM, W)
%
16
17
18
19
20
21
22
% x     input data
% DIM   dimension (default: [])
%       empty DIM sets DIM to first non singleton dimension
% W     weight vector for weighted sum, numel(W) must fit size(x,DIM)
% Y     resulting sum
% N     number of valid (not missing) elements
% SSQ   sum of squares
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
%
% the function FLAG_NANS_OCCURED() returns whether any value in x
%  is a not-a-number (NaN)
%
% features:
% - can deal with NaN's (missing values)
% - implements dimension argument.
% - computes weighted sum
% - compatible with Matlab and Octave
%
% see also: FLAG_NANS_OCCURED, SUM, NANSUM, MEAN, STD, VAR, RMS, MEANSQ,
%      SSQ, MOMENT, SKEWNESS, KURTOSIS, SEM


%    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 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program; If not, see <http://www.gnu.org/licenses/>.

50
%       $Id: sumskipnan.m 9033 2011-11-08 20:58:07Z schloegl $
Stéphane Adjemian's avatar
Stéphane Adjemian committed
51
%
52
%       Copyright (C) 2000-2005,2009,2011 by Alois Schloegl <alois.schloegl@gmail.com>
Stéphane Adjemian's avatar
Stéphane Adjemian committed
53
%       Copyright (C) 2014-2017 Dynare Team
54
55
56
57
58
59
%       This function is part of the NaN-toolbox
%       http://pub.ist.ac.at/~schloegl/matlab/NaN/


global FLAG_NANS_OCCURED;

60
if nargin<2
61
    DIM = [];
62
63
end
if nargin<3
64
    W = [];
65
end
66
67
68
69
70
71
72
73
74
75
76

% an efficient implementation in C of the following lines
% could significantly increase performance
% only one loop and only one check for isnan is needed
% An MEX-Implementation is available in sumskipnan.cpp
%
% Outline of the algorithm:
% for { k=1,o=0,count=0; k++; k<N}
%   if ~isnan(i(k))
%   {   o     += x(k);
%               count += 1;
77
78
79
%               tmp    = x(k)*x(k)
%               o2    += tmp;
%               o3    += tmp.*tmp;
80
81
%       };

82
if isempty(DIM)
83
84
    DIM = find(size(x)>1,1);
    if isempty(DIM), DIM = 1; end
85
end
86
if (DIM<1), DIM = 1; end %% Hack, because min([])=0 for FreeMat v3.5
87
88
89
90

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% non-float data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
91
if  (isempty(W) && (~(isa(x,'float') || isa(x,'double')))) || ~flag_implicit_skip_nan() %%% skip always NaN's
92
93
    if ~isempty(W)
        error('SUMSKIPNAN: weighted sum of integers not supported, yet');
94
    end
95
96
97
98
99
100
101
102
103
104
    x = double(x);
    o = sum(x,DIM);
    if nargout>1
        sz = size(x);
        N  = sz(DIM);
        sz(DIM) = 1;
        count = repmat(N,sz);
        if nargout>2
            x = x.*x;
            SSQ = sum(x,DIM);
105
106
107
108
        end
    end
    return
end
109
110
111
112
113

if (length(size(x))<DIM)
    error('SUMSKIPNAN: DIM argument larger than number of dimensions of x');
elseif ~isempty(W) && (size(x,DIM)~=numel(W))
    error('SUMSKIPNAN: size of weight vector does not match size(x,DIM)');
114
end
115
116
117
118
119
120
121
122
123
124
125
126
127

%% mex and oct files expect double
x = double(x);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% use Matlab-MEX function when available
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%if 1,
try

    %% using sumskipnan_mex.mex

    %% !!! hack: FLAG_NANS_OCCURED is an output argument, reserve memory !!!
128
    if isempty(FLAG_NANS_OCCURED)
129
        FLAG_NANS_OCCURED = logical(0);  % default value
130
    end
131

132
    if (nargout<2)
133
134
135
136
        o = sumskipnan_mex(real(x),DIM,FLAG_NANS_OCCURED,W);
        if (~isreal(x))
            io = sumskipnan_mex(imag(x),DIM,FLAG_NANS_OCCURED,W);
            o  = o + i*io;
137
138
139
        end
        return
    elseif (nargout==2)
140
141
142
143
144
145
146
        [o,count] = sumskipnan_mex(real(x),DIM,FLAG_NANS_OCCURED,W);
        if (~isreal(x))
            [io,icount] = sumskipnan_mex(imag(x),DIM,FLAG_NANS_OCCURED,W);
            if any(count(:)-icount(:))
                error('Number of NaNs differ for REAL and IMAG part');
            else
                o  = o+i*io;
147
148
149
150
            end
        end
        return
    elseif (nargout>=3)
151
152
153
154
155
156
157
158
        [o,count,SSQ] = sumskipnan_mex(real(x),DIM,FLAG_NANS_OCCURED,W);
        if (~isreal(x))
            [io,icount,iSSQ] = sumskipnan_mex(imag(x),DIM,FLAG_NANS_OCCURED,W);
            if any(count(:)-icount(:))
                error('Number of NaNs differ for REAL and IMAG part');
            else
                o  = o+i*io;
                SSQ = SSQ+iSSQ;
159
160
161
162
163
            end
        end
        return
    end
end
164
165
166

if ~isempty(W)
    error('weighted sumskipnan requires sumskipnan_mex');
167
end
168
169
170
171

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% count non-NaN's
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
172
if nargout>1
173
    count = sum(x==x,DIM);
174
    FLAG_NANS_OCCURED = any(count(:)<size(x,DIM));
175
end
176
177
178
179
180
181
182

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% replace NaN's with zero
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x(x~=x) = 0;
o = sum(x,DIM);

183
if nargout>2
184
185
    x = real(x).^2 + imag(x).^2;
    SSQ = sum(x,DIM);
186
end
187
188
189
190
191
192

%!assert(sumskipnan([1,2],1),[1,2])
%!assert(sumskipnan([1,NaN],2),1)
%!assert(sumskipnan([1,NaN],2),1)
%!assert(sumskipnan([nan,1,4,5]),10)
%!assert(sumskipnan([nan,1,4,5]',1,[3;2;1;0]),6)