  We consider long strips of finite width $L \leq 13$ sites of ferromagnetic
Ising spins with random couplings distributed according to the binary
distribution: $P(J_{ij})= {1 \over 2} ( \delta (J_{ij} -J_0) + \delta (J_{ij}
-rJ_0) ) ,\ 0 < r < 1 $. Spin-spin correlation functions $ <\sigma_{0}
\sigma_{R}>$ along the ``infinite'' direction are computed by transfer-matrix
methods, at the critical temperature of the corresponding two-dimensional
system, and their probability distribution is investigated. We show that,
although in-sample fluctuations do not die out as strip length is increased,
averaged values converge satisfactorily. These latter are very close to the
critical correlation functions of the pure Ising model, in agreement with
recent Monte-Carlo simulations. A scaling approach is formulated, which
provides the essential aspects of the $R$-- and $L$-- dependence of the
probability distribution of $\ln <\sigma_{0} \sigma_{R}>$, including the result
that the appropriate scaling variable is $R/L$. Predictions from scaling theory
are borne out by numerical data, which show the probability distribution of
$\ln <\sigma_{0} \sigma_{R}>$ to be remarkably skewed at short distances,
approaching a Gaussian only as $R/L \gg 1$ .
