Is this interpretation of the EM algorithm correct? Suppose we want to estimate the parameter $\theta$ of the distribution $p(z;\theta)$ of a random variable $Z$. If we had the i.i.d samples $z_1, z_2, \dots, z_N$, we could easily do this using maximum likelihood estimation. However, suppose we are only given the samples $x_1,x_2,\dots,x_N$ of another random variable $X$, where $X$ and $Z$ are related via the joint distribution $p(x,z)$. Since $p(x,z) = p(x \mid z) \cdot p(z;\theta)$, then $p(x,z) = p(x,z;\theta)$.
The objective of the EM algorithm is to estimate the samples $z_1, z_2, \dots, z_N$ using the given samples $x_1,x_2,\dots,x_N$, and then use these estimated samples to perform maximum likelihood estimation to estimate $\theta$. However, it turns out that estimating the samples $z_1, z_2, \dots, z_N$ using the samples $x_1,x_2,\dots,x_N$ requires knowledge of $\theta$, which is what we are trying to estimate in the first place. Therefore, we initialize a value for $\theta$, estimate the samples $z_1, z_2, \dots, z_N$, then estimate a new value for $\theta$, and repeat the process until convergence.
More formally, if
$$
p(z_1, z_2, \dots, z_N;\theta)
$$
is the likelihood function for $\theta$, then to maximize this likelihood function, we need to estimate the samples $z_1, z_2, \dots, z_N$ as $\hat{z}_1, \hat{z}_2, \dots, \hat{z}_N$ such that
$$
\hat{z}_i = \mathbb{E}[Z \mid x_i] \quad \forall i \in \{1,2,\dots,N\}
$$
We could then plug these estimates into the likelihood function,
$$
p(\hat{z}_1, \hat{z}_2, \dots, \hat{z}_N;\theta)
$$
to determine $\theta$. However, as each $\hat{z}_i$ depends on $\theta$, we need to use the iterative process mentioned above.
Is this interpretation correct?

The main misconception in my question is that we are not trying to estimate $\theta$ as a parameter for $p(z;\theta)$, but we are trying to estimate $\theta$ as a parameter for $p(x,z;\theta)$. More precisely, we originally want to estimate $\theta$ as a parameter for $p(x;\theta)$, but if some of the data that we need to do this is missing, then we can split $p(x;\theta)$ into $p(x,z;\theta)$, where $z$ represents the missing values. Therefore, the main motivation for the EM algorithm is to deal with the missing data that is needed to estimate $\theta$. See the paper A Gentle Tutorial of the EM Algorithm and its Application to Parameter Estimation for Gaussian Mixture and Hidden Markov Models by Jeff Bilmes, 1998, for more details.

