univgui.m

常用ROBUST STATISTICAL

            [ahat, bhat] = unifit(ud.X(:,ti));
            randmat(:,i) = unifrnd(ahat,bhat,n,1);
        end        
        hf = figure;
        set(hf,'numbertitle','off','name','EDA: QQ Plot - Uniform Distribution')
        % Upon figure close, this should delete from the array.
        set(hf,'CloseRequestFcn',...
            'tg = findobj(''tag'',''univgui''); H = get(tg,''userdata''); H.plots(find(H.plots == gcf)) = []; set(tg,''userdata'',H); delete(gcf)')
        H.plots = [H.plots, hf];
        set(tg,'userdata',H)
    case 7
        % Poisson distribution
        parmhat = poissfit(ud.X(:,dim));
        for i = 1:pp
            ti = dim(i);
            randmat(:,i) = poissrnd(parmhat(i),n,1);
        end
        hf = figure;
        set(hf,'numbertitle','off','name','EDA: QQ Plot - Poisson Distribution')
        % Upon figure close, this should delete from the array.
        set(hf,'CloseRequestFcn',...
            'tg = findobj(''tag'',''univgui''); H = get(tg,''userdata''); H.plots(find(H.plots == gcf)) = []; set(tg,''userdata'',H); delete(gcf)')
        H.plots = [H.plots, hf];
        set(tg,'userdata',H)
end
% Now that we have the distributions - do the qq plots.
for i = 1:pp
    subplot(nr,nc,i);
    ti = dim(i);
    qqplot(ud.X(:,ti),randmat(:,i));
    % get rid of the x and y axis labels
    xlabel(' '), ylabel(' ')
    if isnumeric(ud.varlab)
        title(num2str(ud.varlab(ti)));
    elseif iscell(ud.varlab)
        title(ud.varlab{ti});
    end
end    


%%%%%%%%%%%%%%%%%%%   HELPER FUNCTIONS %%%%%%%%%%%%%%%%%%%%%%%%%
function H = boxprct(X,vw)
%  BOXPRCT Box-Percentile Plot
%   

% Figure out what inputs were given.
if nargin == 1
    % Then just the plain boxplots.
    w = 2;
    varw = 0;
elseif strcmp(vw,'vw')
    varw = 1;
else
    error('Second argument can only be ''vw''')
end
if isnumeric(X) & varw==1
    error('Cannot do variable width boxplots when the sample sizes are the same.')
end
% Find the type of input argument x
if iscell(X)
    p = length(X);
    for i = 1:p
        % Get all of the sample sizes and the index for the median.
        n(i) = length(X{i});
    end
    if length(unique(n))==1 & varw==1
        error('Cannot do variable width boxplots when the sample sizes are the same.')
    end
    % Get the min and the max of the root n.
    maxn = max(sqrt(n));
    minn = min(sqrt(n));
else
    [n,p] = size(X);
end
% Find some of the things needed to plot them.
% Each boxplot will have a maximum width of 2 units.
ctr = 4:4:4*p;
H = figure; hold on
for i = 1:p
    % Extract the data.
    if iscell(X)
        x = sort(X{i});
    else
        x = sort(X(:,i));
        plain = 1;
    end
    if varw
        % Then it is a variable width boxplot.
        ns = sqrt(n(i));
        % Scale between 0.5 and 2.
        w = scale(ns, minn, maxn, 0.5, 2);
    end
    % Get the quartiles.
    q = quartiles(x);
    if length(n) ~= 1
        % Then we have different sample sizes.
        N = n(i);
    else
        N = n;
    end
    qinds = getk(N);
    KK = qinds(2);
    Lside = [ctr(i)-(1:KK)*w/(N+1), ctr(i)-(N+1-(KK+1:N))*w/(N+1)];
    Rside = [ctr(i)+(1:KK)*w/(N+1), ctr(i)+(N+1-(KK+1:N))*w/(N+1)];
    plot(Lside,x,'k', Rside,x,'k')
    % plot the quartiles.
    k1 = qinds(1);
    plot([ctr(i)-k1*w/(N+1),ctr(i)+k1*w/(N+1)],[q(1),q(1)],'k')
    k2 = qinds(2);
    plot([ctr(i)-k2*w/(N+1),ctr(i)+k2*w/(N+1)],[q(2),q(2)],'k')
    k3 = qinds(3);
    plot([ctr(i)-(N+1-(k3))*w/(N+1),ctr(i)+(N+1-(k3))*w/(N+1)],[q(3),q(3)],'k')
