<jupyter_start><jupyter_text>自己組織化<jupyter_code>from scipy.stats import multivariate_normal, cauchy, invgamma
from numpy.random import uniform
import copy


class GaussianNoise(object):
    """多次元ガウス分布
    """
    def __init__(self, covariance, update_target_idx=None):
        self.covariance = covariance
        self.n_state = len(covariance)
        self.update_target_idx = update_target_idx

    def generate(self, n_particles):
        """ノイズの生成
        """
        return (multivariate_normal.rvs(
            np.zeros(self.n_state), self.covariance, n_particles).T)
    
    def pdf(self, y, predicted_obs):
        """確率密度
        """
        r_y = y.reshape(y.size)
        likelihoods = np.empty(predicted_obs.shape[1])
        for i, pred_y in enumerate(predicted_obs.T):
            likelihoods[i] = multivariate_normal.pdf(r_y, 
                                                     pred_y, 
                                                     cov=self.covariance)
        return likelihoods

    def update_parameters(self, particles):
        n_row = len(self.covariance)
        self.covariance[range(n_row), range(n_row)] = \
            np.exp(particles[self.update_target_idx].mean(axis=1))
    

class SelfOrganizationParticleFilter(object):
    def __init__(self, 
                 n_dim_state=None, n_dim_components=None,
                 n_dim_obs=None,
                 transition_func=None, transition_matrix=None,
                 observation_func=None, observation_matrix=None,
                 system_noise=None, observation_noise=None,
                 ar_param_noise=None, n_particles=None,
                 ar_param_idxs=None, ar_param_func=None,
                 Q_param_noise=None, Q_param_idxs=None):
        # 観測データの次元数
        self.n_dim_obs = n_dim_obs
        # 粒子の数
        self.n_particles = n_particles
        # リサンプリングに使用する横櫛（使い回す）
        self.TEETH = np.arange(0, 1, float(1.0)/n_particles)
        # 粒子の単位質量（観測データへの適合度）
        self.weights = np.zeros(n_particles)
        # 粒子
        self.particles = np.zeros((n_dim_state, n_particles))
        # 予測分布（粒子）
        self.predicted_particles = np.zeros((n_dim_state, n_particles))
        self.n_dim_components = n_dim_components

        self.transition_func = transition_func
        self.transition_matrix = transition_matrix
        self.observation_func = observation_func
        self.observation_matrix = observation_matrix
        self.system_noise = system_noise
        self.observation_noise = observation_noise
        self.ar_param_noise = ar_param_noise
        self.ar_param_func = ar_param_func
        self.Q_param_noise = Q_param_noise

        self._init_praticles_distribution()
        self._ar_param_idxs = ar_param_idxs
        self._Q_param_idxs = Q_param_idxs
        self.n_target_params = \
            n_dim_state - len(ar_param_idxs) - len(Q_param_idxs)
        self.start_ar_dim = \
            n_dim_state - len(ar_param_idxs) * 2 - len(Q_param_idxs)
        
    def _init_praticles_distribution(self):
        """粒子の初期化
        x_{0|0}
        """
        self.particles = np.concatenate(
                [self.system_noise.generate(self.n_particles),
                 self.ar_param_noise.generate(self.n_particles),
                 self.Q_param_noise.generate(self.n_particles)],
             axis=0)

    def _update_transition_matrix(self):
        """遷移行列の更新（AR成分の更新）
        """
        for i, ar_idx in enumerate(self._ar_param_idxs):
            self.transition_matrix[self.start_ar_dim, self.start_ar_dim+i] = \
                self.ar_param_func(self.particles[ar_idx].mean())
    
    def update(self, y):
        """リサンプリング
        """
        if isinstance(y, float):
            y = np.array([y])
        self._update_transition_matrix()
        self.system_noise.update_parameters(self.particles)
        self._update_particles()
        self._update_particles_weight(y)
        self._resample()

    def forecast(self, credibility=97.5):
        """長期予測
        """
        self.particles[:self.n_target_params] = \
            self.transition_func(self.particles[:self.n_target_params], 
                                 self.transition_matrix)
        out = self.observation_func(self.particles[:self.n_target_params], 
                                    self.observation_matrix)
        mean_ = out.mean(axis=1)
        # 信用区間の推定
        lower_ = np.percentile(out, 100 - credibility)
        upper_ = np.percentile(out, credibility)
        return mean_, lower_, upper_

    def generate_system_noise(self):
        """v_t"""
        return self.system_noise.generate(self.n_particles)
    
    def _update_one_each_transition_particle(self, added_noise_particles):
        tmp_trans_matrix = copy.deepcopy(self.transition_matrix)
        for ar_coef_vec in self.particles[self._ar_param_idxs].T:
            tmp_trans_matrix[self.start_ar_dim, self.start_ar_dim:] = \
                self.ar_param_func(ar_coef_vec)
            self.predicted_particles[:self.n_target_params] = \
                self.transition_func(added_noise_particles, tmp_trans_matrix)
    
    def _update_particles(self):
        """システムモデルの計算（新たな粒子を生成することで予測分布を計算）
        x_{t|t-1}
        """
        par = self.particles[:self.n_target_params] + self.generate_system_noise()
        self.predicted_particles[:self.n_target_params] = \
            self.transition_func(par, self.transition_matrix)
        self.predicted_particles[self._ar_param_idxs] = \
            self.predicted_particles[self._ar_param_idxs] + \
            self.ar_param_noise.generate(self.n_particles)
        self.predicted_particles[self._Q_param_idxs] = \
            self.predicted_particles[self._Q_param_idxs] + \
            self.Q_param_noise.generate(self.n_particles)
                        
    def _update_particles_weight(self, y):
        """各粒子の重みの計算（予測値が観測値に対してどれだけ適合しているかを計算）
        w_t
        """
        self.weights = self.observation_noise.pdf(
            y, self.observation_func(
                self.predicted_particles[:self.n_target_params], 
                self.observation_matrix))
        
    def _normalize_weights(self):
        """\tild{w}_t"""
        self.weights = self.weights / sum(self.weights)

    def _resample(self):
        """x_{t|t}"""
        self._normalize_weights()

        # accumulate weight
        cum = np.cumsum(self.weights)
        
        # create roulette pointer 
        base = uniform(0, float(1.0)/self.n_particles)
        pointers = self.TEETH + base

        # select particles
        selected_idx = [np.where(cum>=p)[0][0] for p in pointers]
        self.particles = self.predicted_particles[:, selected_idx]

    def predict(self, credibility=97.5):
        """フィルタ系列の推定
        """
        pre_targete_dim = 0
        target_dim = 0
        results = np.empty(len(self.n_dim_components))
        lower = np.empty(len(self.n_dim_components))
        upper = np.empty(len(self.n_dim_components))
        for i, n_dim_each in enumerate(self.n_dim_components):
            target_dim += n_dim_each
            H = self.observation_matrix[:, pre_targete_dim:target_dim]
            matmul = self.observation_func(
                self.predicted_particles[pre_targete_dim:target_dim], 
                H)
            weighted = matmul * self.weights
            mean_ = weighted.sum(axis=1) 

            # 信用区間の推定                        
            lower_ = np.percentile(matmul, 100-credibility)
            upper_ = np.percentile(matmul, credibility)

            results[i] = mean_
            lower[i] = lower_
            upper[i] = upper_
            pre_targete_dim = target_dim
        return results, lower, upper
