It's not that hard. If you have just planar or angular light sources, you can think of them as one light source split into multiple chunks and the only thing to deal with is how to sample this multi-light and how to compute the PDF of the resulting samples.
First, you need to setup the picking probability $P(l)$ for each light source $l$. The picking probabilities can be any non-negative numbers, you just have to make sure that are always non-zero if the contributions of the respective lights are non-zero. However, the closer the probabilities are to the actual (relative) light contribution, the better the overall performance of your Monte Carlo estimator will be.
Contribution estimates may be:
- $1$: naive, uniform probability; good to start with; used by PBRT
- $power$: better
- $power / distance^2$: even better; closer light sources have greater contributions
- …and other improvements: make zero if not visible at all, multiply by the cosine of the inclination angle, etc.
Note that to get the light-picking probabilities, you need to normalise the contribution estimates to the sum of all of them:
You could even use more different picking strategies if you find they work just for certain type of scenes and combine them using multiple importance sampling technique…
When sampling the lights, you just select one light $l$ with the appropriate probability $P(l)$ and then sample the light surface/angle with its own sampling strategy (with PDF $p_l$). The resulting PDF $p$ will be:
p(x)=P(l) * p_l(x)
When computing the light sampling PDF for given direction (as needed for multiple importance sampling), you just find the light source at the direction, evaluate the PDF of the lights own strategy and, again, multiply it by the probability of picking this particular light.
Now you have a sampling strategy which gives you PDF which can be used in the Monte Carlo estimators (e.g. $f(x)/p(x)$) and the resulting estimator will be unbiased.