I am able to reproduce the PSPhot-SDSS bias using simulated data. I generate a collection of fake stars, distributed with a luminosity function: logN ~ alpha*mag. I've set alpha to 0.3 to enhance the number of bright sources in this example. I smooth the image with a PSF, then add noise (gain = 1). I then run psphot on these. For the following plots, I am using a Gaussian PSF, and I set the sky level to the values 10, 30, 100, 300, 1000, 3000 for the sequence of runs. I use the same input star fluxes for each run. I am running psphot with no aperture correction, a Gaussian PSF model, and only using the 'ensemble' psf fits (no non-linear fits).