To add to @ChristianHennig's answer, consider the case where the joint distribution of $X$ and $Z$ is a bivariate normal distribution with parameters $\mu_X,\sigma_X^2,\mu_Z,\sigma_Z^2,$ and $\rho$, where
$$
\rho = \frac{\text{Cov}(X,Z)}{\sqrt{\sigma_X^2 \sigma_Z^2}}
$$
is the correlation coefficient. This means that
$$
p(x,z) = \frac{1}{2\pi \sigma_X \sigma_Z \sqrt{1 - \rho^2}} \cdot \exp\left(-\frac{1}{2(1-\rho^2)}\left[\left(\frac{x - \mu_X}{\sigma_X}\right)^2 + \left(\frac{z - \mu_Z}{\sigma_Z}\right)^2 - 2\rho \frac{(x - \mu_X)(z - \mu_Z)}{\sigma_X \sigma_Z} \right]\right)
$$
In the E-step of the EM algorithm, the conditional expectation
$$
\mathbb{E}[\log p(x,Z) \mid x]
$$
with respect to $Z$ is computed. Since
$$
\log p(X,Z) = -\log \left(2\pi \sigma_X \sigma_Z \sqrt{1 - \rho^2}\right) -\frac{1}{2(1-\rho^2)}\left[\left(\frac{X - \mu_X}{\sigma_X}\right)^2 + \left(\frac{Z - \mu_Z}{\sigma_Z}\right)^2 - 2\rho \frac{(X - \mu_X)(Z - \mu_Z)}{\sigma_X \sigma_Z} \right]
$$
then
$$
\mathbb{E}[\log p(x,Z) \mid x] = -\log \left(2\pi \sigma_X \sigma_Z \sqrt{1 - \rho^2}\right) -\frac{1}{2(1-\rho^2)}\left[\left(\frac{x - \mu_X}{\sigma_X}\right)^2 + \mathbb{E}\left[\left(\frac{Z - \mu_Z}{\sigma_Z}\right)^2 \mid x\right] - 2\rho \frac{(x - \mu_X)(\mathbb{E}[Z \mid x] - \mu_Z)}{\sigma_X \sigma_Z} \right]
$$
Since $(\cdot)^2$ is a convex function, then by Jensen's inequality,
$$
\mathbb{E}\left[\left(\frac{Z - \mu_Z}{\sigma_Z}\right)^2 \mid x\right] \geq \left(\frac{\mathbb{E}[Z \mid x] - \mu_Z}{\sigma_Z}\right)^2
$$
and so
$$
\mathbb{E}[\log p(x,Z) \mid x] \leq -\log \left(2\pi \sigma_X \sigma_Z \sqrt{1 - \rho^2}\right) -\frac{1}{2(1-\rho^2)}\left[\left(\frac{x - \mu_X}{\sigma_X}\right)^2 + \left(\frac{\mathbb{E}[Z \mid x] - \mu_Z}{\sigma_Z}\right)^2 - 2\rho \frac{(x - \mu_X)(\mathbb{E}[Z \mid x] - \mu_Z)}{\sigma_X \sigma_Z} \right]
$$
Note that
$$
\log p(x,\mathbb{E}[Z \mid x]) = -\log \left(2\pi \sigma_X \sigma_Z \sqrt{1 - \rho^2}\right) -\frac{1}{2(1-\rho^2)}\left[\left(\frac{x - \mu_X}{\sigma_X}\right)^2 + \left(\frac{\mathbb{E}[Z \mid x] - \mu_Z}{\sigma_Z}\right)^2 - 2\rho \frac{(x - \mu_X)(\mathbb{E}[Z \mid x] - \mu_Z)}{\sigma_X \sigma_Z} \right]
$$
and so
$$
\mathbb{E}[\log p(x,Z) \mid x] \leq \log p(x,\mathbb{E}[Z \mid x])
$$
which means that $\log p(x,\mathbb{E}[Z \mid x])$ is an upper-bound for $\mathbb{E}[\log p(x,Z) \mid x]$ in the case of a bivariate normal distribution. However, there is no guarantee that $\log p(x,\mathbb{E}[Z \mid x])$ is a lower-bound for $\log p(x)$, and so increasing $\log p(x,\mathbb{E}[Z \mid x])$ does not imply that $\log p(x)$ increases.
Furthermore, the proposed algorithm in the original question can be implemented for the bivariate normal distribution case, but it will not converge. Here is an implementation in Python:
import numpy as np

# initialize random number generator
rng = np.random.default_rng(seed = 42)

# mean vector
mu = np.array([1,2])

# covariance matrix
sigma = np.array([[3,0],
                  [0,5]])

# number of samples
N = 1000

# observed (x) and unobserved (z) samples
x,z = rng.multivariate_normal(mu,sigma,N).T

# initial guesses for the mean vector and covariance matrix. Note that
# sigma_hat is matrix-multiplied by its transpose to ensure it is
# symmetric and positive definite
mu_hat = mu + rng.standard_normal(mu.shape) * 0.1
sigma_hat = sigma + rng.standard_normal(sigma.shape) * 0.1
sigma_hat = sigma_hat.T @ sigma_hat

# the complete dataset. The first column consists of all the observed
# samples of x. The second column consists of initial guesses of z
Y = rng.standard_normal((N,2))
Y[:,0] = x

for _ in range(100):
    # correlation coefficient estimate
    rho_hat = sigma_hat[0,1] / np.sqrt(sigma_hat[0,0] * sigma_hat[1,1])
    
    # estimates of z's
    Y[:,1] = mu_hat[1] + rho_hat * sigma_hat[1,1] * ((Y[:,0] - mu_hat[0])/sigma_hat[0,0])
    
    # estimate of mean vector
    mu_hat = Y.mean(axis = 0)
    
    # estimate of covariance matrix
    sigma_hat = (Y - mu_hat).T @ (Y - mu_hat) / Y.shape[0]

