We now obtain posterior samples for $U_1 := \theta$ and $U_2 := \sigma$ by the Gibbs sampling algorithm (instead of using PyMC). In other words, we need to sample from the target density $f_{\theta, \sigma \mid \text{data}}(t, s)$. Note that, by Bayes rule, this posterior density equals:
\begin{align*}
f_{\theta, \sigma \mid \text{data}}(t, s) = \frac{f_{\theta, \sigma}(t, s) f_{\text{data} \mid \theta = t, \sigma = s}(\text{data})}{\int \int f_{\theta, \sigma}(t, s) f_{\text{data} \mid \theta = t, \sigma = s}(\text{data}) dt ds}
\end{align*}
The numerator above is simply $prior \times likelihood$ and the denominator is the normalizing constant which makes the density (over $t, s$) integrate to one. One often ignores the denominator and writes the above relation using the "proportional to" sign ($\propto$) as:
\begin{align*}
f_{\theta, \sigma \mid \text{data}}(t, s) \propto \text{Prior} \times \text{Likelihood} = f_{\theta, \sigma}(t, s) f_{\text{data} \mid \theta = t, \sigma = s}(\text{data})
\end{align*}
The prior is given by
\begin{align*}
f_{\theta, \sigma}(t, s) &= f_{\theta}(t) \times f_{\sigma}(s), ~~\text{assuming prior independence betweeen } \theta ~ \text{ and } ~\sigma \\
&= \frac{I\{-80 \leq t \leq 80\}}{160} \times f_{\log \sigma}(\log s) \frac{d \log s}{ds}, ~~ \text{ where the identity } f_{\sigma}(s) = f_{\log \sigma}(\log s) \frac{d \log s}{ds} ~ \text{ follows from the Jacobian formula} \\
&= \frac{I\{-80 \leq t \leq 80\}}{160} \times \frac{I\{-10 \leq \log s \leq 10\}}{20} \frac{1}{s} \\
&= \frac{I\{-80 \leq t \leq 80\}}{160} \times \frac{I\{\exp(-10) \leq s\leq \exp(10)\}}{20} \frac{1}{s} \\
&\propto I\{-80 \leq t \leq 80\} I\{\exp(-10) \leq s \leq \exp(10)\} \frac{1}{s}
\end{align*}
where, in the last step, we ignored the constant factors $160$ and $20$ by replacing the equality sign with the proportional sign.
The Likelihood is given by
\begin{align*}
f_{\text{data} \mid \theta = t, \sigma = s}(\text{data}) &= f_{y_1, \dots, y_n \mid \theta = t, \sigma = s}(y_1, \dots, y_n) \\
&= f_{y_1 \mid \theta = t, \sigma = s}(y_1) \times f_{y_2 \mid \theta = t, \sigma = s}(y_2) \times \dots \times f_{y_n \mid \theta = t, \sigma = s}(y_n) \\
&= \frac{1}{s\sqrt{2\pi}} \exp \left(-\frac{(y_1 - t)^2}{2s^2} \right) \times \frac{1}{s\sqrt{2\pi}} \exp \left(-\frac{(y_2 - t)^2}{2s^2} \right) \times \dots \times \frac{1}{s\sqrt{2\pi}} \exp \left(-\frac{(y_n - t)^2}{2s^2} \right) \\
&= \frac{1}{s^n (\sqrt{2 \pi})^n} \exp \left(-\frac{(y_1 - t)^2 + (y_2 - t)^2 + \dots + (y_n - t)^2}{2s^2} \right).
\end{align*}
The constant factor $\frac{1}{(\sqrt{2\pi})^n}$ does not depend on the parameters $t, s$ and so it can be ignored in proportionality leading to
\begin{align*}
f_{\text{data} \mid \theta = t, \sigma = s}(\text{data}) \propto \frac{1}{s^n} \exp \left(-\frac{(y_1 - t)^2 + (y_2 - t)^2 + \dots + (y_n - t)^2}{2s^2} \right).
\end{align*}
Combining the above derived formulae for the prior and likelihood, we obtain the posterior as:
\begin{align*}
& f_{\theta, \sigma \mid \text{data}}(t, s) \\ &\propto \text{Prior} \times \text{Likelihood} \\
&\propto I\{-80 \leq t \leq 80\} I\{\exp(-10) \leq s \leq \exp(10)\} \frac{1}{s} \times \frac{1}{s^n} \exp \left(-\frac{(y_1 - t)^2 + (y_2 - t)^2 + \dots + (y_n - t)^2}{2s^2} \right) \\
&= I\{-80 \leq t \leq 80\} I\{\exp(-10) \leq s \leq \exp(10)\} \frac{1}{s^{n+1}} \exp \left(-\frac{(y_1 - t)^2 + (y_2 - t)^2 + \dots + (y_n - t)^2}{2s^2} \right).
\end{align*}
The above relation is true up to proportionality. If we want to have equality, we have to normalize the right hand side leading to:
\begin{align*}
f_{\theta, \sigma \mid \text{data}}(t, s) = \frac{I\{-80 \leq t \leq 80\} I\{\exp(-10) \leq s \leq \exp(10)\} \frac{1}{s^{n+1}} \exp \left(-\frac{(y_1 - t)^2 + (y_2 - t)^2 + \dots + (y_n - t)^2}{2s^2} \right)}{\int_{-80}^{80} \int_{e^{-10}}^{e^{10}} \frac{1}{s^{n+1}} \exp \left(-\frac{(y_1 - t)^2 + (y_2 - t)^2 + \dots + (y_n - t)^2}{2s^2} \right) dt ds}
\end{align*}
We now use the Gibbs sampler to generate samples $(t^{(0)}, s^{(0)})$, $(t^{(1)}, s^{(1)})$, $\dots$, $(t^{(T)}, s^{(T)})$ from this somewhat complicated looking bivariate density. The main thing to figure out for Gibbs are the conditional distributions $\theta \mid \sigma = s, \text{data}$ and $\sigma \mid \theta = t, \text{data}$. For $f_{\theta \mid \sigma = s, \text{data}}(t)$, write
\begin{align*}
f_{\theta \mid \sigma = s, \text{data}}(t) &\propto f_{\theta, \sigma \mid \text{data}}(t, s) \\
&\propto I\{-80 \leq t \leq 80\} I\{\exp(-10) \leq s \leq \exp(10)\} \frac{1}{s^{n+1}} \exp \left(-\frac{(y_1 - t)^2 + (y_2 - t)^2 + \dots + (y_n - t)^2}{2s^2} \right).
\end{align*}
$f_{\theta \mid \sigma = s, \text{data}}(t)$ is a density in the variable $t$ which means that any multiplicative factor not depending on $t$ can be ignored in proportionality. We thus get
\begin{align*}
f_{\theta \mid \sigma = s, \text{data}}(t) &\propto f_{\theta, \sigma \mid \text{data}}(t, s) \\
&\propto I\{-80 \leq t \leq 80\} \exp \left(-\frac{(y_1 - t)^2 + (y_2 - t)^2 + \dots + (y_n - t)^2}{2s^2} \right).
\end{align*}
Using the elementary equality
\begin{align*}
(y_1 - t)^2 + (y_2 - t)^2 + \dots + (y_n - t)^2 = (y_1 - \bar{y})^2 + (y_2 - \bar{y})^2 + \dots + (y_n - \bar{y})^2 + n(\bar{y} - t)^2,
\end{align*}
we get
\begin{align*}
f_{\theta \mid \sigma = s, \text{data}}(t)
&\propto I\{-80 \leq t \leq 80\} \exp \left(-\frac{(y_1 - \bar{y})^2 + (y_2 - \bar{y})^2 + \dots + (y_n - \bar{y})^2 + n(\bar{y} - t)^2}{2s^2} \right) \\
&= I\{-80 \leq t \leq 80\} \exp \left(-\frac{(y_1 - \bar{y})^2 + (y_2 - \bar{y})^2 + \dots + (y_n - \bar{y})^2}{2 s^2} \right) \exp \left(-\frac{n(\bar{y} - t)^2}{2s^2} \right)
\end{align*}
The multiplicative factor $\exp \left(-\frac{(y_1 - \bar{y})^2 + (y_2 - \bar{y})^2 + \dots + (y_n - \bar{y})^2}{2 s^2} \right)$ does not depend on $t$ and can be ignored in proportionality leading to the simpler expression:
\begin{align*}
f_{\theta \mid \sigma = s, \text{data}}(t) \propto I\{-80 \leq t \leq 80\} \exp \left(-\frac{n(\bar{y} - t)^2}{2s^2} \right) = I\{-80 \leq t \leq 80\} \exp \left(-\frac{(t - \bar{y})^2}{2\frac{s^2}{n}} \right)
\end{align*}
In other words, $\theta \mid \sigma = s, \text{data}$ is the normal distribution $N(\bar{y}, \frac{s^2}{n})$ conditioned to lie between $[-80, 80]$. Generating a random sample from this distribution is easy: sample from $N(\bar{y}, \frac{s^2}{n})$. Keep the sample unless it lies outside $[-80, 80]$. If it lies outside $[-80, 80]$, just throw this sample away and repeat. Note that, in our dataset, $\bar{y} = 17.615$ and $s^2/n$ will be generally small, so the samples from $N(\bar{y}, s^2/n)$ will almost always lie between $-80$ and $80$.
For the other conditional $\sigma \mid \theta = t, \text{data}$, write
\begin{align*}
f_{\sigma \mid \theta = t, \text{data}}(s) &\propto f_{\theta, \sigma \mid \text{data}}(t, s) \\
&\propto I\{-80 \leq t \leq 80\} I\{\exp(-10) \leq s \leq \exp(10)\} \frac{1}{s^{n+1}} \exp \left(-\frac{(y_1 - t)^2 + (y_2 - t)^2 + \dots + (y_n - t)^2}{2s^2} \right).
\end{align*}
Now we can ignore the factor $I\{-80 \leq t \leq 80\}$ which does not depend on $s$. We then get
\begin{align*}
f_{\sigma \mid \theta = t, \text{data}}(s)
&\propto I\{\exp(-10) \leq s \leq \exp(10)\} \frac{1}{s^{n+1}} \exp \left(-\frac{(y_1 - t)^2 + (y_2 - t)^2 + \dots + (y_n - t)^2}{2s^2} \right).
\end{align*}
To parse this density, first of all, we can replace the somewhat complicated indicator function by the simpler $I\{s > 0\}$ because $\exp(-10) \approx 0$ and $\exp(10) \approx \infty$. The rest of the terms involve a power term and an exponential term resembling a Gamma density. However, this seems like a Gamma density in $1/s^2$ as the exponent involves $1/s^2$. We can convert it back to a regular Gamma density by calculating the density of $1/\sigma^2$ using the Jacobian formula:
\begin{align*}
f_{\frac{1}{\sigma^2} \mid \theta = t, \text{data}}(u) = f_{\sigma \mid \theta = t, \text{data}}\left(\sqrt{\frac{1}{u}} \right)\frac{1}{2u^{3/2}} \propto I\{u > 0\} u^{(n-2)/2}\exp \left(-\frac{(y_1 - t)^2 + (y_2 - t)^2 + \dots + (y_n - t)^2}{2} u \right)
\end{align*}
We thus get
\begin{align*}
\frac{1}{\sigma^2} \mid \theta = t, \text{data} \sim \text{Gamma} \left(\text{shape} = \frac{n}{2}, \text{rate} = \frac{(y_1 - t)^2 + (y_2 - t)^2 + \dots + (y_n - t)^2}{2} \right)
\end{align*}
