Skip to content

Instantly share code, notes, and snippets.

@shunsukeaihara
Last active August 14, 2023 06:28
Show Gist options
  • Save shunsukeaihara/8080887 to your computer and use it in GitHub Desktop.
Save shunsukeaihara/8080887 to your computer and use it in GitHub Desktop.
Change finder algorithm
# -*- coding: utf-8 -*-
import numpy as np
import math
def LevinsonDurbin(r, lpcOrder):
"""
from http://aidiary.hatenablog.com/entry/20120415/1334458954
"""
a = np.zeros(lpcOrder + 1,dtype=np.float64)
e = np.zeros(lpcOrder + 1,dtype=np.float64)
a[0] = 1.0
a[1] = - r[1] / r[0]
e[1] = r[0] + r[1] * a[1]
lam = - r[1] / r[0]
for k in range(1, lpcOrder):
lam = 0.0
for j in range(k + 1):
lam -= a[j] * r[k + 1 - j]
lam /= e[k]
U = [1]
U.extend([a[i] for i in range(1, k + 1)])
U.append(0)
V = [0]
V.extend([a[i] for i in range(k, 0, -1)])
V.append(1)
a = np.array(U) + lam * np.array(V)
e[k + 1] = e[k] * (1.0 - lam * lam)
return a, e[-1]
class _SDAR_1Dim(object):
def __init__(self, r, order):
self._r = r
self._mu = np.random.random()
self._sigma = np.random.random()
self._order = order
self._c = np.zeros(self._order+1)
def update(self,x,term):
assert len(term) >= self._order, "term must be order or more"
term = np.array(term)
self._mu = (1 - self._r) * self._mu + self._r * x
for i in range(1,self._order):
self._c[i] = (1-self._r)*self._c[i]+self._r * (x-self._mu) * (term[-i]-self._mu)
self._c[0] = (1-self._r)*self._c[0]+self._r * (x-self._mu)*(x-self._mu)
what,e = LevinsonDurbin(self._c,self._order)
xhat = np.dot(-what[1:],(term[::-1] - self._mu))+self._mu
self._sigma = (1-self._r)*self._sigma + self._r * (x-xhat) * (x-xhat)
return -math.log(math.exp(-0.5 *(x-xhat)**2/self._sigma)/((2 * math.pi)**0.5 * self._sigma**0.5)),xhat
class _ChangeFinderAbstract(object):
def _add_one(self,one,ts,size):
ts.append(one)
if len(ts)==size+1:
ts.pop(0)
def _smoothing(self,ts):
return sum(ts)/float(len(ts))
class ChangeFinder(_ChangeFinderAbstract):
def __init__(self, r = 0.5, order = 1, smooth=7):
assert order > 0, "order must be 1 or more."
assert smooth > 2, "term must be 3 or more."
self._smooth = smooth
self._smooth2 = int(round(self._smooth/2.0))
self._order = order
self._r = r
self._ts = []
self._first_scores = []
self._smoothed_scores = []
self._second_scores = []
self._sdar_first = _SDAR_1Dim(r,self._order)
self._sdar_second = _SDAR_1Dim(r,self._order)
def update(self,x):
score = 0
predict = x
predict2 = 0
if len(self._ts) == self._order:#第一段学習
score,predict = self._sdar_first.update(x,self._ts)
self._add_one(score,self._first_scores,self._smooth)
self._add_one(x,self._ts,self._order)
second_target = None
if len(self._first_scores) == self._smooth:#平滑化
second_target = self._smoothing(self._first_scores)
if second_target and len(self._smoothed_scores) == self._order:#第二段学習
score,predict2 = self._sdar_second.update(second_target,self._smoothed_scores)
self._add_one(score,
self._second_scores,self._smooth2)
if second_target:
self._add_one(second_target,self._smoothed_scores, self._order)
if len(self._second_scores) == self._smooth2:
return self._smoothing(self._second_scores),predict
else:
return 0.0,predict
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment