https://github.com/bastibe/Violinplot-Matlab
Tip revision: d1200471a1c80f406749ebe29b5211c2fba34960 authored by Miguel García García on 16 June 2022, 12:26:15 UTC
Fix #48 (#50)
Fix #48 (#50)
Tip revision: d120047
Violin.m
classdef Violin < handle
% Violin creates violin plots for some data
% A violin plot is an easy to read substitute for a box plot
% that replaces the box shape with a kernel density estimate of
% the data, and optionally overlays the data points itself.
% It is also possible to provide two sets of data which are supposed
% to be compared by plotting each column of the two datasets together
% on each side of the violin.
%
% Additional constructor parameters include the width of the
% plot, the bandwidth of the kernel density estimation, the
% X-axis position of the violin plot, and the categories.
%
% Use <a href="matlab:help('violinplot')">violinplot</a> for a
% <a href="matlab:help('boxplot')">boxplot</a>-like wrapper for
% interactive plotting.
%
% See for more information on Violin Plots:
% J. L. Hintze and R. D. Nelson, "Violin plots: a box
% plot-density trace synergism," The American Statistician, vol.
% 52, no. 2, pp. 181-184, 1998.
%
% Violin Properties:
% ViolinColor - Fill color of the violin area and data points.
% Can be either a matrix nx3 or an array of up to two
% cells containing nx3 matrices.
% Defaults to the next default color cycle.
% ViolinAlpha - Transparency of the violin area and data points.
% Can be either a single scalar value or an array of
% up to two cells containing scalar values.
% Defaults to 0.3.
% EdgeColor - Color of the violin area outline.
% Defaults to [0.5 0.5 0.5]
% BoxColor - Color of the box, whiskers, and the outlines of
% the median point and the notch indicators.
% Defaults to [0.5 0.5 0.5]
% MedianColor - Fill color of the median and notch indicators.
% Defaults to [1 1 1]
% ShowData - Whether to show data points.
% Defaults to true
% ShowNotches - Whether to show notch indicators.
% Defaults to false
% ShowMean - Whether to show mean indicator.
% Defaults to false
% ShowBox - Whether to show the box.
% Defaults to true
% ShowMedian - Whether to show the median indicator.
% Defaults to true
% ShowWhiskers - Whether to show the whiskers
% Defaults to true
% HalfViolin - Whether to do a half violin(left, right side) or
% full. Defaults to full.
% QuartileStyle - Option on how to display quartiles, with a
% boxplot, shadow or none. Defaults to boxplot.
% DataStyle - Defines the style to show the data points. Opts:
% 'scatter', 'histogram' or 'none'. Default is 'scatter'.
%
%
% Violin Children:
% ScatterPlot - <a href="matlab:help('scatter')">scatter</a> plot of the data points
% ScatterPlot2 - <a href="matlab:help('scatter')">scatter</a> second plot of the data points
% ViolinPlot - <a href="matlab:help('fill')">fill</a> plot of the kernel density estimate
% ViolinPlot2 - <a href="matlab:help('fill')">fill</a> second plot of the kernel density estimate
% BoxPlot - <a href="matlab:help('fill')">fill</a> plot of the box between the quartiles
% WhiskerPlot - line <a href="matlab:help('plot')">plot</a> between the whisker ends
% MedianPlot - <a href="matlab:help('scatter')">scatter</a> plot of the median (one point)
% NotchPlots - <a href="matlab:help('scatter')">scatter</a> plots for the notch indicators
% MeanPlot - line <a href="matlab:help('plot')">plot</a> at mean value
% Copyright (c) 2016, Bastian Bechtold
% This code is released under the terms of the BSD 3-clause license
properties (Access=public)
ScatterPlot % scatter plot of the data points
ScatterPlot2 % comparison scatter plot of the data points
ViolinPlot % fill plot of the kernel density estimate
ViolinPlot2 % comparison fill plot of the kernel density estimate
BoxPlot % fill plot of the box between the quartiles
WhiskerPlot % line plot between the whisker ends
MedianPlot % scatter plot of the median (one point)
NotchPlots % scatter plots for the notch indicators
MeanPlot % line plot of the mean (horizontal line)
HistogramPlot % histogram of the data
ViolinPlotQ % fill plot of the Quartiles as shadow
end
properties (Dependent=true)
ViolinColor % fill color of the violin area and data points
ViolinAlpha % transparency of the violin area and data points
MarkerSize % marker size for the median dot
LineWidth % linewidth of the median plot
EdgeColor % color of the violin area outline
BoxColor % color of box, whiskers, and median/notch edges
BoxWidth % width of box between the quartiles in axis space (default 10% of Violin plot width, 0.03)
MedianColor % fill color of median and notches
ShowData % whether to show data points
ShowNotches % whether to show notch indicators
ShowMean % whether to show mean indicator
ShowBox % whether to show the box
ShowMedian % whether to show the median line
ShowWhiskers % whether to show the whiskers
HalfViolin % whether to do a half violin(left, right side) or full
end
methods
function obj = Violin(data, pos, varargin)
%Violin plots a violin plot of some data at pos
% VIOLIN(DATA, POS) plots a violin at x-position POS for
% a vector of DATA points.
%
% VIOLIN(..., 'PARAM1', val1, 'PARAM2', val2, ...)
% specifies optional name/value pairs:
% 'Width' Width of the violin in axis space.
% Defaults to 0.3
% 'Bandwidth' Bandwidth of the kernel density
% estimate. Should be between 10% and
% 40% of the data range.
% 'ViolinColor' Fill color of the violin area
% and data points.Can be either a matrix
% nx3 or an array of up to two cells
% containing nx3 matrices.
% 'ViolinAlpha' Transparency of the violin area and data
% points. Can be either a single scalar
% value or an array of up to two cells
% containing scalar values. Defaults to 0.3.
% 'EdgeColor' Color of the violin area outline.
% Defaults to [0.5 0.5 0.5]
% 'BoxColor' Color of the box, whiskers, and the
% outlines of the median point and the
% notch indicators. Defaults to
% [0.5 0.5 0.5]
% 'MedianColor' Fill color of the median and notch
% indicators. Defaults to [1 1 1]
% 'ShowData' Whether to show data points.
% Defaults to true
% 'ShowNotches' Whether to show notch indicators.
% Defaults to false
% 'ShowMean' Whether to show mean indicator.
% Defaults to false
% 'ShowBox' Whether to show the box
% Defaults to true
% 'ShowMedian' Whether to show the median line
% Defaults to true
% 'ShowWhiskers' Whether to show the whiskers
% Defaults to true
% 'HalfViolin' Whether to do a half violin(left, right side) or
% full. Defaults to full.
% 'QuartileStyle' Option on how to display quartiles, with a
% boxplot or as a shadow. Defaults to boxplot.
% 'DataStyle' Defines the style to show the data points. Opts:
% 'scatter', 'histogram' or 'none'. Default is 'Scatter'.
st = dbstack; % get the calling function for reporting errors
namefun = st.name;
args = obj.checkInputs(data, pos, varargin{:});
if length(data)==1
data2 = [];
data = data{1};
else
data2 = data{2};
data = data{1};
end
if isempty(args.ViolinColor)
Release= strsplit(version('-release'), {'a','b'}); %Check release
if str2num(Release{1})> 2019 || strcmp(version('-release'), '2019b')
C = colororder;
else
C = lines;
end
if pos > length(C)
C = lines;
end
args.ViolinColor = {repmat(C,ceil(size(data,2)/length(C)),1)};
end
data = data(not(isnan(data)));
data2 = data2(not(isnan(data2)));
if numel(data) == 1
obj.MedianPlot = scatter(pos, data, 'filled');
obj.MedianColor = args.MedianColor;
obj.MedianPlot.MarkerEdgeColor = args.EdgeColor;
return
end
hold('on');
%% Calculate kernel density estimation for the violin
[density, value, width] = obj.calcKernelDensity(data, args.Bandwidth, args.Width);
% also calculate the kernel density of the comparison data if
% provided
if ~isempty(data2)
[densityC, valueC, widthC] = obj.calcKernelDensity(data2, args.Bandwidth, args.Width);
end
%% Plot the data points within the violin area
if length(density) > 1
[~, unique_idx] = unique(value);
jitterstrength = interp1(value(unique_idx), density(unique_idx)*width, data, 'linear','extrap');
else % all data is identical:
jitterstrength = density*width;
end
if isempty(data2) % if no comparison data
jitter = 2*(rand(size(data))-0.5); % both sides
else
jitter = rand(size(data)); % only right side
end
switch args.HalfViolin % this is more modular
case 'left'
jitter = -1*(rand(size(data))); %left
case 'right'
jitter = 1*(rand(size(data))); %right
case 'full'
jitter = 2*(rand(size(data))-0.5);
end
% Make scatter plot
switch args.DataStyle
case 'scatter'
if ~isempty(data2)
jitter = 1*(rand(size(data))); %right
obj.ScatterPlot = ...
scatter(pos + jitter.*jitterstrength, data, 'filled');
% plot the data points within the violin area
if length(densityC) > 1
jitterstrength = interp1(valueC, densityC*widthC, data2);
else % all data is identical:
jitterstrength = densityC*widthC;
end
jitter = -1*rand(size(data2));% left
obj.ScatterPlot2 = ...
scatter(pos + jitter.*jitterstrength, data2, args.MarkerSize,'filled');
else
obj.ScatterPlot = ...
scatter(pos + jitter.*jitterstrength, data, args.MarkerSize, 'filled');
end
case 'histogram'
[counts,edges] = histcounts(data, size(unique(data),1));
switch args.HalfViolin
case 'right'
obj.HistogramPlot= plot([pos-((counts')/max(counts))*max(jitterstrength)*2, pos*ones(size(counts,2),1)]',...
[edges(1:end-1)+max(diff(edges))/2; edges(1:end-1)+max(diff(edges))/2],'-','LineWidth',1, 'Color', 'k');
case 'left'
obj.HistogramPlot= plot([pos*ones(size(counts,2),1), pos+((counts')/max(counts))*max(jitterstrength)*2]',...
[edges(1:end-1)+max(diff(edges))/2; edges(1:end-1)+max(diff(edges))/2],'-','LineWidth',1, 'Color', 'k');
otherwise
fprintf([namefun, ' No histogram/bar plot option available for full violins, as it would look overcrowded.\n'])
end
case 'none'
end
%% Plot the violin
halfViol= ones(1, size(density,2));
if isempty(data2) % if no comparison data
switch args.HalfViolin
case 'right'
obj.ViolinPlot = ... % plot color will be overwritten later
fill([pos+density*width halfViol*pos], ...
[value value(end:-1:1)], [1 1 1]);
case 'left'
obj.ViolinPlot = ... % plot color will be overwritten later
fill([halfViol*pos pos-density(end:-1:1)*width], ...
[value value(end:-1:1)], [1 1 1]);
case 'full'
obj.ViolinPlot = ... % plot color will be overwritten later
fill([pos+density*width pos-density(end:-1:1)*width], ...
[value value(end:-1:1)], [1 1 1]);
end
else
% plot right half of the violin
obj.ViolinPlot = ...
fill([pos+density*width pos-density(1)*width], ...
[value value(1)], [1 1 1]);
% plot left half of the violin
obj.ViolinPlot2 = ...
fill([pos-densityC(end)*widthC pos-densityC(end:-1:1)*widthC], ...
[valueC(end) valueC(end:-1:1)], [1 1 1]);
end
%% Plot the quartiles within the violin
quartiles = quantile(data, [0.25, 0.5, 0.75]);
flat= [halfViol*pos halfViol*pos];
switch args.QuartileStyle
case 'shadow'
switch args.HalfViolin
case 'right'
w = [pos+density*width halfViol*pos];
h= [value value(end:-1:1)];
case 'left'
w = [halfViol*pos pos-density(end:-1:1)*width];
h= [value value(end:-1:1)];
case 'full'
w = [pos+density*width pos-density(end:-1:1)*width];
h= [value value(end:-1:1)];
end
w(h<quartiles(1))=flat(h<quartiles(1));
w(h>quartiles(3))=flat((h>quartiles(3)));
obj.ViolinPlotQ = ... % plot color will be overwritten later
fill(w, ...
h, [1 1 1]);
case 'boxplot'
obj.BoxPlot = ... % plot color will be overwritten later
fill(pos+[-1,1,1,-1]*args.BoxWidth, ...
[quartiles(1) quartiles(1) quartiles(3) quartiles(3)], ...
[1 1 1]);
case 'none'
end
%% Plot the data mean
meanValue = mean(data);
if length(density) > 1
[~, unique_idx] = unique(value);
meanDensityWidth = interp1(value(unique_idx), density(unique_idx), meanValue, 'linear','extrap')*width;
else % all data is identical:
meanDensityWidth = density*width;
end
if meanDensityWidth<args.BoxWidth/2
meanDensityWidth=args.BoxWidth/2;
end
switch args.HalfViolin
case 'right'
obj.MeanPlot = plot(pos+[0,1].*meanDensityWidth, ...
[meanValue, meanValue]);
case 'left'
obj.MeanPlot = plot(pos+[-1,0].*meanDensityWidth, ...
[meanValue, meanValue]);
case 'full'
obj.MeanPlot = plot(pos+[-1,1].*meanDensityWidth, ...
[meanValue, meanValue]);
end
obj.MeanPlot.LineWidth = 1;
%% Plot the median, notch, and whiskers
IQR = quartiles(3) - quartiles(1);
lowhisker = quartiles(1) - 1.5*IQR;
lowhisker = max(lowhisker, min(data(data > lowhisker)));
hiwhisker = quartiles(3) + 1.5*IQR;
hiwhisker = min(hiwhisker, max(data(data < hiwhisker)));
if ~isempty(lowhisker) && ~isempty(hiwhisker)
obj.WhiskerPlot = plot([pos pos], [lowhisker hiwhisker]);
end
% Median
obj.MedianPlot = scatter(pos, quartiles(2), args.MarkerSize, [1 1 1], 'filled');
% Notches
obj.NotchPlots = ...
scatter(pos, quartiles(2)-1.57*IQR/sqrt(length(data)), ...
[], [1 1 1], 'filled', '^');
obj.NotchPlots(2) = ...
scatter(pos, quartiles(2)+1.57*IQR/sqrt(length(data)), ...
[], [1 1 1], 'filled', 'v');
%% Set graphical preferences
obj.EdgeColor = args.EdgeColor;
obj.MedianPlot.LineWidth = args.LineWidth;
obj.BoxColor = args.BoxColor;
obj.BoxWidth = args.BoxWidth;
obj.MedianColor = args.MedianColor;
obj.ShowData = args.ShowData;
obj.ShowNotches = args.ShowNotches;
obj.ShowMean = args.ShowMean;
obj.ShowBox = args.ShowBox;
obj.ShowMedian = args.ShowMedian;
obj.ShowWhiskers = args.ShowWhiskers;
if not(isempty(args.ViolinColor))
if size(args.ViolinColor{1},1) > 1
ViolinColor{1} = args.ViolinColor{1}(pos,:);
else
ViolinColor{1} = args.ViolinColor{1};
end
if length(args.ViolinColor)==2
if size(args.ViolinColor{2},1) > 1
ViolinColor{2} = args.ViolinColor{2}(pos,:);
else
ViolinColor{2} = args.ViolinColor{2};
end
else
ViolinColor{2} = ViolinColor{1};
end
else
% defaults
if args.scpltBool
ViolinColor{1} = obj.ScatterPlot.CData;
else
ViolinColor{1} = [0 0 0];
end
ViolinColor{2} = [0 0 0];
end
obj.ViolinColor = ViolinColor;
if not(isempty(args.ViolinAlpha))
if length(args.ViolinAlpha{1})>1
error('Only scalar values are accepted for the alpha color channel');
else
ViolinAlpha{1} = args.ViolinAlpha{1};
end
if length(args.ViolinAlpha)==2
if length(args.ViolinAlpha{2})>1
error('Only scalar values are accepted for the alpha color channel');
else
ViolinAlpha{2} = args.ViolinAlpha{2};
end
else
ViolinAlpha{2} = ViolinAlpha{1}/2; % default unless specified
end
else
% default
ViolinAlpha = {1,1};
end
obj.ViolinAlpha = ViolinAlpha;
end
%% SET METHODS
function set.EdgeColor(obj, color)
if ~isempty(obj.ViolinPlot)
obj.ViolinPlot.EdgeColor = color;
obj.ViolinPlotQ.EdgeColor = color;
if ~isempty(obj.ViolinPlot2)
obj.ViolinPlot2.EdgeColor = color;
end
end
end
function color = get.EdgeColor(obj)
if ~isempty(obj.ViolinPlot)
color = obj.ViolinPlot.EdgeColor;
end
end
function set.MedianColor(obj, color)
obj.MedianPlot.MarkerFaceColor = color;
if ~isempty(obj.NotchPlots)
obj.NotchPlots(1).MarkerFaceColor = color;
obj.NotchPlots(2).MarkerFaceColor = color;
end
end
function color = get.MedianColor(obj)
color = obj.MedianPlot.MarkerFaceColor;
end
function set.BoxColor(obj, color)
if ~isempty(obj.BoxPlot)
obj.BoxPlot.FaceColor = color;
obj.BoxPlot.EdgeColor = color;
obj.WhiskerPlot.Color = color;
obj.MedianPlot.MarkerEdgeColor = color;
obj.NotchPlots(1).MarkerFaceColor = color;
obj.NotchPlots(2).MarkerFaceColor = color;
elseif ~isempty(obj.ViolinPlotQ)
obj.WhiskerPlot.Color = color;
obj.MedianPlot.MarkerEdgeColor = color;
obj.NotchPlots(1).MarkerFaceColor = color;
obj.NotchPlots(2).MarkerFaceColor = color;
end
end
function color = get.BoxColor(obj)
if ~isempty(obj.BoxPlot)
color = obj.BoxPlot.FaceColor;
end
end
function set.BoxWidth(obj,width)
if ~isempty(obj.BoxPlot)
pos=mean(obj.BoxPlot.XData);
obj.BoxPlot.XData=pos+[-1,1,1,-1]*width;
end
end
function width = get.BoxWidth(obj)
width=max(obj.BoxPlot.XData)-min(obj.BoxPlot.XData);
end
function set.ViolinColor(obj, color)
obj.ViolinPlot.FaceColor = color{1};
obj.ScatterPlot.MarkerFaceColor = color{1};
obj.MeanPlot.Color = color{1};
if ~isempty(obj.ViolinPlot2)
obj.ViolinPlot2.FaceColor = color{2};
obj.ScatterPlot2.MarkerFaceColor = color{2};
end
if ~isempty(obj.ViolinPlotQ)
obj.ViolinPlotQ.FaceColor = color{1};
end
for idx = 1: size(obj.HistogramPlot,1)
obj.HistogramPlot(idx).Color = color{1};
end
end
function color = get.ViolinColor(obj)
color{1} = obj.ViolinPlot.FaceColor;
if ~isempty(obj.ViolinPlot2)
color{2} = obj.ViolinPlot2.FaceColor;
end
end
function set.ViolinAlpha(obj, alpha)
obj.ViolinPlotQ.FaceAlpha = .65;
obj.ViolinPlot.FaceAlpha = alpha{1};
obj.ScatterPlot.MarkerFaceAlpha = 1;
if ~isempty(obj.ViolinPlot2)
obj.ViolinPlot2.FaceAlpha = alpha{2};
obj.ScatterPlot2.MarkerFaceAlpha = alpha{2};
end
end
function alpha = get.ViolinAlpha(obj)
alpha{1} = obj.ViolinPlot.FaceAlpha;
if ~isempty(obj.ViolinPlot2)
alpha{2} = obj.ViolinPlot2.FaceAlpha;
end
end
function set.ShowData(obj, yesno)
if yesno
obj.ScatterPlot.Visible = 'on';
for idx = 1: size(obj.HistogramPlot,1)
obj.HistogramPlot(idx).Visible = 'on';
end
else
obj.ScatterPlot.Visible = 'off';
for idx = 1: size(obj.HistogramPlot,1)
obj.HistogramPlot(idx).Visible = 'off';
end
end
if ~isempty(obj.ScatterPlot2)
obj.ScatterPlot2.Visible = obj.ScatterPlot.Visible;
end
end
function yesno = get.ShowData(obj)
if ~isempty(obj.ScatterPlot)
yesno = strcmp(obj.ScatterPlot.Visible, 'on');
end
end
function set.ShowNotches(obj, yesno)
if ~isempty(obj.NotchPlots)
if yesno
obj.NotchPlots(1).Visible = 'on';
obj.NotchPlots(2).Visible = 'on';
else
obj.NotchPlots(1).Visible = 'off';
obj.NotchPlots(2).Visible = 'off';
end
end
end
function yesno = get.ShowNotches(obj)
if ~isempty(obj.NotchPlots)
yesno = strcmp(obj.NotchPlots(1).Visible, 'on');
end
end
function set.ShowMean(obj, yesno)
if ~isempty(obj.MeanPlot)
if yesno
obj.MeanPlot.Visible = 'on';
else
obj.MeanPlot.Visible = 'off';
end
end
end
function yesno = get.ShowMean(obj)
if ~isempty(obj.BoxPlot)
yesno = strcmp(obj.BoxPlot.Visible, 'on');
end
end
function set.ShowBox(obj, yesno)
if ~isempty(obj.BoxPlot)
if yesno
obj.BoxPlot.Visible = 'on';
else
obj.BoxPlot.Visible = 'off';
end
end
end
function yesno = get.ShowBox(obj)
if ~isempty(obj.BoxPlot)
yesno = strcmp(obj.BoxPlot.Visible, 'on');
end
end
function set.ShowMedian(obj, yesno)
if ~isempty(obj.MedianPlot)
if yesno
obj.MedianPlot.Visible = 'on';
else
obj.MedianPlot.Visible = 'off';
end
end
end
function yesno = get.ShowMedian(obj)
if ~isempty(obj.MedianPlot)
yesno = strcmp(obj.MedianPlot.Visible, 'on');
end
end
function set.ShowWhiskers(obj, yesno)
if ~isempty(obj.WhiskerPlot)
if yesno
obj.WhiskerPlot.Visible = 'on';
else
obj.WhiskerPlot.Visible = 'off';
end
end
end
function yesno = get.ShowWhiskers(obj)
if ~isempty(obj.WhiskerPlot)
yesno = strcmp(obj.WhiskerPlot.Visible, 'on');
end
end
end
methods (Access=private)
function results = checkInputs(~, data, pos, varargin)
isscalarnumber = @(x) (isnumeric(x) & isscalar(x));
p = inputParser();
p.addRequired('Data', @(x)isnumeric(vertcat(x{:})));
p.addRequired('Pos', isscalarnumber);
p.addParameter('Width', 0.3, isscalarnumber);
p.addParameter('Bandwidth', [], isscalarnumber);
iscolor = @(x) (isnumeric(x) & size(x,2) == 3);
p.addParameter('ViolinColor', [], @(x)iscolor(vertcat(x{:})));
p.addParameter('MarkerSize', 36, @isnumeric);
p.addParameter('LineWidth', 0.75, @isnumeric);
p.addParameter('BoxColor', [0.5 0.5 0.5], iscolor);
p.addParameter('BoxWidth', 0.01, isscalarnumber);
p.addParameter('EdgeColor', [0.5 0.5 0.5], iscolor);
p.addParameter('MedianColor', [1 1 1], iscolor);
p.addParameter('ViolinAlpha', {0.3,0.15}, @(x)isnumeric(vertcat(x{:})));
isscalarlogical = @(x) (islogical(x) & isscalar(x));
p.addParameter('ShowData', true, isscalarlogical);
p.addParameter('ShowNotches', false, isscalarlogical);
p.addParameter('ShowMean', false, isscalarlogical);
p.addParameter('ShowBox', true, isscalarlogical);
p.addParameter('ShowMedian', true, isscalarlogical);
p.addParameter('ShowWhiskers', true, isscalarlogical);
validSides={'full', 'right', 'left'};
checkSide = @(x) any(validatestring(x, validSides));
p.addParameter('HalfViolin', 'full', checkSide);
validQuartileStyles={'boxplot', 'shadow', 'none'};
checkQuartile = @(x)any(validatestring(x, validQuartileStyles));
p.addParameter('QuartileStyle', 'boxplot', checkQuartile);
validDataStyles = {'scatter', 'histogram', 'none'};
checkStyle = @(x)any(validatestring(x, validDataStyles));
p.addParameter('DataStyle', 'scatter', checkStyle);
p.parse(data, pos, varargin{:});
results = p.Results;
end
end
methods (Static)
function [density, value, width] = calcKernelDensity(data, bandwidth, width)
if isempty(data)
error('Empty input data');
end
[density, value] = ksdensity(data, 'bandwidth', bandwidth);
density = density(value >= min(data) & value <= max(data));
value = value(value >= min(data) & value <= max(data));
value(1) = min(data);
value(end) = max(data);
value = [value(1)*(1-1E-5), value, value(end)*(1+1E-5)];
density = [0, density, 0];
% all data is identical
if min(data) == max(data)
density = 1;
value= mean(value);
end
width = width/max(density);
end
end
end