np.random.seed(555)

n_dim_trend = 2
n_dim_obs = 1
# 季節成分の
n_dim_series = 12
# AR成分の次数
n_dim_ar = 2
# 数値予測の際に成分を分けるため各成分の要素数を格納
n_dim_components = (n_dim_trend, n_dim_series-1, n_dim_ar)
# 粒子数
n_particles = 1000
# システムノイズの分散共分散行列の各要素の値
system_sigma2 = 1500

n_dim_state, F, H, Q = FGHset(n_dim_trend, 
                              n_dim_obs, 
                              n_dim_series, 
                              n_dim_ar, 
                              Q_sigma2=system_sigma2)

# 自己組織化するハイパーパラメタに該当する粒子のindex
ar_param_idxs = [idx for idx in range(n_dim_state, 
                                      n_dim_state + n_dim_ar)]
Q_param_idxs = [idx for idx in range(ar_param_idxs[-1] + 1, 
                                     ar_param_idxs[-1] + 1 + len(Q))]
n_dim_state += len(ar_param_idxs)
n_dim_state += len(Q_param_idxs)

# 観測ノイズの分散共分散行列の各要素の値
# 小さければ観測データに忠実にフィットしやすい
observation_sigma2 = 400

# AR係数はtanhで[-1, 1]の範囲で推定する
# そのためランダムウォークの分散はtanhの逆関数の出力とする
inv_tanh = lambda x: 0.5 * np.log((1 + x) / (1 - x))
ar_params_sigma2 = inv_tanh(0.25)

