function [pfss, gRe, cdf, zRe] = SS(Fun, n, paras)
% Algorithm parameters
N = paras.NumSam;
p0 = paras.CondPro;
nt = round(N*p0);
ns = ceil(1/p0-1);
%% Step 1 (Monte Carlo simulation)
z = normrnd(0,1,N,n); % Generate standard normal samples
% Calculate LSF
for i=1:N;
g(i) = feval(Fun,z(i,:));
end;
% Sort samples and LSF values
[gSort,index]=sort(g);
Z = z(index,:);
% Record the calculated results
pfss = nt/N;
gRe = gSort;
cdf = (N-[N:-1:1]+1)/N;
zRe = Z;
%% Step 2 (Markov Chain Monte Carlo)
% preparation for MCMC
sigma = std(Z(1:nt,:),0,1);
stopFlag = 0;
iter = 1;
while (stopFlag == 0)
w = Z(1:nt,:);
g = gSort(1:nt);
iter = iter+1;
len = nt+1;
for i=1:nt
seed = w(i,:);
seed_g = g(i);
for j = 1:ns % A Markov chain
for k = 1:n % Component by component
% Generate an candidate
u(k) = seed(k)+(2*rand()-1)*sigma(k);
% Calculate the acceptance probability
pdf2 = exp(-0.5*seed(k).^2);
pdf1 = exp(-0.5*u(k).^2);
alpha = pdf1/pdf2;
% Accept u(k) with a probability of alpha
if alpha > rand();
w(len,k)= u(k);
else
w(len,k)= seed(k);
end;
end
gTemp = feval(Fun,w(len,:));
% Accept or reject
if gTemp <= gSort(nt+1);
g(len) = gTemp;
seed = w(len,:);
seed_g = g(len);
else
g(len) = seed_g;
w(len,:) = seed;
end
len = len+1;
% terminate the simulation of a Markov chain
if len > N break; end
end
if len > N break; end
end
% Sort samples and LSF values
[gSort,index]=sort(g);
Z = w(index,:);
% Record the calculated results
gRe = [gSort,gRe(nt+1:end)];
cdf = [(N-[N:-1:1]+1)/N*(nt/N)^(iter-1),cdf(nt+1:end)];
zRe = [Z;zRe(nt+1:end,:)];
% Terminate simulation
if iter >= paras.MaxTry || gSort(nt+1) <= 0
pfss = pfss*(size(find(g<0),2)/N);
stopFlag = 1;
break;
else
pfss = pfss*(nt/N);
end;
end;