Here is how I was able to solve it in the end (if anybody wants to do a similar display):
stats=[3.0000 26.7150 43.0300 53.6300;
4.0000 17.5450 33.4050 80.7700;
5.0000 32.3300 39.6800 109.7800;
15.0000 179.5300 179.5300 179.5300;
16.0000 21.5550 96.1550 156.43003]; % the row size of this matrix changes depending on previous calculations
facnrs=stats(:,1);
P10=stats(:,2);
P50=stats(:,3);
P90=stats(:,4);
neg=P50-P10;
pos=P90-P50;
% plot the results
errorbar(stats(:,1),P50,neg,pos,'o')
xsmin=min(stats(:,1))-1;
xsmax=max(stats(:,1))+1;
set(gca,'xlim',[xsmin xsmax])
set(gca,'ydir','reverse');
xlabel('Facies')
ylabel('Depth (m) from top for P10, P50 & P90')
ylh = get(gca,'ylabel');
ylp = get(ylh, 'Position');
set(ylh, 'Rotation',270, 'Position',ylp, 'VerticalAlignment','middle', 'HorizontalAlignment','center')
% label the number of occurrences for each facies on top of errorbars;
% start with counting the nr of occurrences per facies
for i=1:length(facnrs)
x=find(facies==facnrs(i));
noccur(i)=length(x);
end
noccur=noccur';
% compile string for number of occurrences
for k = 1:length(facnrs)
bartext(k) = compose('$\\stackrel{n = %d}{\\downarrow}$', noccur(k));
end
for m=1:length(facnrs)
textXpos(m)=facnrs(m)-0.5;
end
textXpos=textXpos';
for j=1:length(facnrs)
textYpos(j)=P10(j)-10;
end
textYpos=textYpos';
% label the bars
text(textXpos,textYpos,bartext,'interpreter','latex','fontsize',20);
Plot looks like this: