  We present and assess a Bayesian method to interpret gravitational wave
signals from binary black holes. Our method directly compares gravitational
wave data to numerical relativity simulations. This procedure bypasses
approximations used in semi-analytical models for compact binary coalescence.
In this work, we use only the full posterior parameter distribution for generic
nonprecessing binaries, drawing inferences away from the set of NR simulations
used, via interpolation of a single scalar quantity (the marginalized
log-likelihood, $\ln {\cal L}$) evaluated by comparing data to nonprecessing
binary black hole simulations. We also compare the data to generic simulations,
and discuss the effectiveness of this procedure for generic sources. We
specifically assess the impact of higher order modes, repeating our
interpretation with both $l\le2$ as well as $l\le3$ harmonic modes. Using the
$l\le3$ higher modes, we gain more information from the signal and can better
constrain the parameters of the gravitational wave signal. We assess and
quantify several sources of systematic error that our procedure could
introduce, including simulation resolution and duration; most are negligible.
We show through examples that our method can recover the parameters for equal
mass, zero spin; GW150914-like; and unequal mass, precessing spin sources. Our
study of this new parameter estimation method demonstrates we can quantify and
understand the systematic and statistical error. This method allows us to use
higher order modes from numerical relativity simulations to better constrain
the black hole binary parameters.
