Starting from the general effective hamiltonian relevant to the b→s transitions, we derive the expressions for the full angular distributions of the B→K(∗)ℓ1ℓ2 decay modes, as well as for B(Bs→ℓ1ℓ2) (ℓ1=ℓ2). We point out the differences in the treatment of the lepton flavor violating modes with respect to the lepton flavor conserving ones. Concerning the relevant Wilson coefficients we evaluate them in two different scenarios: (i) The (pseudo-)scalar coefficients are obtained using the (pseudo-)scalar coupling extracted from the experimental non-zero value of B(h→μτ), (ii) Revisiting a Z′-model in which the flavor changing neutral couplings are allowed. We provide the numerical estimates of the branching fractions of the above-mentioned modes in both scenarios.