m=100000;
z=1:100;
p=0.25;
alpha=4;
beta=150;
n1=(1-p)*length(z); # limit for tail threshold, as measured by the observations ordered
z=sort(z);
U=runif(m);
ind=floor(1+U*n1); # m random indices, up to tail threshold, for mixed sampling further down
Y=z[n1]+(U**(-1/alpha)-1)*beta; # with probability p the claim comes from z_(n1) + Pareto(alpha, beta)
L=runif(m)<1-p; # the part of the Monte Carlo sample to be sampled non parametrically
Z=L*z[ind]+(1-L)*Y;
mean(Z)