0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

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

Last updated at Posted at 2021-08-01

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

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

image.png

最終的に完全な直線に乗るとは限らないので、最小二乗法で直線近似をする。基本的な考え方としては、数の系列$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)

image.png

その2に続く。

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?