折返しのある系列を巻き戻す

こんな感じで直線（的なもの）が折り返されているときに、折返しを巻き戻したい。

最終的に完全な直線に乗るとは限らないので、最小二乗法で直線近似をする。基本的な考え方としては、数の系列$x_{k}$と、折返し何周目なのかを表す系列$y_{k}$を使い、折返しを元に戻したときにできるだけ直線に近くなるようにする。図では連続的に点があるが、途中の点が抜けてもいいように考慮する。

系列$x_{1},x_{2},\ldots,x_{N}$について、非減少系列$y_{1},y_{1},\ldots,y_{N}(y_{k} \leq y_{k + 1},\ y_{k} \in Z)$があったとき、重み$w_{1},\ldots,w_{N} \in \left\{ 0,1\right\}$を既知として

$$w_{n}(x_{n} + \alpha y_{n})$$

ができるだけ等差数列に近くなるように$\alpha$を決めたい。

最終目標となる等差数列を

$$a + bn$$

とおく。この時の二乗誤差は

$$e^{2} = \sum_{i = 1}^{n}{w_{n}\left( x_{n} + \alpha y_{n} - a - bn \right)}^{2}\quad　(1)$$

$$\frac{\partial e^{2}}{\partial\alpha} = \sum_{n = 1}^{N}{2w_{n}y_{n}(x_{n} + \alpha y_{n} - a - bn)} = 0\quad　(2)$$

$$\frac{\partial e^{2}}{\partial a} = - \sum_{n = 1}^{N}{2w_{n}(x_{n} + \alpha y_{n} - a - bn)} = 0\quad　(3)$$

$$\frac{\partial e^{2}}{\partial b} = - \sum_{n = 1}^{N}{2nw_{n}(x_{n} + \alpha y_{n} - a - bn)} = 0\quad　(4)$$

ここで次のような数を考える。

$$\sum_{n = 1}^{N}{w_{n}x_{n}} = S_{X}$$

$$\sum_{n = 1}^{N}{w_{n}y_{n}} = S_{Y}$$

$$\sum_{n = 1}^{N}{w_{n}x_{n}y_{n}} = S_{\text{XY}}$$

$$\sum_{n = 1}^{N}{w_{n}y_{n}^{2}} = S_{YY}$$

$$\sum_{n = 1}^{N}{w_{n}nx_{n}} = M_{X}$$

$$\sum_{n = 1}^{N}{w_{n}ny_{n}} = M_{Y}$$

$$\sum_{n = 1}^{N}w_{n} = W$$

$$\sum_{n = 1}^{N}{w_{n}n} = \Sigma$$

$$\sum_{n = 1}^{N}{w_{n}n^{2}} = \Sigma_{2}$$

(3)より、

$$S_{X} + \alpha S_{Y} - a W - b\Sigma = 0 \quad(5)$$

(2)より

$$S_{XY} + \alpha S_{YY} - a S_{Y} - b M_{Y} = 0\quad　(6)$$

(4)より

$$M_{X} + \alpha M_{Y} - a\Sigma - b\Sigma_{2} = 0\quad　(7)$$

以上より

$$
\begin{pmatrix}S_{XY} \\ S_{X} \\ M_{X}\end{pmatrix} +
\begin{pmatrix}
S_{YY} & - S_{Y} & - M_{Y} \\
S_{Y} & - W & - \Sigma \\
M_{Y} & - \Sigma & - \Sigma_{2}
\end{pmatrix}\begin{pmatrix}
\alpha \\
a \\
b
\end{pmatrix} = 0$$

これを解いて$\alpha$を求めるとよい。プログラムで書けばこんな感じ。

findlinear <- function(pos,x,y) {
  if (max(y) == 0) {
    # only one segment
    model <- lm(x~pos,data=data.frame(x=x,pos=pos))
    return(c(0,model$coefficients))
  }
  sx <- sum(x)
  sy <- sum(y)
  sxy <- sum(x*y)
  syy <- sum(y^2)
  mx <- sum(pos*x)
  my <- sum(pos*y)
  v <- -c(sxy,sx,mx)
  n <- length(pos)  
  sn <- sum(pos)
  sn2 <- sum(pos^2)
  m <- matrix(c(syy,-sy,-my,sy,-n,-sn,my,-sn,-sn2),
              nrow=3,ncol=3,byrow=TRUE)
  solve(m,v)
}

やってみるとこんな感じ

> x
 [1]   0   3   6   9  12 -14 -11  -8  -5  -2   1   4   7  10  13 -13 -10  -7  -4
[20]  -1   2   5   8  11  14 -12  -9  -6  -3   0   3   6   9  12
> y
 [1] 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3
> n
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
[26] 26 27 28 29 30 31 32 33 34
> findlinear(n,x,y)
[1] 29 -3  3
> plot(x+29*y)

その2に続く。