  Transfer-matrix methods are used to study the probability distributions of
spin-spin correlation functions $G$ in the two-dimensional random-field Ising
model, on long strips of width $L = 3 - 15$ sites, for binary field
distributions at generic distance $R$, temperature $T$ and field intensity
$h_0$. For moderately high $T$, and $h_0$ of the order of magnitude used in
most experiments, the distributions are singly-peaked, though rather
asymmetric. For low temperatures the single-peaked shape deteriorates, crossing
over towards a double-$\delta$ ground-state structure. A connection is obtained
between the probability distribution for correlation functions and the
underlying distribution of accumulated field fluctuations. Analytical
expressions are in good agreement with numerical results for $R/L \gtrsim 1$,
low $T$, $h_0$ not too small, and near G=1. From a finite-size {\it ansatz} at
$T=T_c (h_0=0)$, $h_0 \to 0$, averaged correlation functions are predicted to
scale with $L^y h_0$, $y =7/8$. From numerical data we estimate y=0.875 \pm
0.025$, in excellent agreement with theory. In the same region, the RMS
relative width $W$ of the probability distributions varies for fixed $R/L=1$ as
$W \sim h_0^{\kappa} f(L h_0^u)$ with $\kappa \simeq 0.45$, $u \simeq 0.8$ ;
$f(x)$ appears to saturate when $x \to \infty$, thus implying $W \sim
h_0^{\kappa}$ in $d=2$.