end
ax = axis;
axis([ctr(1)-2 ctr(end)+2 ax(3:4)])
set(gca,'XTickLabel',' ')
hold off

function K = getk(n)
% This function gets the index K for the median of a sample of size n.
K = zeros(1,3);
if rem(n,2)==0
    % then its even
    K(2) = n/2;
    ptrs = K(2)+1:n; 
else
    K(2) = (n+1)/2;
    ptrs = K(2):n;
end
% Get the index to the first quartile.
if rem(K(2),2) ==0
    % If the median is at an even index, then the halves will have even
    % sizes.
    K(1) = K(2)/2;
    K(3) = ptrs(K(1));
else
    K(1) = (K(2)+1)/2;
    K(3) = ptrs(K(1));
end

function H = boxp(X, vw)

%  BOXP  Boxplot and variations
%
%   EXAMPLES
%           boxp(x)      % Plots the plain boxplot for each col/cell of x
%           boxp(x,'vw') % Plots the variable width boxplots
%           boxp(x,'hp') % Plots the histplot

% Figure out what inputs were given.
if nargin == 1
    % Then just the plain boxplots.
    plain = 1;
    varw = 0;
    histp = 0;
    w = 2;
elseif strcmp(vw,'vw')
    plain = 0;
    varw = 1;
    histp = 0;
elseif strcmp(vw,'hp')
    plain = 0;
    varw = 0;
    histp = 1;
else
    error('Second argument can only be ''vw'' or ''hp''')
end
% Find the type of input argument x
if iscell(X)
    p = length(X);
    for i = 1:p
        % Get all of the sample sizes.
        n(i) = length(X{i});
    end
    % Get the min and the max of the root n.
    maxn = max(sqrt(n));
    minn = min(sqrt(n));
else
    [n,p] = size(X);