# システムノイズの分散共分散行列の要素は非負であり、推定時にはexp(x)とする
# そのためランダムウォークの分散は自然対数とする
# 大きくし過ぎると発散しやすい
Q_params_sigma2 = np.log(25)

# システム関数
transition_func = lambda x, F: F.dot(x)
# 観測関数
observation_func = lambda x, H: H.dot(x)
# システムノイズオブジェクトのインスタンス生成
system_noise = GaussianNoise(Q, Q_param_idxs)
# 観測ノイズオブジェクトのインスタンス生成
observation_noise = GaussianNoise(
    np.eye(n_dim_obs) * observation_sigma2)
# AR成分のノイズオブジェクトのインスタンス生成
ar_param_noise = GaussianNoise(
    np.eye(len(ar_param_idxs)) * ar_params_sigma2)
# AR係数を推定する際の復元関数
ar_param_func = lambda x: np.tanh(x)
# システムノイズ共分散行列の要素のノイズオブジェクトのインスタンス生成
Q_param_noise = GaussianNoise(
    np.eye(len(Q_param_idxs)) * Q_params_sigma2)

pf = SelfOrganizationParticleFilter(n_dim_state=n_dim_state, 
                                    n_dim_components=n_dim_components,
                                    n_dim_obs=n_dim_obs,
                                    transition_func=transition_func,
                                    transition_matrix=F,
                                    observation_func=observation_func,
                                    observation_matrix=H,
                                    system_noise=system_noise,
                                    observation_noise=observation_noise,
                                    ar_param_noise=ar_param_noise,
                                    n_particles=n_particles,
                                    ar_param_idxs=ar_param_idxs,
                                    ar_param_func=ar_param_func,
                                    Q_param_noise=Q_param_noise,
                                    Q_param_idxs=Q_param_idxs)


n_train = 120
predicted_trend = np.empty(n_train)
predicted_seasonal = np.empty(n_train)
predicted_ar = np.empty(n_train)
p_lower = np.empty(n_train)
p_upper = np.empty(n_train)
for i, d in enumerate(y[:n_train]):
    pf.update(d)
    results, lower, upper = pf.predict()
    predicted_trend[i] = results[0]
    predicted_seasonal[i] = results[1]
    predicted_ar[i] = results[2]
    
    p_lower[i] = np.sum(lower)
    p_upper[i] = np.sum(upper)
    
rng = range(n_train)
plt.plot(rng, y[:n_train], label="observation")
plt.plot(rng, predicted_trend, '--', label="trend")
plt.plot(rng, predicted_seasonal, ':', label="seasonal")
plt.plot(rng, predicted_ar, '+-', label="ar")
rng = range(n_train)
predicted = predicted_trend + predicted_seasonal + predicted_ar
plt.plot(rng, y[:n_train],label="observation")
plt.plot(rng, predicted, '--', label="predict")
plt.fill_between(rng, p_upper, p_lower, alpha=0.5, label="credible interval")
plt.ylim(0, 550)
infered_value = np.empty(len(y[n_train:]))
inf_lower = np.empty(len(y[n_train:]))
inf_upper = np.empty(len(y[n_train:]))

for i, d in enumerate(y[n_train:]):
    infered_value[i], inf_lower[i], inf_upper[i] = pf.forecast()

rng = range(len(y))
plt.plot(rng, y.values, label="observation")
plt.plot(rng, 
         np.concatenate((predicted, infered_value)), 
         '--', label="forecast")

all_lower = np.concatenate((p_lower, inf_lower))
all_upper = np.concatenate((p_upper, inf_upper))
plt.fill_between(rng, all_upper, all_lower, alpha=0.5, label="credible interval")
plt.ylim(-100, 750)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10006687844842, 'fraction_numerical': 0.008025413810399598, 'mean_word_length': 2.3110988098533074, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '248520', 'license': 'Apache-2.0', 'n_tokens_mistral': 4518, 'n_tokens_neox': 4360, 'n_words': 560, 'path': 'chapter03/3.3.6_particle_filter_so_cp.ipynb', 'repo_name': 'TaigaHasegawa/time-series-analysis-and-anomaly-detection', 'size': '11962'}