When solving dissipative quantum dynamics using quantum trajectories, it is assumed that for averaging over a sufficiently high number of $N$ trajectories, the full density matrix can be recovered as

$$
\newcommand{\Op}[1]{\hat{#1}}
\newcommand{\ket}[1]{\vert #1 \rangle}
\newcommand{\bra}[1]{\langle #1 \vert}
\newcommand{\Avg}[1]{\langle{#1}\rangle}
\newcommand{\var}{\operatorname{var}}
\Op{\rho}_N(t) = \frac{1}{N} \sum_{i=1}^N \ket{\Psi_i(t)}\bra{\Psi_i(t)}\,.
$$

Any expectation value for an observable $\Op{A}$ with respect to the full density matrix can be calculated as the average of the expectation values from the individual trajectories:

$$
\Avg{\Op{A}}_{\rho_N}
= \frac{1}{N} \sum_{i=1}^N \Avg{\Op{A}}_i\,;
\qquad
\Avg{\Op{A}}_i \equiv \Avg{\Op{A}}_{\Psi_i}\,.
$$

The variance of $\Op{A}$ must be calculated as

$$
\var(\Op{A})_{\rho_N}
= \Avg{\Op{A}^2}_{\rho_N} - \Avg{\Op{A}}_{\rho_N}^2\,.
$$

Note that the variance for the ensemble is in general *not* the average of the
variances of all trajectories, unless the expectation value in each trajectory
is approximately equal. See this pdf file for the full details (tex).