A formulation of the anisotropic pressure magnetohydrodynamic equilibrium problem for three-dimensional plasmas with imposed nested magnetic surfaces is developed based on a bi-Maxwellian model of the distribution function for the energetic particle species. The hot particle distribution function satisfies the constraint B . del F-h = 0. Large parallel and perpendicular anisotropy factors can be explored within the model through the choice of the hot particle perpendicular to parallel temperature ratio T-perpendicular to/T-parallel to. A fixed boundary version of the VMEC code has been adapted to numerically compute three-dimensional anisotropic pressure equilibria. Applications to a 10-field period Heliotron device and a 2-field quasiaxisymmetric stellarator demonstrate that the pressures do not vary significantly around the magnetic surfaces when the total parallel pressure p(parallel to), is larger than its perpendicular counterpart p(perpendicular to). For off-axis hot particle deposition with p(perpendicular to) > p(parallel to), p(perpendicular to) concentrates in the region where the energetic particles are generated. On the other hand, p(parallel to) is distributed roughly uniformly around the flux surfaces in the Heliotron but is localized on the low field side in the quasiaxisymmetric machine. The hot particle density structure correlates more closely with the corresponding perpendicular rather than with parallel pressure. The specific form for the definition of that best correlates with the Shafranov shift is identified.