  The distributions $P(X)$ of singular thermodynamic quantities in an ensemble
of quenched random samples of linear size $l$ at the critical point $T_c$ are
studied by Monte Carlo in two models. Our results confirm predictions of
Aharony and Harris based on Renormalization group considerations. For an
Ashkin-Teller model with strong but irrelevant bond randomness we find that the
relative squared width, $R_X$, of $P(X)$ is weakly self averaging. $R_X\sim
l^{\alpha/\nu}$, where $\alpha$ is the specific heat exponent and $\nu$ is the
correlation length exponent 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 study
the distribution of the pseudo-critical temperatures $T_c(i,l)$ of the ensemble
defined as the temperatures of the maximum susceptibility of each sample. We
find that its variance scales as $(\delta T_c(l))^2 \sim l^{-2/\nu}$ and NOT as
$\sim l^{-d}. We find that $R_\chi$ is reduced by a factor of $\sim 70$ with
respect to $R_\chi (T_c)$ by measuring $\chi$ of each sample at $T_c(i,l)$. We
analyze correlations between the magnetization at criticality $m_i(T_c,l)$ and
the pseudo-critical temperature $T_c(i,l)$ in terms of a sample independent
finite size scaling function of a sample dependent reduced temperature
$(T-T_c(i,l))/T_c$. This function is found to be universal and to behave
similarly to pure systems.
