Smirnov test for 2 distributions [H,prob,d] = smirnov(x1 , x2 , alpha, iflag ) Written by Marco Ratto Joint Research Centre, The European Commission, (http://eemc.jrc.ec.europa.eu/), marco.ratto@jrc.it Reference: M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
0001 function [H,prob,d] = smirnov(x1 , x2 , alpha, iflag ) 0002 % Smirnov test for 2 distributions 0003 % [H,prob,d] = smirnov(x1 , x2 , alpha, iflag ) 0004 % 0005 % Written by Marco Ratto 0006 % Joint Research Centre, The European Commission, 0007 % (http://eemc.jrc.ec.europa.eu/), 0008 % marco.ratto@jrc.it 0009 % 0010 % Reference: 0011 % M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006. 0012 0013 % Copyright (C) 2012 Dynare Team 0014 % 0015 % This file is part of Dynare. 0016 % 0017 % Dynare is free software: you can redistribute it and/or modify 0018 % it under the terms of the GNU General Public License as published by 0019 % the Free Software Foundation, either version 3 of the License, or 0020 % (at your option) any later version. 0021 % 0022 % Dynare is distributed in the hope that it will be useful, 0023 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0024 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0025 % GNU General Public License for more details. 0026 % 0027 % You should have received a copy of the GNU General Public License 0028 % along with Dynare. If not, see <http://www.gnu.org/licenses/>. 0029 0030 if nargin<3 0031 alpha = 0.05; 0032 end 0033 if nargin<4, 0034 iflag=0; 0035 end 0036 0037 % empirical cdfs. 0038 xmix= [x1;x2]; 0039 bin = [-inf ; sort(xmix) ; inf]; 0040 0041 ncount1 = histc (x1 , bin); 0042 ncount1 = ncount1(:); 0043 ncount2 = histc (x2 , bin); 0044 ncount2 = ncount2(:); 0045 0046 cum1 = cumsum(ncount1)./sum(ncount1); 0047 cum1 = cum1(1:end-1); 0048 0049 cum2 = cumsum(ncount2)./sum(ncount2); 0050 cum2 = cum2(1:end-1); 0051 0052 n1= length(x1); 0053 n2= length(x2); 0054 n = n1*n2 /(n1+n2); 0055 0056 % Compute the d(n1,n2) statistics. 0057 0058 if iflag==0, 0059 d = max(abs(cum1 - cum2)); 0060 elseif iflag==-1 0061 d = max(cum2 - cum1); 0062 elseif iflag==1 0063 d = max(cum1 - cum2); 0064 end 0065 % 0066 % Compute P-value check H0 hypothesis 0067 % 0068 0069 lam = max((sqrt(n) + 0.12 + 0.11/sqrt(n)) * d , 0); 0070 if iflag == 0 0071 j = [1:101]'; 0072 prob = 2 * sum((-1).^(j-1).*exp(-2*lam*lam*j.^2)); 0073 0074 prob=max(prob,0); 0075 prob=min(1,prob); 0076 else 0077 prob = exp(-2*lam*lam); 0078 end 0079 0080 H = (alpha >= prob);