Can I smoothen the cdfplot() result?

Hi, suppose I have a sample data x, and i want to plot the the cdf function on the data and plot it, ofcoz I will use cdfplot(x) and I will a graph similar to the one attached below. My question is, can i smoothen it, is there a way I can make the plot smoother and not look like a a step funtion plot?

[Alexander]

Can you supply the code and data for your figure above?

[Mathieu NOE]

even with smoothing , I doubt you can make those large flats disappear

wuld probably rather try to fit either a polynomial or exponential function

as mentionned by @Alexander, sharing the data would be a plus...

[Mathieu NOE]

hello

so have you solved your problem ?

Answer 1

[Alexander]

@Mathieu NOE, thank you for remaindig me. This question kept nagging me, but I just forgotten it. Here are some rudimentary aproches. As the data were artificially generated, with real measurement the solution are only a rough direction what can be done and needed a lot of refinement.

commandwindow;

y = [];

for ii = 1:20

y = [y ones(1,100*ii)*ii];

end

n = length(y);

x = linspace(1,25,n);

plot(x,y);grid minor;

% 1. Polynomial Fit, 2. Order

% As suggested by @Mathieu NOE

P = polyfit(x,y,2);

yNew = polyval(P,x);

plot(x,y,x,yNew);grid minor; title('Polyfit');pause;

% 2. Logarithmic Fit

% Could be done much easier with lsqcurvefit, if you have optimization

% toobox. I don't have it, hence, it's handcrafted.

yLogSum = sum(y.*log(x)); ySum = sum(y); sumLnx = sum(log(x));

Zaeler = yLogSum - (ySum*sumLnx)/n;

Nenner = sum(log(x).^2) - sumLnx^2/n;

b = Zaeler/Nenner;

a = ySum/n - b*sumLnx/n;

yLogRes = a + b*log(x);

plot(x,y,x,yLogRes);grid minor; title('Logarithmic Fit');pause;

% 3. Filtering

% I just used a butterworth filter. As more longer the steps are, as higher

% is the weaviness of the output. The filter parameters has to be adjusted

% to fit the input curve. Another filter type might fit better.

[B,A] = butter(2,0.0005);

yFF = filtfilt(B,A,y);

plot(x,y,x,yFF);grid minor; title('Butterworth filtered');pause;

% 4. Step/Edge Evaluation

% The input was generated. With measured data the following lines must be

% adjusted.

% Also the data has to be interpolated between the x-values. yDiff has now

% as many data values as there are steps in the measurement. Maybe with

% interp().

yDiff = diff(y);

yFind = find(yDiff > 0.5);

Edge = 1; %Upper Edge = 1; lower edge = 0;

plot(x,y,x(yFind+Edge),y(yFind+Edge));grid minor; title('Edge Evaluation');

legend('Original',['Edge ' num2str(Edge)], 'location','northwest')

[Mathieu NOE]

hello again

