I think this is one of the most important part when you write a physically-based renderer.
The solution I came up with is to put a sphere to the scene with a constant background radiance. Solve the rendering equation analytically for constant radiance which you have chosen. There will be no recursion in the integral since the rays will not bounce more than one. After that, compare result of the rendered image and the analytic solution. Does this make sense?
Are there any ways to this more effectively?