new utilities to avoid Statistics Toolbox
git-svn-id: file:///home/sebastien/dynare/gsa_dyn@28 f1850c17-3b45-254b-b221-fcb05880fee1time-shift
rattoma
parent
d00551791c
commit
2cdd54e173
|
|
@ -0,0 +1,38 @@
|
|||
function x = beta_inv(p, a, b)
|
||||
% PURPOSE: inverse of the cdf (quantile) of the beta(a,b) distribution
|
||||
%--------------------------------------------------------------
|
||||
% USAGE: x = beta_inv(p,a,b)
|
||||
% where: p = vector of probabilities
|
||||
% a = beta distribution parameter, a = scalar
|
||||
% b = beta distribution parameter b = scalar
|
||||
% NOTE: mean [beta(a,b)] = a/(a+b), variance = ab/((a+b)*(a+b)*(a+b+1))
|
||||
%--------------------------------------------------------------
|
||||
% RETURNS: x at each element of p for the beta(a,b) distribution
|
||||
%--------------------------------------------------------------
|
||||
% SEE ALSO: beta_d, beta_pdf, beta_inv, beta_rnd
|
||||
%--------------------------------------------------------------
|
||||
|
||||
% Anders Holtsberg, 18-11-93
|
||||
% Copyright (c) Anders Holtsberg
|
||||
% documentation modified by LeSage to
|
||||
% match the format of the econometrics toolbox
|
||||
|
||||
if (nargin ~= 3)
|
||||
error('Wrong # of arguments to beta_inv');
|
||||
end
|
||||
|
||||
if any(any((a<=0)|(b<=0)))
|
||||
error('beta_inv parameter a or b is nonpositive');
|
||||
end
|
||||
if any(any(abs(2*p-1)>1))
|
||||
error('beta_inv: A probability should be 0<=p<=1');
|
||||
end
|
||||
|
||||
x = a ./ (a+b);
|
||||
dx = 1;
|
||||
while any(any(abs(dx)>256*eps*max(x,1)))
|
||||
dx = (betainc(x,a,b) - p) ./ beta_pdf(x,a,b);
|
||||
x = x - dx;
|
||||
x = x + (dx - x) / 2 .* (x<0);
|
||||
end
|
||||
|
||||
|
|
@ -0,0 +1,44 @@
|
|||
function x = gamm_inv(p,a)
|
||||
% PURPOSE: returns the inverse of the cdf at p of the gamma(a) distribution
|
||||
%---------------------------------------------------
|
||||
% USAGE: x = gamm_inv(p,a)
|
||||
% where: p = a vector of probabilities
|
||||
% a = a scalar parameter gamma(a)
|
||||
%---------------------------------------------------
|
||||
% RETURNS:
|
||||
% a vector x of the quantile at each element of p of the gamma(a) distribution
|
||||
% --------------------------------------------------
|
||||
% SEE ALSO: gamm_d, gamm_pdf, gamm_rnd, gamm_cdf
|
||||
%---------------------------------------------------
|
||||
|
||||
% Anders Holtsberg, 18-11-93
|
||||
% Copyright (c) Anders Holtsberg
|
||||
% documentation modified by LeSage to fit the format
|
||||
% of the econometrics toolbox
|
||||
|
||||
if (nargin ~= 2)
|
||||
error('Wrong # of arguments to gamm_inv');
|
||||
end
|
||||
|
||||
|
||||
if any(any(abs(2*p-1)>1))
|
||||
error('gamm_inv: a probability should be 0<=p<=1')
|
||||
end
|
||||
if any(any(a<=0))
|
||||
error('gamma_inv: parameter a is wrong')
|
||||
end
|
||||
|
||||
x = max(a-1,0.1);
|
||||
dx = 1;
|
||||
while any(any(abs(dx)>256*eps*max(x,1)))
|
||||
dx = (gamm_cdf(x,a) - p) ./ gamm_pdf(x,a);
|
||||
x = x - dx;
|
||||
x = x + (dx - x) / 2 .* (x<0);
|
||||
end
|
||||
|
||||
I0 = find(p==0);
|
||||
x(I0) = zeros(size(I0));
|
||||
I1 = find(p==1);
|
||||
x(I1) = zeros(size(I0)) + Inf;
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,154 @@
|
|||
|
||||
function s = myboxplot (data,notched,symbol,vertical,maxwhisker)
|
||||
|
||||
|
||||
% % % % endif
|
||||
if nargin < 5, maxwhisker = 1; end
|
||||
if nargin < 4, vertical = 1; end
|
||||
if nargin < 3, symbol = ['+','o']; end
|
||||
if nargin < 2, notched = 0; end
|
||||
|
||||
if length(symbol)==1, symbol(2)=symbol(1); end
|
||||
|
||||
if notched==1, notched=0.25; end
|
||||
a=1-notched;
|
||||
|
||||
% ## figure out how many data sets we have
|
||||
if iscell(data),
|
||||
nc = length(data);
|
||||
else
|
||||
% if isvector(data), data = data(:); end
|
||||
nc = size(data,2);
|
||||
end
|
||||
|
||||
% ## compute statistics
|
||||
% ## s will contain
|
||||
% ## 1,5 min and max
|
||||
% ## 2,3,4 1st, 2nd and 3rd quartile
|
||||
% ## 6,7 lower and upper confidence intervals for median
|
||||
s = zeros(7,nc);
|
||||
box = zeros(1,nc);
|
||||
whisker_x = ones(2,1)*[1:nc,1:nc];
|
||||
whisker_y = zeros(2,2*nc);
|
||||
outliers_x = [];
|
||||
outliers_y = [];
|
||||
outliers2_x = [];
|
||||
outliers2_y = [];
|
||||
|
||||
for i=1:nc
|
||||
% ## Get the next data set from the array or cell array
|
||||
if iscell(data)
|
||||
col = data{i}(:);
|
||||
else
|
||||
col = data(:,i);
|
||||
end
|
||||
% ## Skip missing data
|
||||
% % % % % % % col(isnan(col) | isna (col)) = [];
|
||||
col(isnan(col)) = [];
|
||||
|
||||
% ## Remember the data length
|
||||
nd = length(col);
|
||||
box(i) = nd;
|
||||
if (nd > 1)
|
||||
% ## min,max and quartiles
|
||||
% s(1:5,i) = statistics(col)(1:5);
|
||||
s(1,i)=min(col);
|
||||
s(5,i)=max(col);
|
||||
s(2,i)=myprctilecol(col,25);
|
||||
s(3,i)=myprctilecol(col,50);
|
||||
s(4,i)=myprctilecol(col,75);
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
% ## confidence interval for the median
|
||||
est = 1.57*(s(4,i)-s(2,i))/sqrt(nd);
|
||||
s(6,i) = max([s(3,i)-est, s(2,i)]);
|
||||
s(7,i) = min([s(3,i)+est, s(4,i)]);
|
||||
% ## whiskers out to the last point within the desired inter-quartile range
|
||||
IQR = maxwhisker*(s(4,i)-s(2,i));
|
||||
whisker_y(:,i) = [min(col(col >= s(2,i)-IQR)); s(2,i)];
|
||||
whisker_y(:,nc+i) = [max(col(col <= s(4,i)+IQR)); s(4,i)];
|
||||
% ## outliers beyond 1 and 2 inter-quartile ranges
|
||||
outliers = col((col < s(2,i)-IQR & col >= s(2,i)-2*IQR) | (col > s(4,i)+IQR & col <= s(4,i)+2*IQR));
|
||||
outliers2 = col(col < s(2,i)-2*IQR | col > s(4,i)+2*IQR);
|
||||
outliers_x = [outliers_x; i*ones(size(outliers))];
|
||||
outliers_y = [outliers_y; outliers];
|
||||
outliers2_x = [outliers2_x; i*ones(size(outliers2))];
|
||||
outliers2_y = [outliers2_y; outliers2];
|
||||
elseif (nd == 1)
|
||||
% ## all statistics collapse to the value of the point
|
||||
s(:,i) = col;
|
||||
% ## single point data sets are plotted as outliers.
|
||||
outliers_x = [outliers_x; i];
|
||||
outliers_y = [outliers_y; col];
|
||||
else
|
||||
% ## no statistics if no points
|
||||
s(:,i) = NaN;
|
||||
end
|
||||
end
|
||||
% % % % if isempty(outliers2_y)
|
||||
% % % % outliers2_y=
|
||||
% ## Note which boxes don't have enough stats
|
||||
chop = find(box <= 1);
|
||||
|
||||
% ## Draw a box around the quartiles, with width proportional to the number of
|
||||
% ## items in the box. Draw notches if desired.
|
||||
box = box*0.23/max(box);
|
||||
quartile_x = ones(11,1)*[1:nc] + [-a;-1;-1;1;1;a;1;1;-1;-1;-a]*box;
|
||||
quartile_y = s([3,7,4,4,7,3,6,2,2,6,3],:);
|
||||
|
||||
% ## Draw a line through the median
|
||||
median_x = ones(2,1)*[1:nc] + [-a;+a]*box;
|
||||
% median_x=median(col);
|
||||
median_y = s([3,3],:);
|
||||
|
||||
% ## Chop all boxes which don't have enough stats
|
||||
quartile_x(:,chop) = [];
|
||||
quartile_y(:,chop) = [];
|
||||
whisker_x(:,[chop,chop+nc]) = [];
|
||||
whisker_y(:,[chop,chop+nc]) = [];
|
||||
median_x(:,chop) = [];
|
||||
median_y(:,chop) = [];
|
||||
% % % %
|
||||
% ## Add caps to the remaining whiskers
|
||||
cap_x = whisker_x;
|
||||
cap_x(1,:) =cap_x(1,:)- 0.05;
|
||||
cap_x(2,:) =cap_x(2,:)+ 0.05;
|
||||
cap_y = whisker_y([1,1],:);
|
||||
|
||||
% #quartile_x,quartile_y
|
||||
% #whisker_x,whisker_y
|
||||
% #median_x,median_y
|
||||
% #cap_x,cap_y
|
||||
%
|
||||
% ## Do the plot
|
||||
|
||||
mm=min(min(col));
|
||||
MM=max(max(col));
|
||||
|
||||
if vertical
|
||||
plot (quartile_x, quartile_y, 'b', ...
|
||||
whisker_x, whisker_y, 'b--', ...
|
||||
cap_x, cap_y, 'k', ...
|
||||
median_x, median_y, 'r', ...
|
||||
outliers_x, outliers_y, [symbol(1),'r'], ...
|
||||
outliers2_x, outliers2_y, [symbol(2),'r']);
|
||||
set(gca,'XTick',1:nc);
|
||||
set(gca, 'XLim', [0.5, nc+0.5]);
|
||||
set(gca, 'YLim', [mm-(MM-mm)*0.05, MM+(MM-mm)*0.05]);
|
||||
|
||||
else
|
||||
% % % % % plot (quartile_y, quartile_x, "b;;",
|
||||
% % % % % whisker_y, whisker_x, "b;;",
|
||||
% % % % % cap_y, cap_x, "b;;",
|
||||
% % % % % median_y, median_x, "r;;",
|
||||
% % % % % outliers_y, outliers_x, [symbol(1),"r;;"],
|
||||
% % % % % outliers2_y, outliers2_x, [symbol(2),"r;;"]);
|
||||
end
|
||||
|
||||
% % % endfunction
|
||||
|
|
@ -0,0 +1,20 @@
|
|||
function y = myprctilecol(x,p);
|
||||
xx = sort(x);
|
||||
[m,n] = size(x);
|
||||
|
||||
if m==1 | n==1
|
||||
m = max(m,n);
|
||||
if m == 1,
|
||||
y = x*ones(length(p),1);
|
||||
return;
|
||||
end
|
||||
n = 1;
|
||||
q = 100*(0.5:m - 0.5)./m;
|
||||
xx = [min(x); xx(:); max(x)];
|
||||
else
|
||||
q = 100*(0.5:m - 0.5)./m;
|
||||
xx = [min(x); xx; max(x)];
|
||||
end
|
||||
|
||||
q = [0 q 100];
|
||||
y = interp1(q,xx,p);
|
||||
|
|
@ -0,0 +1,44 @@
|
|||
function invp = norm_inv(x, m, sd)
|
||||
% PURPOSE: computes the quantile (inverse of the CDF)
|
||||
% for each component of x with mean m, standard deviation sd
|
||||
%---------------------------------------------------
|
||||
% USAGE: invp = norm_inv(x,m,v)
|
||||
% where: x = variable vector (nx1)
|
||||
% m = mean vector (default=0)
|
||||
% sd = standard deviation vector (default=1)
|
||||
%---------------------------------------------------
|
||||
% RETURNS: invp (nx1) vector
|
||||
%---------------------------------------------------
|
||||
% SEE ALSO: norm_d, norm_rnd, norm_inv, norm_cdf
|
||||
%---------------------------------------------------
|
||||
|
||||
% Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Oct 26, 1994
|
||||
% Copyright Dept of Probability Theory and Statistics TU Wien
|
||||
% Converted to MATLAB by JP LeSage, jpl@jpl.econ.utoledo.edu
|
||||
|
||||
if nargin > 3
|
||||
error ('Wrong # of arguments to norm_inv');
|
||||
end
|
||||
|
||||
[r, c] = size (x);
|
||||
s = r * c;
|
||||
|
||||
if (nargin == 1)
|
||||
m = zeros(1,s);
|
||||
sd = ones(1,s);
|
||||
end
|
||||
|
||||
|
||||
|
||||
x = reshape(x,1,s);
|
||||
m = reshape(m,1,s);
|
||||
sd = reshape(sd,1,s);
|
||||
|
||||
invp = zeros (1,s);
|
||||
|
||||
invp = m + sd .* (sqrt(2) * erfinv(2 * x - 1));
|
||||
|
||||
|
||||
|
||||
invp = reshape (invp, r, c);
|
||||
|
||||
Loading…
Reference in New Issue