rng(5);
data = chi2rnd(5,1000,1);
mc1 = medcouple(data);
stddev = std(data);
figure;
histogram(data);
title(['Medcouple = ' num2str(mc1) ' -- Standard deviation = ' num2str(stddev)],'FontSize',16);

rng(5);
method = 1;
mc2 = medcouple(chi2rnd(5,1000,1),method);

z = chi2rnd(5,25,1);
mcm = 1; % use naive method (n is small ...)
medoutR = 0.5*(-medcouple(repmat(max(z),size(z,1),1)-z , mcm) + medcouple(z,mcm));
medout  = medcouple(z,mcm);
disp('Application to a non problematic dataset:');
disp(['result with reflected option = ' num2str(medoutR) ' --- result on z only = ' num2str(medoutR)]);

rng(5);
%orientation = 'horizontal'; 
orientation = 'vertical'; 
% normal data with an outlier
datain1 = 15+randn(200,1);     
datain1 = [datain1 ; 35];      
% log-normal data with an outlier
datain2 = lognrnd(0,1,200,1);  
datain2 = [datain2 ; 35]; 
k  = 1.5;
% limits on normal data + one outlier
MC_1  = medcouple(datain1);
Q1_1  = quantile(datain1,0.25);
Q3_1  = quantile(datain1,0.75);
IQR_1 = Q3_1 - Q1_1;
DataLimLower_1 = Q1_1 - k * exp(-3.5 * MC_1) * IQR_1;
DataLimUpper_1 = Q3_1 + k * exp(4 * MC_1)    * IQR_1;
% limits on log-normal data + one outlier
MC_2  = medcouple(datain2);
Q1_2  = quantile(datain2,0.25);
Q3_2  = quantile(datain2,0.75);
IQR_2 = Q3_2 - Q1_2;
DataLimLower_2 = Q1_2 - k * exp(-3.5 * MC_2) * IQR_2;
DataLimUpper_2 = Q3_2 + k * exp(4 * MC_2)    * IQR_2;
% boxplot
boxplot([datain1,datain2],'Whisker',k,'Orientation',orientation,...
'OutlierSize',10,'Symbol','*r', ...
'labels',{'Normal(15,1)+outlier','Log-normal(0,1)+outlier'});
title({'Boxplot adjusted for skewness' , 'with medcouple (red limits)'}, 'FontSize' , 20)
set(gca,'FontSize',16);
switch orientation
case 'vertical'
% add the limits (for vertical boxplot)
xx  = xlim(gca);
xx2 = xx/1.1;   xx2(1) = xx(2)/1.5;
xx1 = xx/2; xx1(1) = xx(2)/3.5;
line(xx1,[DataLimLower_1 , DataLimLower_1] , 'Color','red');
line(xx1,[DataLimUpper_1 , DataLimUpper_1] , 'Color','red');
line(xx2,[DataLimLower_2 , DataLimLower_2] , 'Color','red');
line(xx2,[DataLimUpper_2 , DataLimUpper_2] , 'Color','red');
case 'horizontal'
% add the limits (for horizontal boxplot)
yy  = ylim(gca);
yy2 = yy/1.1;   yy2(1) = yy(2)/1.5;
yy1 = yy/2; yy1(1) = yy(2)/3.5;
line([DataLimLower_1 , DataLimLower_1] , yy1,'Color','red');
line([DataLimUpper_1 , DataLimUpper_1] , yy1,'Color','red');
line([DataLimLower_2 , DataLimLower_2] , yy2,'Color','red');
line([DataLimUpper_2 , DataLimUpper_2] , yy2,'Color','red');
end