  We use recent 36 observational Hubble data (OHD) in the redshift range
$0.07\leq z\leq 2.36$, latest \textgravedbl joint light curves\textacutedbl
(JLA) sample, comprised of 740 type Ia supernovae (SNIa) in the redshift range
$0.01\leq z \leq 1.30$, and their joint combination datasets to constrain
anisotropic Bianchi type I (BI) dark energy (DE) model. To estimate model
parameters, we apply Hamiltonian Monte Carlo technique. We also compute the
covariance matrix for BI dark energy model by considering different datasets to
compare the correlation between model parameters. To check the acceptability of
our fittings, all results are compared with those obtained from 9 year WMAP as
well as Planck (2015) collaboration. Our estimations show that at 68\%
confidence level the dark energy equation of state (EOS) parameter for OHD or
JLA datasets alone varies between quintessence and phantom regions whereas for
OHD+JLA dataset this parameter only varies in phantom region. It is also found
that the current cosmic anisotropy is of order $\sim10^{-3}$ which imply that
the OHD and JLA datasets do not put tight constraint on this parameter.
Therefore, to constraint anisotropy parameter, it is necessary to use high
redshif dataset namely cosmic microwave background (CMB). Moreover, from the
calculation of $p$-value associated with $\chi^{2}$ statistic we observed that
non of the $\omega \mbox{BI}$ and flat $\omega\mbox{CDM}$ models rule out by
OHD or JLA datasets. The deceleration parameter is obtained as $q=-0.46^{+0.89
+0.36}_{-0.41 -0.37}$, $q=-0.619^{+0.12 +0.20}_{-0.095 -0.24}$, and
$q=-0.52^{+0.080 +0.014}_{-0.046 -0.15}$ for OHD, SNIa, and OHD+SNIa data
respectively.
