Are you constrained to using Euler angles? If you can decide freely what representation to use, quaternions would be preferable. In that representation, it's very easy to find which of two quaternions representing a rotation (which differ only by a sign) is closer to a given one.
With Euler angles, if you want the representation to be continuous, you can let the angles range over $\mathbb R$ instead of restricting them to a single (half-)period. However, choosing the nearest of the equivalent representations will then be more complicated than with quaternions. Quite generally, a lot of things that are nasty, complicated and potentially numerically unstable with Euler angles become nicer and easier when you use quaternions.
[Edit in response to the comments:]
There are three equivalences, one obvious, another less obvious and a third only applicable in certain circumstances.
The obvious one is that you can always add multiples of $2\pi$ to any of the angles; if you let them range over $\mathbb R$, which you must if you want to get continuous curves, this corresponds to using $\mathbb R^3$ as the parameter space instead of the quotient $(\mathbb R/2\pi\mathbb Z)^3$. This equivalence is easy to handle since you can change the three angles independently, that is, if you change one of them by a multiple of $2\pi$, you directly get the same rotation without changing the other two parameters.
What's less obvious is that (referring to this image) the transformation $(\alpha,\beta,\gamma)\rightarrow(\alpha+\pi,\pi -\beta,\gamma+\pi)$ leads to the same rotation. (This is why, in order to get unique angles, $\beta$ has to be limited to an interval of length $\pi$, not $2\pi$.)
A third equivalence comes into play only if $\beta\equiv0\pmod{\pi/2}$, since in this case $\alpha$ and $\gamma$ apply to the same axis and changing $\alpha+\gamma$ doesn't change the rotation. If your rotations are arbitrary and have no reason to have $\beta\equiv0$, you won't need to consider this case, though it may cause numerical problems if you get close to $\beta\equiv0$ (which is one good reason to use quaternions instead of Euler angles).
These three transformations generate all values of the Euler angles that are equivalent to each other. Remember that you also have to consider combinations of them, e.g. you can add multiples of $2\pi$ in $(\alpha+\pi,-\beta,\gamma+\pi)$ to get further equivalent angles.
Okay, so you have an East-North-Down coordinate system as the fixed frame, and you have the phone's Up-Left-Out (of the screen) set of axes, and the quaternion is given in terms of the fixed coordinate system's basis vectors. I'll call the fixed axes $\hat x, \hat y, \hat z$ and the phone's axes $\hat x', \hat y', \hat z'$.
If you're using the phone as a camera, then I imagine it's the phone's Out-of-screen direction that is most like the direction the camera is facing (or rather, the negative of this). So the vector $-\hat z'$ should be the direction that the camera is facing. You should then project this vector onto the $xy$-plane by zeroing out its $\hat z$ component: $P_{xy}(\hat z') = \hat z' - \hat z (\hat z \cdot \hat z')$. You can then renormalize the result and find the angle in the $xy$-plane by basic trig. This should give you the compass direction that the camera is pointing (as an angle relative to east).
Angle off the horizon is easier. Take $-\hat z' \cdot \hat z = \cos \theta$. $\theta$ is then the angle off the vertical. $\theta - \pi/2$ should then give the angle you want: if $\theta = \pi/2$, then the camera is level. If it's larger, then the camera is tilted upwards. That checks out.
I admit, I'm not sure why you want to know if the device is horizontal to the ground; it's not clear to me how this relates to the direction of the camera, and so I'm not sure what answer to give for that part.
Best Answer
Empirically, by trying combinations of yaw, pitch, and roll angles ranging independently between $-30$ and $30$ degrees, and converting each sequence of Euler rotations to an equivalent quaternion using the same algorithm as in How to convert Euler angles to Quaternions and get the same Euler angles back from Quaternions?, the real part of the quaternion ranges from approximately $0.88388347648318455$ to $1$ and the other parts each range from approximately $-0.30618621784789724$ to $0.30618621784789724.$
Note that $5/\sqrt{32} \approx 0.88388347648318441$ and $\sqrt{3/32} \approx 0.30618621784789726,$ so I'll use those numbers as the extreme values of the parts of the quaternion and blame the inexactness of numerical computation for the differences with the values in the previous paragraph.
In principle, if we could sample uniformly over all unit-magnitude quaternions in these ranges of values, and reject the ones that convert back to Euler angles outside your desired range, that would be an answer. The tricky part is the uniform sampling over only quaternions with unit magnitude. It might be nice to get those by rejection sampling within the volume $$\left[-\sqrt{\tfrac{3}{32}},\sqrt{\tfrac{3}{32}}\right] \times \left[-\sqrt{\tfrac{3}{32}},\sqrt{\tfrac{3}{32}}\right] \times \left[-\sqrt{\tfrac{3}{32}},\sqrt{\tfrac{3}{32}}\right] \times \left[\tfrac{5}{\sqrt{32}},1\right],$$ that is, letting each of the imaginary parts of the quaternion range from $-\sqrt{\tfrac{3}{32}}$ to $\sqrt{\tfrac{3}{32}}$ and letting the real part range from $\tfrac{5}{\sqrt{32}}$ to $1,$ but since mathematically the unit quaternions occupy zero percent of this volume, that's not a good approach. Instead, you might try to get a sample over a set that includes all the quaternions of the form $r\mathbf q$, where $0 < r \leq 1$ and $\mathbf q$ is one of the unit quaternions producing your desired set of rotations. (Since $r\mathbf q$ describes the same rotation as $\mathbf q$ when $r$ is real and non-zero, this gives the same distribution of rotations as if we sampled only over the suitable unit quaternions.) You can get this larger set of quaternions by sampling over $$\left[-\sqrt{\tfrac{3}{32}},\sqrt{\tfrac{3}{32}}\right] \times \left[-\sqrt{\tfrac{3}{32}},\sqrt{\tfrac{3}{32}}\right] \times \left[-\sqrt{\tfrac{3}{32}},\sqrt{\tfrac{3}{32}}\right] \times \left[0,1\right]$$ (similar to the previous volume but letting the real part range all the way from $0$ to $1$) and rejecting any quaternion whose magnitude is greater than $1$ (less than $1$ is OK) or whose Euler angles are not within the desired range. This will be better than letting all the parts of the quaternion range over $[-1,1],$ because you'll be able to accept a much larger percentage of the random quaternions within that region; but I think you may do even better by accepting only quaternions in a smaller range of magnitudes, perhaps $\frac12$ to $1,$ using an initial sample over a smaller region.
So I might suggest something like this: set $r_\min$ to a suitable value such as $\frac12.$ Independently choose three numbers each in the range $\left[-\sqrt{\tfrac{3}{32}},\sqrt{\tfrac{3}{32}}\right]$ and one number in the range $\left[\tfrac{5}{\sqrt{32}}r_\min,1\right].$ Assign the first three numbers to the imaginary parts of a quaternion and the last number to the real part. If the resulting quaternion has a magnitude less than $r_\min$ or greater than $1,$ reject it and start over. Convert the resulting quaternion to Euler angles; if any of those angles is less than $-30$ degrees or greater than $30$ degrees, reject the quaternion and start over. But if the quaternion passes all these tests, use it as your next random rotation.
I think any method similar to this is going to end up rejecting a large percentage of the quaternions, but as long as you can get a reasonable percentage to be in the "accept" region then it can be a practical way to get a uniform distribution over the desired set of rotations.
For reference, here is the python code I used to explore the range of unit quaternions that you might use:
The last two lines of the output were