Probability of randomly designated subsets cover the universe
There is no polynomial-time algorithm that computes an exact solution (i.e. the probability expressed as a ratio of two nonnegative integers $a/b$) to your problem unless $\textbf{P} = \textbf{NP}$.
I will show that such an exact solution would provide a polynomial time algorithm for set-cover.

Let $(U, S \subseteq \mathscr{P}(U), k)$ be an instance of set-cover. Let $m = |S|$.

Let
$$
 \Gamma = \left\{ K \in \mathscr{P}(S) : U =\bigcup_{X \in K} X \right\}
$$

For each $h \in \mathbb{N}$ let:

$$
 \nu_h = |\{K \in \Gamma : |K|=h\}|
$$

Observe that our instance of set-cover is satisfiable if and only if $\nu_k \neq 0$.

Now, let $T$ be a set, constructed as in the problem statement. $T$ covers $U$ if and only if $T \in \Gamma$. Therefore:

$$
\mathbb{Pr}\left[ U =\bigcup_{X \in T} X \right] = \mathbb{Pr}[T \in \Gamma] = \sum_{X\in\Gamma}p^{|X|}(1-p)^{m-|X|} = \sum_{h=0}^m \nu_h p^h(1-p)^{m-h}
$$

The rightmost expression is a polynomial over the ring $\mathbb{Q}[p]$, which means that its coefficients can be found by computing its value at $m-1$ points and applying an algorithm such as the finite differences method.