Tutorial 13: 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));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]))#include "line/lang/prior.h"
Network model("UncertainServiceModel");
Source source(model, "Source");
Queue queue(model, "Queue", SchedStrategy::FCFS);
Sink sink(model, "Sink");
OpenClass jobClass(model, "Jobs");
source.set_arrival(jobClass, Exp(0.5));
const std::size_t numAlternatives = 30;
const double priorMean = 1.3, priorStd = 0.4;
std::vector<Distrib<double>> alternatives;
std::vector<double> priorProbs(numAlternatives);
double norm = 0.0;
for (std::size_t i = 0; i < numAlternatives; ++i) {
const double rate = 0.7 + (2.5 - 0.7) * double(i) /
double(numAlternatives - 1);
alternatives.push_back(Exp(rate));
const double z = (rate - priorMean) / priorStd;
priorProbs[i] = std::exp(-0.5 * z * z);
norm += priorProbs[i];
}
for (double& p : priorProbs) p /= norm;
queue.set_service(jobClass, prior_discrete(alternatives, priorProbs));
Routing P;
P.set(source, queue, 1.0);
P.set(queue, sink, 1.0);
model.link(P);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);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)#include "line/solvers/mva/solver_mva_runner.h"
#include "line/solvers/uq/solver_uq.h"
const uq::UqStageSolver<double> stage =
[](const qn::NetworkStruct<double>& sn) {
mva::MvaOptions options;
return mva::solver_mva_run_analyzer(
sn, options, Matrix<double>());
};
uq::UqOptions options;
const uq::UqSolution<double> post =
uq::solver_uq_run_analyzer(model, stage, options);
// post.avg contains prior-weighted means; post.points contains the
// result at each alternative, paired with the weights in post.weights.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 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);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}')// queue and jobClass are 1-based station/class indices in the C++ API.
const auto respDist = uq::uq_posterior_cdf(post, "R", queue, jobClass);
const auto moments = uq::uq_moments(post, "R", queue, jobClass);
const double expectedR = moments.mean;
double medianR = respDist.values.front();
for (std::size_t i = 0; i < respDist.cdf.size(); ++i) {
if (respDist.cdf[i] >= 0.5) {
medianR = respDist.values[i];
break;
}
}
std::cout << "E[R] = " << expectedR
<< ", median(R) = " << medianR << '\n';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.