I slightly modified your code (first I reduced the data size by factor 10 to run it faster, but that's really minor) , then I think using the upper and lower edge approach we can build a good solution

would like to continue with real data

y = [];
for ii = 1:20
    y = [y  ones(1,10*ii)*ii];
end
n = length(y);
x = linspace(1,25,n);

% 1. Polynomial Fit, 2. Order
% As suggested by @Mathieu NOE
P = polyfit(x,y,6);
yNew = polyval(P,x);
figure,plot(x,y,x,yNew);grid minor; title('Polyfit');

% 2. Logarithmic Fit
% Could be done much easier with lsqcurvefit, if you have optimization
% toobox. I don't have it, hence, it's handcrafted.
yLogSum    = sum(y.*log(x)); ySum = sum(y); sumLnx = sum(log(x));
Zaeler = yLogSum - (ySum*sumLnx)/n;
Nenner = sum(log(x).^2) - sumLnx^2/n;
b = Zaeler/Nenner;
a = ySum/n - b*sumLnx/n;
yLogRes = a + b*log(x);
figure,plot(x,y,x,yLogRes);grid minor; title('Logarithmic Fit');

% 3. Filtering
% I just used a butterworth filter. As more longer the steps are, as higher
% is the weaviness of the output. The filter parameters has to be adjusted
% to fit the input curve. Another filter type might fit better.
[B,A] = butter(2,0.005);
yFF = filtfilt(B,A,y);
plot(x,y,x,yFF);grid minor; title('Butterworth filtered');

% 4. Step/Edge Evaluation
% The input was generated. With measured data the following lines must be
% adjusted.
% Also the data has to be interpolated between the x-values. yDiff has now
% as many data values as there are steps in the measurement. Maybe with
% interp().
yDiff = diff(y);
yFind = find(yDiff > 0.5);

Edge = 0; %Upper Edge = 1; lower edge = 0;
xl = x(yFind+Edge);
yl = y(yFind+Edge);

Edge = 1; %Upper Edge = 1; lower edge = 0;
xu = x(yFind+Edge);
yu = y(yFind+Edge);

% Median curve 
xm = (xu+xl)/2;
ym = (yu+yl)/2;

figure,plot(x,y,xl,yl,xu,yu,xm,ym);grid minor; title('Edge Evaluation');
legend('Original','Lower Edge','Upper Edge','Median curve' , 'location','northwest')

[Alexander]

Thanks for the refinement, it looks good. And I agree, real data would be fine.

[okoth ochola]

Thank you all fro these wonderful answers. You dont how much have been helped

[Alexander]

Cheers and upvote Mathieu also.

Answer 2

[Mathieu NOE]

some more code to look at - I kept only the best solutions , even though there are probably lot of other solutions

also I introduce some randomness in the data to see how robust is the code

enjoy !

y = [];
for ii = 1:20
    y = [y  ones(1,round(10*ii+25*rand))*ii];
end
n = length(y);
x = linspace(1,25,n);
% add some noise on y 
y = y + 0.05*randn(size(y));

% 1. Polynomial Fit,
% As suggested by @Mathieu NOE
P = polyfit(x,y,6);
yNew = polyval(P,x);
figure,plot(x,y,x,yNew);grid minor; title('Polyfit');

% 2. Step/Edge Evaluation
% The input was generated. With measured data the following lines must be
% adjusted.
% Also the data has to be interpolated between the x-values. yDiff has now
% as many data values as there are steps in the measurement. Maybe with
% interp().

yDiff = diff(y);
yFind = find(yDiff > max(yDiff)/2);

Edge = 0; %Upper Edge = 1; lower edge = 0;
xl = x(yFind+Edge);
yl = y(yFind+Edge);

Edge = 1; %Upper Edge = 1; lower edge = 0;
xu = x(yFind+Edge);
yu = y(yFind+Edge);

% Median curve 
xm = (xu+xl)/2;
ym = (yu+yl)/2;

figure,plot(x,y,xl,yl,xu,yu,xm,ym);grid minor; title('Edge Evaluation');
legend('Original','Lower Edge','Upper Edge','Median curve' , 'location','northwest')

% 3. Envelop (secant method)
N = 200;
[yu] = env_secant(x, y, N, 'top'); % (see function attached)
[yl] = env_secant(x, y, N, 'bottom'); % (see function attached)
% Median curve
ym = (yu+yl)/2;

figure,plot(x,y,x,yl,x,yu,x,ym);grid minor; title('Envelop (secant method)');
legend('Original','Lower Envelop','Upper Envelop','Median curve' , 'location','northwest')

%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [env] = env_secant(x_data, y_data, view, side)
% Function call: env_secant(x_data, y_data, view, side)
% Calculates the top envelope of data <y_data> over <x_data>.
% Method used: 'secant-method'
% env_secant() observates the max. slope of about <view> points,
% and joints them to the resulting envelope.
% An interpolation over original x-values is done finally.
% <side> ('top' or 'bottom') defines which side to evolve.
% Author: Andreas Martin, Volkswagen AG, Germany

    if nargin == 0
        test( false );
        test( true );
        return
    end

    if nargin < 4
        error( '%s needs at least 4 input arguments!', mfilename );
    end

    assert( isnumeric(view) && isscalar(view) && view > 1, ...
           'Parameter <view> must be a value greater than 1!' );
    assert( isvector(x_data) && isnumeric(x_data) && all( isfinite( x_data ) ), ...
           'Parameter <x_data> has to be of vector type, holding finite numeric values!' );
    assert( isvector(y_data) && isnumeric(y_data) && all( isfinite( y_data ) ), ...
           'Parameter <y_data> has to be of vector type, holding finite numeric values!' );
    assert( isequal( size(x_data), size(y_data) ), ...
           'Parameters <x_data> and <y_data> must have same size and dimension!' );
    assert( ischar(side), ...
            'Parameter <side> must be ''top'' or ''bottom''!' );

    switch lower( side )
        case 'top'
            side = 1;
        case 'bottom'
            side = -1;
        otherwise
            error( 'Parameter <side> must be ''top'' or ''bottom''!' );
    end

    sz = size( x_data );

    x_data = x_data(:);
    x_diff = diff( x_data );
    x_diff = [min(x_diff), max(x_diff)];
    assert( x_diff(1) > 0, '<x_data> must be monotonic increasing!' );

    y_data = y_data(:);
    data_len = length( y_data );
    x_new = [];
    y_new = [];

    if diff( x_diff ) < eps( max(x_data) ) + eps( min(x_data) )
        % x_data is equidistant
        search_fcn = @( y_data, ii, i ) ...
                     max( ( y_data(ii) - y_data(i) ) ./ (ii-i)' * side );
    else
        % x_data is not equidistant
        search_fcn = @( y_data, ii, i ) ...
                     max( ( y_data(ii) - y_data(i) ) ./ ( x_data(ii) - x_data(i) ) * side );
    end


    i = 1;
    while i < data_len;
        ii = i+1:min( i + view, data_len );
        [ m, idx ] = search_fcn( y_data, ii, i );

        % New max. slope: store new "observation point"
        i = i + idx;
        x_new(end+1) = x_data(i);
        y_new(end+1) = y_data(i);
    end;

    env = interp1( x_new, y_new, x_data, 'linear', 'extrap' );
    env = reshape( env, sz );




function test( flagMonotonic )
    npts = 100000;
    y_data = cumsum( randn( npts, 1 ) ) .* cos( (1:npts)/50 )' + 100 * cos( (1:npts)/6000 )';
    if flagMonotonic
        x_data = (1:npts)' + 10;
    else
        x_diff = rand( size( y_data ) );
        x_data = cumsum( x_diff );
    end

    view = ceil( npts * 0.01 ); % 1 Percent of total length
    env_up = env_secant( x_data, y_data, view, 'top' );
    env_lo = env_secant( x_data, y_data, view, 'bottom' );

    figure
    plot( x_data, y_data, '-', 'Color', 0.8 * ones(3,1) );
    hold all
    h(1) = plot( x_data, env_up, 'b-', 'DisplayName', 'top' );
    h(2) = plot( x_data, env_lo, 'g-', 'DisplayName', 'bottom' );

    grid
    legend( 'show', h )
end

end