  We calculate the differential branching fractions and angular distributions
of the rare decays $B^0 \to K^{*0} \mu^+ \mu^-$ and $B_s^0 \to \phi \mu^+
\mu^-$, using for the first time form factors from unquenched lattice QCD. We
focus on the kinematic region where the $K^*$ or $\phi$ recoils softly; there
the newly available form factors are most precise and the nonlocal matrix
elements can be included via an operator product expansion. Our results for the
differential branching fractions calculated in the Standard Model are higher
than the experimental data. We consider the possibility that the deviations are
caused by new physics, and perform a fit of the Wilson coefficients $C_9$ and
$C_9^\prime$ to the experimental data for multiple $B^0 \to K^{*0} \mu^+ \mu^-$
and $B_s^0 \to \phi \mu^+ \mu^-$ observables. In agreement with recent results
from complementary studies, we obtain $C_9-C_9^{SM} = -1.0 \pm 0.6$ and
$C_9^\prime = 1.2 \pm 1.0$, whose deviations from zero would indicate the
presence of non-standard fundamental interactions.
