Tutorial 12: Posterior Analysis with Uncertain Parameters
Propagate parameter uncertainty through a queueing model. The Posterior solver attaches a Prior distribution to a model parameter, evaluates the model under each alternative, and aggregates the results into a posterior distribution of the chosen performance metrics. The example uses an M/M/1 queue with arrival rate λ = 0.5 and an uncertain service rate μ drawn from a Gaussian-shaped prior over [0.7, 2.5] centred at μ = 1.3.
Block 1: Build the model with a Prior service rate
model = Network('UncertainServiceModel');
source = Source(model, 'Source');
queue = Queue(model, 'Queue', SchedStrategy.FCFS);
sink = Sink(model, 'Sink');
jobClass = OpenClass(model, 'Jobs');
source.setArrival(jobClass, Exp(0.5));
numAlternatives = 30;
serviceRates = linspace(0.7, 2.5, numAlternatives);
priorMean = 1.3; priorStd = 0.4;
priorProbs = exp(-0.5 * ((serviceRates - priorMean) / priorStd).^2);
priorProbs = priorProbs / sum(priorProbs);
alternatives = cell(1, numAlternatives);
for i = 1:numAlternatives
alternatives{i} = Exp(serviceRates(i));
end
queue.setService(jobClass, Prior(alternatives, priorProbs));
model.link(Network.serialRouting(source, queue, sink));Network model = new Network("UncertainServiceModel");
Source source = new Source(model, "Source");
Queue queue = new Queue(model, "Queue", SchedStrategy.FCFS);
Sink sink = new Sink(model, "Sink");
OpenClass jobClass = new OpenClass(model, "Jobs", 0);
source.setArrival(jobClass, new Exp(0.5));
int numAlternatives = 30;
double priorMean = 1.3, priorStd = 0.4;
List<Distribution> alternatives = new ArrayList<>();
double[] priorProbs = new double[numAlternatives];
double sum = 0.0;
for (int i = 0; i < numAlternatives; i++) {
double rate = 0.7 + (2.5 - 0.7) * i / (double) (numAlternatives - 1);
alternatives.add(new Exp(rate));
double z = (rate - priorMean) / priorStd;
priorProbs[i] = Math.exp(-0.5 * z * z);
sum += priorProbs[i];
}
for (int i = 0; i < numAlternatives; i++) priorProbs[i] /= sum;
queue.setService(jobClass, new Prior(alternatives, priorProbs));
model.link(Network.serialRouting(source, queue, sink));val model = Network("UncertainServiceModel")
val source = Source(model, "Source")
val queue = Queue(model, "Queue", SchedStrategy.FCFS)
val sink = Sink(model, "Sink")
val jobClass = OpenClass(model, "Jobs", 0)
source.setArrival(jobClass, Exp(0.5))
val numAlternatives = 30
val priorMean = 1.3; val priorStd = 0.4
val alternatives = (0 until numAlternatives).map {
Exp(0.7 + (2.5 - 0.7) * it / (numAlternatives - 1).toDouble())
}
val rawProbs = (0 until numAlternatives).map {
val rate = 0.7 + (2.5 - 0.7) * it / (numAlternatives - 1).toDouble()
val z = (rate - priorMean) / priorStd
Math.exp(-0.5 * z * z)
}
val sum = rawProbs.sum()
val priorProbs = rawProbs.map { it / sum }.toDoubleArray()
queue.setService(jobClass, Prior(alternatives, priorProbs))
model.link(Network.serialRouting(source, queue, sink))from line_solver import *
import numpy as np
model = Network('UncertainServiceModel')
source = Source(model, 'Source')
queue = Queue(model, 'Queue', SchedStrategy.FCFS)
sink = Sink(model, 'Sink')
jobClass = OpenClass(model, 'Jobs')
source.setArrival(jobClass, Exp(0.5))
num_alternatives = 30
service_rates = np.linspace(0.7, 2.5, num_alternatives)
prior_mean, prior_std = 1.3, 0.4
prior_probs = np.exp(-0.5 * ((service_rates - prior_mean) / prior_std) ** 2)
prior_probs = prior_probs / np.sum(prior_probs)
alternatives = [Exp(rate) for rate in service_rates]
queue.setService(jobClass, Prior(alternatives, prior_probs.tolist()))
model.link(Network.serial_routing([source, queue, sink]))Block 2: Run Posterior over a back-end solver
Posterior wraps an arbitrary back-end solver (here MVA), runs it once per alternative, and aggregates the per-alternative results into prior-weighted means and per-alternative tables.
post = Posterior(model, @SolverMVA);
post.runAnalyzer();
avgTable = post.getAvgTable() % prior-weighted means
postTable = post.getPosteriorTable() % per-alternative results
respDist = post.getPosteriorDist('R', queue, jobClass);
qlenDist = post.getPosteriorDist('Q', queue, jobClass);
utilDist = post.getPosteriorDist('U', queue, jobClass);Posterior post = new Posterior(model, m -> new MVA(m));
post.runAnalyzer();
post.getAvgTable().print(); // prior-weighted means
post.getPosteriorTable().print(); // per-alternative results
SolverPosterior.EmpiricalCDF respDist = post.getPosteriorDist("R", queue, jobClass);
SolverPosterior.EmpiricalCDF qlenDist = post.getPosteriorDist("Q", queue, jobClass);
SolverPosterior.EmpiricalCDF utilDist = post.getPosteriorDist("U", queue, jobClass);val post = Posterior(model) { m -> MVA(m) }
post.runAnalyzer()
post.getAvgTable().print()
post.getPosteriorTable().print()
val respDist = post.getPosteriorDist("R", queue, jobClass)
val qlenDist = post.getPosteriorDist("Q", queue, jobClass)
val utilDist = post.getPosteriorDist("U", queue, jobClass)post = Posterior(model, MVA)
post.run_analyzer()
avg_table = post.get_avg_table() # prior-weighted means
post_table = post.get_posterior_table() # per-alternative results
resp_dist = post.get_posterior_dist('R', queue, jobClass)
qlen_dist = post.get_posterior_dist('Q', queue, jobClass)
util_dist = post.get_posterior_dist('U', queue, jobClass)Block 3: Extract statistics from the posterior CDFs
Each posterior distribution returned by getPosteriorDist/get_posterior_dist is an empirical CDF. In MATLAB and Python the underlying data is exposed as an n×2 matrix (cumulative probability, metric value); in Java and Kotlin the same information is reachable via parallel values, probabilities, and cdf arrays sorted by value. Standard summary statistics follow directly:
respValues = respDist.data(:, 2);
respCdf = respDist.data(:, 1);
respProbs = [respCdf(1); diff(respCdf)];
expectedR = sum(respValues .* respProbs);
medianR = respValues(find(respCdf >= 0.5, 1, 'first'));
fprintf('E[R] = %.4f, median(R) = %.4f\n', expectedR, medianR);double expectedR = respDist.getMean();
double medianR = respDist.values[0];
for (int i = 0; i < respDist.cdf.length; i++) {
if (respDist.cdf[i] >= 0.5) { medianR = respDist.values[i]; break; }
}
System.out.printf("E[R] = %.4f, median(R) = %.4f%n", expectedR, medianR);val expectedR = respDist.getMean()
val medianR = respDist.values[respDist.cdf.indexOfFirst { it >= 0.5 }]
println("E[R] = $expectedR, median(R) = $medianR")resp_values = resp_dist.data[:, 1]
resp_cdf = resp_dist.data[:, 0]
resp_probs = np.concatenate([[resp_cdf[0]], np.diff(resp_cdf)])
expected_R = np.sum(resp_values * resp_probs)
median_R = resp_values[np.searchsorted(resp_cdf, 0.5)]
print(f'E[R] = {expected_R:.4f}, median(R) = {median_R:.4f}')The same recipe applies to queue length, utilization, throughput, or any other metric exposed by the back-end solver. Plotting respValues against respProbs yields a posterior PDF; plotting against respCdf yields the posterior CDF.