end
% Find some of the things needed to plot them.
% Each boxplot will have a maximum width of 2 units.
N = 4 + p-1 + 3*p;
ctr = 4:4:4*p;
Le = ctr - 1;
Re = ctr + 1;
H = figure; hold on
for i = 1:p
    % Extract the data.
    if iscell(X)
        x = X{i};
    else
        x = X(:,i);
    end
    % Get the quartiles.
    q = quartiles(x);
    % Get the outliers and adjacent values.
    [adv,outs] = getout(q,x);
    % Get the maximum width of the boxplot.
    if varw
        if length(n) == 1
            error('All sample sizes are equal. Do not do variable width boxplot.')
        end
        % Then it is a variable width boxplot.
        ns = sqrt(n(i));
        % Scale between 0.5 and 2.
        w = scale(ns, minn, maxn, 0.5, 2);
    end
    % Draw the boxes.
    if plain
        % If we have plain boxplots.
        % Draw the quartiles.
        plot([Le(i) Re(i)],[q(1),q(1)])
        plot([Le(i) Re(i)],[q(2),q(2)])
        plot([Le(i) Re(i)],[q(3),q(3)])
        % Draw the sides of the box
        plot([Le(i) Le(i)],[q(1),q(3)])
        plot([Re(i) Re(i)],[q(1),q(3)])
        % Draw the whiskers.
        plot([ctr(i) ctr(i)],[q(1),adv(1)], [ctr(i)-.25 ctr(i)+.25], [adv(1) adv(1)]) 
        plot([ctr(i) ctr(i)],[q(3),adv(2)], [ctr(i)-.25 ctr(i)+.25], [adv(2) adv(2)])
    elseif varw
        % If we have variable width boxplots.
        % Draw the quartiles.
        plot([ctr(i)-w ctr(i)+w],[q(1),q(1)])
        plot([ctr(i)-w ctr(i)+w],[q(2),q(2)])
        plot([ctr(i)-w ctr(i)+w],[q(3),q(3)])
        % Draw the sides of the box
        plot([ctr(i)-w ctr(i)-w],[q(1),q(3)])
        plot([ctr(i)+w ctr(i)+w],[q(1),q(3)])
        % Draw the whiskers.
        ww = scale(ns, minn, maxn, 0.1, 0.25);
        plot([ctr(i) ctr(i)],[q(1),adv(1)], [ctr(i)-ww ctr(i)+ww], [adv(1) adv(1)]) 
        plot([ctr(i) ctr(i)],[q(3),adv(2)], [ctr(i)-ww ctr(i)+ww], [adv(2) adv(2)])
    else
        % We must have the histplot. Plot widths at quartiles proportional
        % to the density estimate there.
        % Draw the quartiles - first get estimates of the density at each.
        for j = 1:3
            fhat(j) = cskern1d(x,q(j));
        end
        w1 = scale(fhat(1),min(fhat),max(fhat),0.5,2);
        w1 = w1/2;
        plot([ctr(i)-w1 ctr(i)+w1],[q(1),q(1)])
        w2 = scale(fhat(2),min(fhat),max(fhat),0.5,2);
        w2 = w2/2;
        plot([ctr(i)-w2 ctr(i)+w2],[q(2),q(2)])
        w3 = scale(fhat(3),min(fhat),max(fhat),0.5,2);
        w3 = w3/2;
        plot([ctr(i)-w3 ctr(i)+w3],[q(3),q(3)])
        % Plot the sides.
        plot([ctr(i)-w1 ctr(i)-w2],[q(1),q(2)])
        plot([ctr(i)+w1 ctr(i)+w2],[q(1),q(2)])
        plot([ctr(i)-w2 ctr(i)-w3],[q(2),q(3)])
        plot([ctr(i)+w2 ctr(i)+w3],[q(2),q(3)])
        % Draw the whiskers.
        plot([ctr(i) ctr(i)],[q(1),adv(1)], [ctr(i)-.25 ctr(i)+.25], [adv(1) adv(1)]) 
        plot([ctr(i) ctr(i)],[q(3),adv(2)], [ctr(i)-.25 ctr(i)+.25], [adv(2) adv(2)])
    end
    % Draw the outliers.
    plot(ctr(i)*ones(size(outs)), outs,'o')
    
    
end
ax = axis;
axis([Le(1)-2 Re(end)+2 ax(3:4)])
set(gca,'XTickLabel',' ')
hold off


function [adv,outs] = getout(q,x)
% This helper function returns the adjacent values and
% outliers.
x = sort(x);
n = length(x);
% Get the upper and lower limits.
iq = q(3) - q(1);
UL = q(3) + iq*1.5;
LL = q(2) - iq*1.5;
% Find any outliers.
ind = [find(x > UL); find(x < LL)];
outs = x(ind);
% Get the adjacent values. Find the
% points that are NOT outliers.
inds = setdiff(1:n,ind);
% Get their min and max.
adv = [x(inds(1)) x(inds(end))];

function nx = scale(x, a, b, c, d)
% This function converts a value x that is orignally between a and b to
% one that is between c and d.
nx = (d - c)*(x - a)/(b - a) + c;

function fhat = cskern1d(data,x)
n = length(data);
fhat = zeros(size(x));
h = 1.06*n^(-1/5);
for i=1:n
   f=exp(-(1/(2*h^2))*(x-data(i)).^2)/sqrt(2*pi)/h;
   fhat = fhat+f/(n);
end

function q = quartiles(x)

%   QUARTILES   Finds the three sample quartiles
%
%   Q = quartiles(X)
%   This returns the three sample quartiles as defined by Tukey,
%   Exploratory Data Analysis, 1977.

% First sort the data.
x = sort(x);
% Get the median.
q2 = median(x);
% First find out if n is even or odd.
n = length(x);
if rem(n,2) == 1
    odd = 1;
else
    odd = 0;
end
if odd
    q1 = median(x(1:(n+1)/2));
    q3 = median(x((n+1)/2:end));
else
    q1 = median(x(1:n/2));
    q3 = median(x(n/2:end));
end
q(1) = q1;
q(2) = q2;
q(3) = q3;