First, check whether or not the output you are creating is phase sensitive or not. If it's not (e.g. you only use it to control operations and are going to uncompute it later), then you can prepare it much more efficiently by using classical sampling methods and making them reversible. For example, you can use a quantum variant of alias sampling (as explained in Section III of https://arxiv.org/abs/1805.03662) which involves a QROM read circuit, a comparison, and a conditional swap. The rest of this answer assumes you can't do that.
Basically what you want to do is iteratively prepare each qubit to have the right magnitude when conditioned on each possible value of the previous qubits. For example, suppose that when you look at the distribution you're trying to prepare, the total probability of all cases where qubit 0 is ON is 30%. That is to say, $P(b_0) = 0.3$. Then you want to start by rotating qubit 0 until it has a 30% chance of being ON if measured in the Z basis. Now further suppose that, when qubit 0 is OFF, qubit 1 has a 40% chance of being ON and when qubit 0 is ON then qubit 1 has a 70% chance of being on. That is to say, $P(b_1 | !b_0) = 0.4$ and $P(b_1 | b_0) = 0.7$. Then you want to rotate qubit 1 by enough so it has a 40% of being ON, then rotate by an additional delta that increases this chance to 70% but conditioned on qubit 0 being ON.
You continue on in this fashion, rotating each successive qubit by an amount that is conditioned on each possible value of the previous qubits and chosen to make its conditional probabilities match the distribution you're trying to prepare. The controlled rotations have to be decomposed into your gate set. If you correctly optimize the controlled rotations, so that they share controls as much as possible (read: use unary iteration as covered in https://arxiv.org/abs/1805.03662 to generate qubits for all the cases you want to condition on) you will find that this has a total non-Clifford gate count of $O(N \cdot \lg \frac{1}{\epsilon})$ where $N$ is the number of cases in the distribution. If you don't use unary iteration, don't share combined controls between cases, the cost will increase to $O(N \cdot \lg N \cdot \lg \frac{1}{\epsilon})$. For scale, if you were not phase sensitive you'd have only needed $O(N + \lg \frac{1}{\epsilon})$ and the constant factors hiding behind the asymptotic notation would be better.