Weighted least squares with angular data 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.
 A: 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.
$$
