MaxIter =10,000;
Initialize the flux v(0) with random value generator;
for k =1:MaxIter
v(k) = v(k -1);
for l=1:n
Calculate μi and Σi according to equation (10);
Generate the sample θ from the Gaussian distribution
N(μii) within the constraint Φ(∙);
Update vi(k)  with the sample θ, vi(k)=θ ;
end
end
Table 1: Basic procedure for Gibbs sampling algorithm.