The forward-backward algorithm
This presentation of the forward-backward algorithm is similar to the textbook treatment found in Greenberg (2013) Introduction to Bayesian Econometrics (pp. 193-195).
Let $Y_t = (y_1, \ldots, y_t)$.
The forward part of the algorithm (sometimes called filtering) computes $p(x_t|Y_t)$ from $p(x_{t-1}|Y_{t-1})$ for $t = 1, \ldots, T$.
This forward recursion is initialized with $p(x_0|Y_0) = p(x_0)$, which is the distribution for $x_0$ absent any observations.
The recursion can be decomposed into two steps; namely, the prediction step
\begin{equation}
p(x_t|Y_{t-1}) = \int_S p(x_t|x_{t-1})\,p(x_{t-1}|Y_{t-1})\,d x_{t-1} ,
\end{equation}
followed by the observation step
\begin{equation}
p(x_t|Y_t) \propto p(y_t|x_t)\,p(x_t|Y_{t-1}) .
\end{equation}
The prediction step produces a prior distribution for the observation step in which the observation distribution is the likelihood.
The filtered distributions $p(x_t|Y_t)$ play a central role in the backward part of the algorithm, which where the draw of $X$ is made.
The backward part of the algorithm (sometimes called smoothing) is based on the following factorization of the joint posterior distribution for the state variable:
\begin{equation}
p(X|Y) = p(x_T|Y) \prod_{t=1}^{T-1} p(x_t|Y,X^{t+1}) ,
\end{equation}
where $X^{t+1} = (x_{t+1}, \ldots, x_T)$.
The backward part of the algorithm draws states sequentially beginning with $p(x_T|Y)$ and then using $p(x_{t}|Y,X^{t+1})$ repeatedly. The result is a random draw from $p(X|Y)$.
In order to sample from $p(x_t|Y,X^{t+1})$ we simplify it and express it
in a form that relies on the filtered distributions computed in the forward part.
Begin by breaking up the conditioning variables in the distribution as follows:
\begin{equation}
p(x_t|Y,X^{t+1}) \equiv p(x_t|Y_t,Y^{t+1},x_{t+1},X^{t+2}) .
\end{equation}
As can be seen in the Figure, information contained in either $Y^{t+1}$ or $X^{t+2}$ must pass through $x_{t+1}$ to reach $x_t$.
Therefore, $x_t$ is independent of both $Y^{t+1}$ and $X^{t+2}$ as long as $x_{t+1}$ is in the conditioning list. Consequently,
\begin{equation}
p(x_{t}|Y,X^{t+1}) = p(x_t|Y_t,x_{t+1}) .
\end{equation}
Finish by using Bayes' rule to write
\begin{equation}
p(x_t|Y_t,x_{t+1}) \propto p(x_{t+1}|x_{t})\,p(x_{t}|Y_{t}) .
\end{equation}
The right-hand side provides the key to sampling from $p(X|Y)$. Note in passing it is possible, given the draw of $x_1$, to go one step farther with the backward part of the algorithm and draw $x_0 \sim p(x_0|x_1) \propto p(x_1|x_0)\,p(x_0)$.
Observable covariates
The second example given in the question requires a slight generalization of the setup because the observation distribution depends on an observable covariate, $z_t$. (In the example, $z_t = t$.) It is straightforward to include observable covariates in the observation distribution:
\begin{equation}
p(y_t|x_t,z_t) ,
\end{equation}
where $z_t$ is an observable variable (or a vector of observable variables).
Then $p(Y|X,Z) = \prod_{t=1}^T p(y_t|x_t,z_t)$ and $p(X|Y,Z) \propto p(Y|X,Z)\,p(X)$.
The forward and backward equations are characterized by
\begin{align}
p(x_t|Y_t,Z_t) &\propto p(y_t|x_t,z_t) \int_S p(x_t|x_{t-1})\,p(x_{t-1}|Y_{t-1},Z_{t-1})\,dx_{t-1} \\
p(x_t|Y_t,Z_t,x_{t+1}) &\propto p(x_{t+1}|x_{t})\,p(x_t|Y_t,Z_t) ,
\end{align}
where $Z_t = (z_1, \ldots, z_t)$.
Illustration I: Gaussian model
To derive the forward recursion, begin by assuming
\begin{equation}
p(x_{t-1}|Y_{t-1}) = \textsf{N}(x_{t-1}|\mu_{t-1}, \sigma_{t-1}^2) ,
\end{equation}
where $\mu_{t-1}$ and $\sigma_{t-1}^2$ are sufficient statistics for $Y_{t-1}$.
Then the prediction step implies
\begin{equation}
p(x_t|Y_{t-1}) = \textsf{N}(x_t|M_{t-1}, V_{t-1})
\end{equation}
where
\begin{equation}
M_{t-1} = a+b\,\mu_{t-1}
\qquad\text{and}\qquad
V_{t-1} = s^2 + b^2\,\sigma_{t-1}^2 .
\end{equation}
Consequently, the observation step delivers
\begin{equation}
p(x_t|Y_t) = \textsf{N}(x_t|\mu_t,\sigma_t^2) ,
\end{equation}
where
\begin{equation}
\mu_t = \sigma_t^2 \left(\frac{M_{t-1}}{V_{t-1}} + \frac{(y_t - \alpha)/\beta}{\tau^2/\beta^2}\right)
\qquad\text{and}\qquad
\sigma_t^2 = \left(\frac{1}{V_{t-1}} + \frac{1}{\tau^2/\beta^2}\right)^{-1} .
\end{equation}
The forward recursion is initialized with $(\mu_0,\sigma_0^2)$ and produces $\{(\mu_t,\sigma_t^2)\}_{t=1}^T$.
This recursion is an example of the Kalman filter.
The smoothing distribution is given by
\begin{equation}
p(x_t|Y_t,x_{t+1}) = \textsf{N}(x_t| \mathcal{M}_t(x_{t+1}),\mathcal{V}_t) ,
\end{equation}
where
\begin{equation}
\mathcal{M}_t(x_{t+1}) = \mathcal{V}_t \left(\frac{\mu_t}{\sigma_t^2} + \frac{(x_{t+1}-a)/b}{s^2/b^2}\right)
\qquad\text{and}\qquad
\mathcal{V}_t = \left(\frac{1}{\sigma_t^2} + \frac{1}{s^2/b^2}\right)^{-1} .
\end{equation}
Let $X^{(r)} = (x_1^{(r)},\ldots, x_T^{(r)})$ denote the $r$-th draw of $X$.
The backward part begins with
\begin{equation}
x_T^{(r)} \sim \textsf{N}(x_T|\mu_T,\sigma_T^2)
\end{equation}
followed by
\begin{equation}
x_{T-1}^{(r)} \sim \textsf{N}(x_{T-1}|\mathcal{M}_{T-1}(x_{T}^{(r)}),\mathcal{V}_{T-1})
\end{equation}
and so forth until finally
\begin{equation}
x_1^{(r)} \sim \textsf{N}(x_1|\mathcal{M}_1(x_2^{(r)}),\mathcal{V}_1) .
\end{equation}
Side note: Special example
Let $(a,b,s^2) = (0,1,0)$ so that $x_t = x_0$ for all $t$.
In addition let $(\alpha,\beta,\tau^2) = (0,1,\tau^2)$ so that $p(y_t|x_t) = \textsf{N}(y_t|x_0,\tau^2)$. Evidently $x_0$ is the common (unknown) mean of the observations. (The variance of the observations, $\tau^2$, is known.)
In this case the prediction step produces $M_{t-1} = \mu_{t-1}$ and $V_{t-1} = \sigma_{t-1}^2$, and the observation step delivers $p(x_0|Y_t) = \textsf{N}(x_0|\mu_t,\sigma_t^2)$, where
\begin{equation}
\mu_t = \sigma_t^2 \left(\frac{\mu_{t-1}}{\sigma_{t-1}^2} + \frac{y_t}{\tau^2}\right)
\qquad\text{and}\qquad
\sigma_t^2 = \left(\frac{1}{\sigma_{t-1}^2} + \frac{1}{\tau^2}\right)^{-1} .
\end{equation}
If we further assume $(\mu_0,\sigma_0^2) = (\mu_0,\infty)$, so that there is effectively no prior information regarding $x_0$, then
\begin{equation}
\mu_t = \frac{1}{t} \sum_{s=1}^t y_s
\qquad\text{and}\qquad
\sigma_t^2 = \frac{\tau^2}{t} .
\end{equation}
The initial step of the backward recursion is given by $p(x_0|Y_T) = \textsf{N}(x_0|\mu_T,\sigma_T^2)$. For additional steps, one may show (by taking the limit as $s^2 \to 0$) that $\mathcal{M}_t(x_{t+1}) = x_{t+1}$ and $\mathcal{V}_t = 0$.
Illustration II: Hidden Markov model
I begin with two preliminary comments.
First, given the discrete nature of $S$, some integrals become summations, such as
\begin{equation}
\int_S p(x_t|x_{t-1})\,p(x_{t-1})\,dx_{t-1} \quad\rightarrow\quad \sum_{x_{t-1}\in S} p(x_t|x_{t-1})\,p(x_{t-1}) .
\end{equation}
Second, the observations depend on an observable covariate, $z_t = t$, which is made explicit here:
\begin{equation}
p(y_t|x_t = j, z_t = t) = \textsf{N}(y_t|\alpha_j + \beta_j\,t,\sigma^2) .
\end{equation}
The forward-backward algorithm can be expressed compactly in this case.
Let $L$ and $G$ denote $T \times K$ matrices where
\begin{align}
L_{tj} &= p(y_t|x_t = j, z_t = t) \\
G_{tj} &= p(x_t=j|Y_{t},Z_t) .
\end{align}
Let $G_0$ denote the (row) vector of probabilities for $x_0$ absent any observations, so that $G_{0j} = p(x_0=j)$. The input for the algorithm consists of $(L,Q,G_0)$.
The forward part of the algorithm can be expressed as the recursion
\begin{equation}
G_{t} \propto L_t \circ (G_{t-1} Q) ,
\end{equation}
where $G_t$ is the $t$-th row of $G$, $L_t$ is the $t$-th row of $L$, and "$\circ$" denotes the Hadamard product (component-by-component).
Note the prediction step is embodied in the expression $G_{t-1} Q$.
Turning to the backward part, let $(Q^\top)_k$ denote the $k$-th row of the transpose of $Q$ and let $H_t^k$ denote the row vector of probabilities where
\begin{equation}
H_t^k \propto (Q^\top)_k \circ G_t .
\end{equation}
The probabilities for the backward part are given by
\begin{equation}
p(x_T=j|Y) = G_{Tj} ,
\end{equation}
and
\begin{equation}
p(x_t = j|Y_t,Z_t,x_{t+1} = k) = H_{tj}^k .
\end{equation}
What to do if the parameters are not known?
Thus far we have assumed the parameters of the transition and observation distributions were known. How can one proceed if the parameters are not known? The answer involves using the Gibbs sampler.
We can include the parameters in the notation as follows:
\begin{equation}
p(x_t|x_{t-1},\theta_x)
\qquad\text{and}\qquad
p(y_t|x_t,\theta_y) .
\end{equation}
In the Gaussian example, $\theta_x = (a,b,s^2)$ and $\theta_y = (\alpha,\beta,\tau^2)$. In the HMM example, $\theta_x = Q$ and $\theta_y = (\phi,\sigma^2)$, where $\phi = (\phi_1, \ldots, \phi_K)$ and $\phi_j = (\alpha_j, \beta_j)$.
It is convenient to assume prior independence, where $p(\theta_x,\theta_y) = p(\theta_x)\,p(\theta_y)$. Then the posterior distribution for the complete set of unknowns is
\begin{equation}
p(X,\theta_x,\theta_y|Y) \propto p(Y|X,\theta_y)\,p(X|\theta_x)\,p(\theta_x)\,p(\theta_y) .
\end{equation}
Note that $p(x_1|\theta_x) = \int p(x_1|x_0,\theta_x)\,p(x_0|\theta_x)\,dx_0$ where the marginal distribution for $x_0$ may or may not depend on $\theta_x$.
The joint distribution for the unobserved state variable and the unknown parameters can be characterized by the following collection of full conditional distributions:
\begin{align}
p(X|Y,\theta_x,\theta_y) &\propto p(Y|X,\theta_y)\,p(X|\theta_x) \\
p(\theta_x|X,Y,\theta_y) &\propto p(X|\theta_x)\,p(\theta_x) \\
p(\theta_y|X,Y,\theta_x) &\propto p(Y|X,\theta_y)\,p(\theta_y) .
\end{align}
The Gibbs sampler produces draws from the joint posterior distribution by sequentially drawing from the full conditional distributions. Given $(X^{(r-1)}, \theta_x^{(r-1)}, \theta_y^{(r-1)})$, the next draw in the chain can be computed as follows:
\begin{align}
X^{(r)} &\sim p(X|Y,\theta_x^{(r-1)}, \theta_y^{(r-1)}) \\
\theta_x^{(r)} &\sim p(\theta_x|Y,X^{(r)},\theta_y^{(r-1)}) \\
\theta_y^{(r)} &\sim p(\theta|Y,X^{(r)},\theta_x^{(r)}) .
\end{align}
The draws from $X^{(r)}$ are made using the forward-backward algorithm described above, treating $(\theta_x^{(r-1)},\theta_y^{(r-1)})$ as known.
The way in which draws from the other two distributions are made will depend on the specific case.
Gibbs sampler for the Gaussian example
Recall $\theta_x = (a,b,s^2)$ and $\theta_y = (\alpha,\beta,\tau^2)$. Let $p(\theta_x) \propto 1/s^2$ and let $p(\theta_y) \propto 1/\tau^2$.
It is convenient to augment the list of unknowns with $x_0$. In this case we have
\begin{equation}
p(X, x_0,\theta_x,\theta_y|Y) \propto p(Y|X,\theta_y)\, p(X|x_0,\theta_x)\,p(x_0)\,p(\theta_x)\,p(\theta_y) .
\end{equation}
We have already discussed the possibility of using the forward-backward algorithm to draw from
\begin{equation}
p(X,x_0|-) \propto p(Y|X,\theta_y)\, p(X|x_0,\theta_x)\,p(x_0)
\end{equation}
by taking one additional backward step.
Then
\begin{equation}
p(\theta_x|-) \propto p(X|x_0,\theta_x)\,p(\theta_x) = \frac{1}{s^2} \prod_{t=1}^T \textsf{N}(x_t|a+b\,x_{t-1}, s^2) .
\end{equation}
If either $\mu_0$ or $\sigma_0^2$ depends on $\theta_x$, then the factor $p(x_0|\theta_x)$ must be included in the kernel for $p(\theta_x|-)$. Finally,
\begin{equation}
p(\theta_y|-) \propto p(Y|X,\theta_y)\,p(\theta_y) = \frac{1}{\tau^2} \prod_{t=1}^T \textsf{N}(y_t|\alpha + \beta\,x_t, \tau^2) .
\end{equation}
One can make standard Gibbs-sampler draws of $\theta_x$ and $\theta_y$.
Gibbs sampler for HMM example
Recall $\theta_x = Q$ and $\theta_y = (\phi, \sigma^2)$.
We begin with the draws of $Q$. Let $T_j$ denote the number of times $x_t = j$ for $t = 1, \ldots, T-1$. Let $C$ denote a $K\times K$ matrix where $C_{jk}$ is the number of times $(x_t, x_{t+1}) = (j,k)$ for $t = 1, \ldots, T-1$, and let $C_j$ denote the $j$-th row of $C$.
Then
\begin{equation}
p(C_j|Q_j,T_j) = \textsf{Multinomial}(C_j|T_j,Q_j) .
\end{equation}
If $p(Q_j) = \textsf{Dirichlet}(Q_j|\alpha_j)$, then
\begin{equation}
p(\theta_x|-) = \prod_{j=1}^K p(Q_j|C_j) ,
\end{equation}
where
\begin{equation}
p(Q_j|C_j) = \textsf{Dirichlet}(Q_j|\alpha_j + C_j) .
\end{equation}
We now turn to the draws of $\theta_y$. Let
\begin{equation}
\mathcal{Y}_j = \{y_t : x_t = j\} .
\end{equation}
(Note $\cup_{j=1}^K \mathcal{Y}_j = Y$.)
Let
\begin{equation}
p(\mathcal{Y}_j|\phi_j,\sigma^2) = \prod_{y_t \in \mathcal{Y}_j} p(y_t|\phi_j,\sigma^2) .
\end{equation}
Then
\begin{equation}
p(\phi|-) = \prod_{j=1}^K p(\phi_j|-)
\end{equation}
where
\begin{equation}
p(\phi_j|-) \propto p(\phi_j) \prod_{t\in\tau_j} p(y_t|x_t,\theta) .
\end{equation}
For $\sigma^2$ we have
\begin{equation}
p(\sigma^2|-) \propto p(Y|X,\theta_y)\,p(\sigma^2) .
\end{equation}
