We evaluate by Monte Carlo simulations various singular thermodynamic quantities X, for an ensemble of quenched random Ising and Ashkin-Teller models. The measurements are taken at T c and we study how the distributions P (X) (and, in particular, their relative squared width, R X ) over the ensemble depend on the system size l. The Ashkin-Teller model was studied in the regime where bond randomness is irrelevant and we found weak self averaging; R X ∼ l α ν → 0, where α < 0 and ν are the exponents (of the pure model fixed point) governing the transition. For the site dilute Ising model on a cubic lattice, known to be governed by a random fixed point, we find that R X tends to a universal constant independent of the amount of dilution (no self averaging). However this constant is different for canonical and grand canonical disorder. We identify the pseudo-critical temperatures T c (i, l) of each sample i, as the temperature at which the susceptibility reaches its maximal value. The distribution of these T c (i, l) over the ensemble was investigated; we found that its variance scales as (δT c (l)) 2 ∼ l − 2 ν . These results are in agreement with the recent predictions of Aharony and Harris. Our previously proposed finite size scaling ansatz for disordered systems was tested and found to hold. When we fit the data obtained for many samples of different sizes by a sample-independent form, the resulting scaling function was found to be universal and to behave similarly to pure systems. We did observe that to describe deviations from this universal function we do need sample-dependent scaling functions. These deviations are, however, relatively small and this led us to an interesting side result: sample-to-sample fluctuations of χ max , the susceptibility measured at T c (i, l), are smaller by a factor of 70 than those of χ(T c ). This indicates that to obtain a fixed statistical error it may be more computationally efficient to measure χ max . 05.50.+q, 75.10Nr, 75.40Mg, 75.50.Lk