Wavelet-based regression analysis is widely used mostly for equally-spaced designs. For such designs wavelets are superior to other traditional orthonormal bases because of their versatility and ability to parsimoniously describe irregular functions. If the regression design is random, an automatic solution is not available. Given the observations (X_i, Y_i), i = 1,..., n, we estimate the regression function m(x)=E(Y|X=x) as a series \sum_k \hat c_{jk} \phi_{jk}(x) where \{ \phi_{jk}(x), ~k \in Z \} are scaling functions spanning the multiresolution subspace V_j. We propose a method that utilizes a probabilistic model on X_i's in defining the empirical coefficients \hat c_{jk}. The paper deals with both theoretical and practical aspects of the proposed estimator. We explore MSE convergence rates of the estimator. The performance of the estimator is compared to that of some traditional regression methods. The manuscript is available in PostScript format.