[step:Prove the tensorization inequality by induction over the number of coordinates]
We prove the theorem for all $n\in\mathbb{N}$ by induction. For $n=1$, the assertion is an equality:
\begin{align*}
\operatorname{Ent}_{\mu_1}(f)=\int_{E_{\varnothing}}\operatorname{Ent}_{\mu_1,1}(f)\,d\mu_{\varnothing}.
\end{align*}
Here $E_{\varnothing}$ is a one-point space and $\mu_{\varnothing}$ is its unit mass.
Assume the theorem holds for $n-1$ factors. Define the product space $E'$ and product probability measure $\mu'$ by
\begin{align*}
E' := \prod_{i=1}^{n-1}E_i.
\end{align*}
\begin{align*}
\mu' := \bigotimes_{i=1}^{n-1}\mu_i.
\end{align*}
Define the marginal function $g:E'\to[0,\infty]$ by
\begin{align*}
g(x'):=\int_{E_n} f(x',z_n)\,d\mu_n(z_n).
\end{align*}
Applying the two-factor decomposition with $(S,\rho)=(E_n,\mu_n)$ and $(T,\sigma)=(E',\mu')$ gives
\begin{align*}
\operatorname{Ent}_\mu(f)=\int_{E'}\operatorname{Ent}_{\mu_n,n}(f)(x')\,d\mu'(x')+\operatorname{Ent}_{\mu'}(g).
\end{align*}
The first term on the right is non-negative and $\operatorname{Ent}_\mu(f)<\infty$ by hypothesis, hence $\operatorname{Ent}_{\mu'}(g)<\infty$. Therefore the induction hypothesis applies to $g$ on $E'$, and gives
\begin{align*}
\operatorname{Ent}_{\mu'}(g)\leq \sum_{i=1}^{n-1}\int_{E'_{-i}}\operatorname{Ent}_{\mu_i,i}(g)(y_{-i})\,d\mu'_{-i}(y_{-i}),
\end{align*}
where $E'_{-i}:=\prod_{1\leq j\leq n-1,\ j\ne i}E_j$ and $\mu'_{-i}:=\bigotimes_{1\leq j\leq n-1,\ j\ne i}\mu_j$.
Since $\operatorname{Ent}_\mu(f)<\infty$, the mass of $f$ is finite:
\begin{align*}
\int_E f\,d\mu<\infty.
\end{align*}
Fix $i\in\{1,\dots,n-1\}$. By Tonelli's theorem, for $\mu'_{-i}$-almost every $y_{-i}\in E'_{-i}$ the two-coordinate slice has finite total mass:
\begin{align*}
\int_{E_i}\int_{E_n} f(z_i,y_{-i},z_n)\,d\mu_n(z_n)\,d\mu_i(z_i)<\infty.
\end{align*}
Fix such a $y_{-i}$. For each $z_n\in E_n$, define the slice map $u_{z_n}:E_i\to[0,\infty]$ by
\begin{align*}
u_{z_n}(z_i):=f(z_i,y_{-i},z_n).
\end{align*}
The $i$-th slice of $g$ at $y_{-i}$ is the mixture map $g_{y_{-i}}:E_i\to[0,\infty]$ defined by
\begin{align*}
g_{y_{-i}}(z_i):=\int_{E_n}u_{z_n}(z_i)\,d\mu_n(z_n).
\end{align*}
By the convexity estimate from the first step,
\begin{align*}
\operatorname{Ent}_{\mu_i,i}(g)(y_{-i})\leq \int_{E_n}\operatorname{Ent}_{\mu_i}(u_{z_n})\,d\mu_n(z_n).
\end{align*}
Integrating this inequality with respect to $d\mu'_{-i}(y_{-i})$ gives
\begin{align*}
\int_{E'_{-i}}\operatorname{Ent}_{\mu_i,i}(g)(y_{-i})\,d\mu'_{-i}(y_{-i})\leq \int_{E_{-i}}\operatorname{Ent}_{\mu_i,i}(f)(x_{-i})\,d\mu_{-i}(x_{-i}).
\end{align*}
Here the exceptional $y_{-i}$ set has $\mu'_{-i}$-measure zero, and the identification of the right-hand side is Tonelli's theorem for the non-negative coordinate entropy terms together with product-measure associativity.
Combining the last estimates with the decomposition of $\operatorname{Ent}_\mu(f)$ yields
\begin{align*}
\operatorname{Ent}_\mu(f)\leq \sum_{i=1}^{n}\int_{E_{-i}}\operatorname{Ent}_{\mu_i,i}(f)(x_{-i})\,d\mu_{-i}(x_{-i}).
\end{align*}
This is the desired tensorization inequality.
[/step]