  Let $m,n\ge 1$ and $g_{\lambda_1,\lambda_2}^*$ be the bi-parameter
Littlewood-Paley $g_{\lambda}^{*}$-function defined by $$
g_{\lambda_1,\lambda_2}^*(f)(x)= \bigg(\iint_{\R^{m+1}_{+}} \big(\frac{t_2}{t_2
+ |x_2 - y_2|}\big)^{m \lambda_2} \iint_{\R^{n+1}_{+}} \big(\frac{t_1}{t_1 +
|x_1 - y_1|}\big)^{n \lambda_1}|\theta_{t_1,t_2} f(y_1,y_2)|^2 \frac{dy_1
dt_1}{t_1^{n+1}} \frac{dy_2 dt_2}{t_2^{m+1}} \bigg)^{1/2}, \lambda_1>1,\quad
\lambda_2>1 $$ where $\theta_{t_1,t_2} f$ is a non-convolution kernel defined
on $\mathbb{R}^{m+n}$. In this paper, we showed that the bi-parameter
Littlewood-Paley function $g_{\lambda_1,\lambda_2}^*$ was bounded from
$L^2(\R^{n+m})$ to $L^2(\R^{n+m})$. This was done by means of probabilistic
methods and by using a new averaging identity over good double Whitney regions.
