Suppose I have a system whose state is $\Theta=(\theta_1,\theta_2,\ldots,\theta_n)$, where $\theta_i\in[-\pi,\pi)$ (i.e., they are angles). I'd like to determine the most likely estimate of $\Theta$ given $m$ measurements $z_{ij}$ of the difference between two of the angles, where $z_{ij}\in[-\pi,\pi)$. In other words,
$$ z_{ij} = \begin{cases} (\theta_i-\theta_j)\texttt{ % }2\pi & (\theta_i-\theta_j)\texttt{ % }2\pi < \pi, \\ [(\theta_i-\theta_j)\texttt{ % }2\pi]-2\pi & (\theta_i-\theta_j)\texttt{ % }2\pi \geq \pi. \end{cases} $$ There are more measurements then angles in $\Theta$ (i.e., $m>n$), and each measurement has variance $\sigma^2_{ij}$.
For example, given the measurements $Z = (z_{01}, z_{12}, z_{12}, z_{13}, z_{14}, z_{23}, z_{24}, z_{34}, z_{34})$, what is the least squares estimate of $\Theta=(\theta_1,\theta_2,\theta_3,\theta_4)$? Note that I've assumed that multiple measurements between two angles are possible, and $z_{01}$ is a direct measurement of $\theta_1$, needed to make the solution unique.
EDITED THIS IS INCORRECT Tecnically this is straightforward. You specify the equations $$ z_{ij}=\theta_i-\theta_j + \sigma_{ij} \epsilon_{ij} $$ and try to minimize the error with respect to $\theta$s $$ \sum_{ij}\epsilon_{ij}^2=\sum_{ij}\frac{1}{\sigma_{ij}^2}(z_{ij}-\theta_i+\theta_j)^2 $$ with your contsraints on $\theta$s.
EDITED The formula above is incorrect since $z_{ij}=\theta_i-\theta_j$ by modulo $2\pi$.This is more difficult. How we can fix it? One way to avoid this to take $\sin$ from both parts and mininmize sum of squares of $$ \sum (\sin (z_{ij})-\sin(\theta_i-\theta_j))^2 $$ But this is awkward, nonlinear and $\sigma_{ij}$ are not properly specified here. But I think I see a more elegant way to solve it. I will try to give it later.
This problem can be converted into standard linear weighted least squares by making the variable substitution $y_{ij} = z_{ij} + y_{(i-1)i}$ for every $z_{ij}$. For the stated example, $Y=(y_{01}, y_{12}, \ldots, y_{34})$, where $$ \begin{align} y_{01} &= z_{01}, \\ y_{12} &= z_{12} + z_{01}, \\ y_{12} &= z_{12} + z_{01}, \\ y_{13} &= z_{13} + z_{01}, \\ y_{14} &= z_{13} + z_{01}, \\ y_{23} &= z_{23} + z_{12} + z_{01}, \\ y_{24} &= z_{24} + z_{12} + z_{01}, \\ y_{34} &= z_{34} + z_{34} + z_{12} + z_{01}, \\ y_{34} &= z_{34} + z_{34} + z_{12} + z_{01}, \end{align} $$ and the variances of $y_{ij}$ are (for example) $\sigma^2_{y_{23}} = \sigma^2_{z_{23}} + \sigma^2_{z_{01}} + \sigma^2_{z_{01}}$, etc., assuming each $z_{ij}$ is an independent measurement. One can then adjust each $y_{ij}$ using the modulo operator (as defined in the question) to ensure $y_{ij}\in[-\pi,\pi)$.
The maximum likelihood estimate of $\Theta$ is equivalent to minimizing the negative log-likelihood of all the measurements. In other words, it is the $\Theta$ that minimizes $$ F = (Y - H\Theta)^T\Omega(Y - H\Theta), $$ where $\Omega=$diag$(\sigma^2_{y_{01}},\sigma^2_{y_{12}},\ldots)^{-1}$ and $$ H = \begin{bmatrix} -1 & 0 & 0 & 0 \\ 0 & -1 & 0 & 0 \\ 0 & -1 & 0 & 0 \\ 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & -1 \\ 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & -1 \\ 0 & 0 & 0 & -1 \\ 0 & 0 & 0 & -1 \end{bmatrix} $$ The solution is obtained by setting $\frac{\partial F}{\partial\Theta}=0$ and solving for $\Theta$, yielding $$ \Theta = (H^T\Omega H)^{-1}H^T\Omega Y. $$