  In this paper we give a simple proof of the fact that the average over all
dyadic lattices of the dyadic $H^1$-norm of a function gives an equivalent
$H^1$-norm. The proof we present works for both one-parameter and
multi-parameter Hardy spaces.
  The results of such type are known. The first result (for one-parameter Hardy
spces) belongs to Burgess Davis (1980). Also, by duality, such results are
equivalent to the "BMO from dyadic BMO" statements proved by
Garnett-Jones(1982} for one parameter case, and by Pipher-Ward (2008) for
two-parameter case.
  While the paper generalizes these results to the multi-parameter setting,
this is not its main goal. The purpose of the paper is to present an approach
leading to a simple proof, which works in both one-parameter and
multi-parameter cases.
  The main idea of treating square function as a Calderon--Zygmind operator is
a commonplace in harmonic analysis; the main observation, on which the paper is
based, is that one can treat the random dyadic square function this way. After
that, all is proved by using the standard and well-known results about
Calderon--Zygmind operators in the Hilbert-space-valued setting.
  As an added bonus, we get a simple proof of the (equivalent by duality)
inclusion $\text{BMO}\subset \text{BMO}_d$, $H^1_d \subset H^1$ in the
multi-parameter case. Note, that unlike the one-parameter case, the inclusions
in the general situation are far from trivial.
