時系列データにおける異常検知

2018/01/17
このエントリーをはてなブックマークに追加


時系列データにおける異常検知

はじめに

カブクで機械学習エンジニアをしている大串正矢です。今回は時系列データにおける異常検知について書きます。

背景

時系列データの異常検知の手法は多種多様に存在していますがウェブ上にまとまった情報が日本語でないため記述することにしました。以前に記述した異常検知の基礎に異常検知の基礎的な内容が記述されているのでその内容を踏まえた上で読んで頂けると理解がしやすいと思います。

時系列データにおける異常検知手法

通常の異常検知手法との違いは時系列データなのでデータが独立ではなくデータ間に依存性が存在することです。前回紹介した手法はデータが独立であることを仮定しているため使用できません。
異なる手法が必要になります。
代表的な手法を下記に記述します。

  • 動的時間伸縮法
    • 利点:シンプル
    • 欠点:ノイズに弱い
  • 特異スペクトル変換
    • 利点:ノイズに頑健
    • 欠点:計算コストが高い
  • 自己回帰モデル
    • 利点:周期性があるデータに強い。異常検知以外にも時系列予測、変化点検知に使用可能
    • 欠点:次数が決めうちのため状態が一定なものでないと適用が難しい
  • 状態空間モデル
    • 利点:状態が変化しても対応可能かつノイズに強い。自己回帰モデルと同様に時系列予測、変化点検知に使用可能
    • 欠点:周期性があり安定しているデータには対応しづらい

下記の点から今回は自己回帰モデルと状態空間モデルにフォーカスして記述します。

  • 自己回帰モデルによる異常検知は予測モデルを自己回帰モデルから他のモデルに置き換えても使用可能なため拡張性が高い
  • 状態空間モデルは内部状態が変化しても対応可能なため、状態が変化する機器の異常検知にも対応可能。例:車のセンサーなど

自己回帰モデル

  • モデルの推定
    • モデルを定義し、正常なデータから学習
  • 異常度の定義
    • 実波形と予測波形の差
  • 閾値の設定

データについて

下記の図のようにデータをwindowで区切って扱います。このようにすることで行列としてデータを扱うことができ、行列計算などの処理が容易になります。ここで区切った各ウィンドウをεと置きます。



適切なwindowの幅を設定してデータを分割します。適切なwindowの幅はAICと呼ばれる手法で決めれますがこの点は後述します。

自己回帰モデルの学習

下記に簡単な図を示します。学習用の波形と予測用の波形を用意しています。この予測波形をうまく模倣できるように自己回帰モデルを推定します。



この時のwindowの幅が次数rになります。

\begin{align}
y^{t} \equiv \epsilon^{(t+r)}
\end{align}

\begin{align}
\boldsymbol{x}^{t} \equiv \left( \begin{array}{c} \epsilon^{(t)} \ \epsilon^{(t+1)}
\ \vdots \ \epsilon^{(t+r-1)} \end{array} \right)
\end{align}

\begin{align}
N \equiv T – r
\end{align}

\begin{align}
y^{t} \approx \alpha_1\epsilon^{(t-r)} + \alpha_2\epsilon^{(t-r+1)} + … + \alpha_{r-1}\epsilon^{(t-2)} + \alpha_{r}\epsilon^{(t-1)}
\end{align}

このモデルの予測値と実測値の予測誤差を正規分布でモデル化します。



\begin{align}
p(y^{n} | \boldsymbol{\alpha}, \sigma) = \mathcal{N}(y^{n} | \boldsymbol{\alpha}^T\boldsymbol{x}, \sigma^2)
\end{align}

上式の平均を下記で置き換え

\begin{align}
\mu=\boldsymbol{\alpha}^{T}\boldsymbol{x}^{(n)}
\end{align}

出力値を下記で置き換え

\begin{align}
x’=y^{n}
\end{align}

この正規分布に対して対数尤度関数の式を立てます。下記の式は正規分布で使用した対数尤度です。

\begin{align}
\frac{1}{2\sigma^2}(x’ – \mu)^2 + \frac{1}{2}\ln(2\pi\sigma^2)
\end{align}

ここで元の時系列データを下記のように変更します。

\begin{align}
D = {(\boldsymbol{x}^{(1)}, y^{(1)}), …, (\boldsymbol{x}^{(N)}, y^{(N)}) }
\end{align}

N個分のデータを考慮する必要があるので予測誤差は下記の式になります。

\begin{align}
L(\boldsymbol{\alpha}, \sigma^2 | D) = -\frac{N}{2}\ln(2\pi\sigma^2) – \frac{1}{2\sigma^2}\sum^N_{n=1}[y^{n} – \boldsymbol{\alpha}^T\boldsymbol{x}^{(n)}]^{2}
\end{align}

ではこの予測誤差の分布が最小になれば予測と実測の値が最小になるので最小な値を取れるようなαとσを導出します。
αで微分して0と等価とすると

\begin{align}
0 = -\frac{1}{2\sigma^2}(-\boldsymbol{X}\boldsymbol{y} + \alpha \boldsymbol{X}\boldsymbol{X}^T)
\end{align}

\begin{align}
\boldsymbol{\alpha} = [\boldsymbol{X}\boldsymbol{X}^T]^{-1}\boldsymbol{X}\boldsymbol{y}
\end{align}

σで微分して0と等価とすると

\begin{align}
0 = -\frac{N4\pi\sigma}{4\pi\sigma^2} + \frac{2}{2\sigma^{3}}\sum^N_{n=1}[y^{n} – \boldsymbol{\alpha}^T\boldsymbol{x}^{(n)}]^{2}
\end{align}

\begin{align}
\sigma = \frac{1}{N}\sum^N_{n=1}[y^{n} – \boldsymbol{\alpha}^T\boldsymbol{x}^{(n)}]^{2}
\end{align}

次数rの決定

モデルの複雑性がハイパーパラメータである次数rで定まります。クロスバリデーションなどの手法でrを決めるのは時間的な順序があるため使えないのでこの複雑さを決める基準としてAIC(赤池情報量基準)を用います。

\begin{align}
\ln|\sigma| + \frac{2r}{N}
\end{align}

上記の数式が最小になるようにモデルを選択します。

Pythonによる実装

自己回帰モデルを実装します
今回は統計的なモデルが利用できるライブラリである”StatsModels”を使用してモデルを実現します。

データはオーストラリアの公開されている気温のデータを使用しました。ただしデータの中に?などが混じっており数値データ以外の項目があるので適宜削除してください。


#!/usr/bin/env python3


import pylab
from statsmodels.tsa.ar_model import AR
import numpy as np
from pandas import Series

# Reference
#   data : https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/master/csv/texmex/nidd.csv  # noqa


def train_autoregression(train_data, test_data):
    # train autoregression
    print(train_data)
    model = AR(train_data)
    model_fit = model.fit()
    window = model_fit.k_ar
    coef = model_fit.params
    print('Lag: %s' % model_fit.k_ar)
    print('Coefficients: %s' % model_fit.params)
    # make predictions
    print('aic {}'.format(model_fit.aic))
    history = train_data[len(train_data) - window:]
    history = [h for h in history]
    predictions = []

    for index, test in enumerate(test_data):
        length = len(history)
        lag = [history[i] for i in range(length - window, length)]
        yhat = coef[0]
        for d in range(window):
            yhat += coef[d+1] * lag[window - d - 1]
        obs = test
        predictions.append(yhat)
        history.append(obs)

    return predictions


def main():
    wave_data = Series.from_csv('../data/raw/daily-minimum-temperatures-in-me.csv', header=0)  # noqa
    X = np.array(wave_data.values, dtype=float)
    print(len(X))
    train_data_size = 3000
    train_data, test_data = X[1:train_data_size], X[train_data_size:]
    predict = train_autoregression(train_data, test_data)

    pylab.figure(figsize=(14, 8))
    pylab.subplot(211)
    pylab.xlabel('time')
    pylab.ylabel('temperature')
    pylab.plot(test_data, label='real')
    pylab.plot(predict, '--', label='predict')
    pylab.legend(loc='upper right')
    pylab.show()


if __name__ == '__main__':
    main()



コードの下記の部分で自己回帰モデルの出力を実現しています。


       for d in range(window):
           yhat += coef[d+1] * lag[window - d - 1]



StatsModelsのライブラリがAICに基づく最適な次数を設定してくれています。
オレンジが予測波形、青が実波形で波形を模倣できていることが確認できます。

異常度の計算

ホテリング理論を用いて下記のように定義します。
\begin{align}
a(\epsilon^{(t)}) = \frac{1}{\sigma^2}[\epsilon^{(t)} - \sum^r_{l=1}{\alpha_{l}\epsilon^{(t-l)}}]^{2}
\end{align}

αlεt-ll=1r

の部分は自己回帰モデルで導出した値ですが、この部分は他の予測モデルで置き換え可能です。

閾値の設定

前回の手法と同様に閾値は学習データの分位点を基準にするかもしくは自由度1のカイ2乗分布を当てはめて導出します。

ここで異常度は各点において計算されるので異常と判定する場合に何度、閾値を超えたかも設定する必要があります。
その代わり各点における異常度が計算できているため変化点の検知も可能になっています。

自己回帰モデルの派生

詳細は記述しませんが今回紹介したモデルの派生系として下記のようなモデルがあります。

  • MAモデル(移動平均モデル)
  • ARMAモデル(自己回帰移動平均モデル)
  • ARIMAモデル(自己回帰和文移動平均モデル)
  • SARIMAモデル(季節自己回帰和文移動平均モデル)

異常度の導出、閾値の導出は自己回帰モデルと同様なため最初に拡張性が高いと記述したのはそのためです。

状態空間モデル

定常的な時系列データを扱うのに便利な自己回帰モデルを紹介しました。ここからは状態が変化する時系列を扱える状態空間モデルについて紹介します
イメージを掴みやすい動画を見つけたのでこちらをまず紹介します。

見えない状態空間の影響を受けてなおかつそれが時系列データの場合に有効な手法になります。正常な状態だが微妙に状態が変化する事象は現実的な問題として多数存在します。そのようなケースにおいて有効です。


上記の図のように状態空間と観測空間に分けます。このようなモデルを状態空間モデルと呼びます。
下記の数式で記述できます。

\begin{align}
p(\boldsymbol{x}^{t}|\boldsymbol{z}^{t}) = \mathcal{N}(\boldsymbol{x}^{t}| \boldsymbol{C}\boldsymbol{z}^{t}, \boldsymbol{R})
\end{align}

\begin{align}
p(\boldsymbol{z}^{t}|\boldsymbol{z}^{t-1}) = \mathcal{N}(\boldsymbol{z}^{t}| \boldsymbol{A}\boldsymbol{z}^{t - 1}, \boldsymbol{Q})
\end{align}

  • A:観測空間を表現する行列
  • C:状態空間を表現する行列
  • R:観測のばらつき
  • Q:内部状態のばらつき

本記事ではこれらの未知のパラメータの導出には言及しません。入門 機械学習による異常検知―Rによる実践ガイド-井手-剛のp.214からp.220に導出があります。

状態空間モデル

観測されたデータから状態変数の確率分布を求めます。

\begin{align}
p(\boldsymbol{z}^{t}|\boldsymbol X_t) = \mathcal{N}(\boldsymbol{z}^{t}| \boldsymbol{\mu}_{t}, \boldsymbol V_t)
\end{align}

下記の結果用います。本記事ではこれらの未知のパラメータの導出には言及しません。入門 機械学習による異常検知―Rによる実践ガイド-井手-剛のp.221からp.224に導出があります。

t=0の時の状態変数の平均を

\begin{align}
\boldsymbol{A} \boldsymbol \mu_{0} = \boldsymbol z_{0}
\end{align}

として定義し、t=1,...Tに対して次の計算を繰り返します。

\begin{align}
\boldsymbol K_t = \boldsymbol Q_{t-1}\boldsymbol{C}^{T}(\boldsymbol{R} + \boldsymbol{C}\boldsymbol Q_{t-1}\boldsymbol{C}^T)^{-1}
\end{align}

\begin{align}
\boldsymbol \mu_t = \boldsymbol{A}\mu_{t-1}+ \boldsymbol K_{t}(\boldsymbol{x}^{(t)} - \boldsymbol{C}\boldsymbol{A}\mu_{t-1})
\end{align}

\begin{align}
\boldsymbol V_t = (\boldsymbol I_m - \boldsymbol K_{t}\boldsymbol{C})\boldsymbol Q_{t-1}
\end{align}

\begin{align}
\boldsymbol Q_t = \boldsymbol{Q} + \boldsymbol{A}\boldsymbol V_t\boldsymbol{A}^{T}
\end{align}

この時のKtはカルマンゲインと呼ばれ観測値とのばらつきに対する反応の敏感さを表現しています。
状態変数における平均と分散が導出できます。

ここでt-1までの観測データが与えられた時のxの分布を与えます。tの時の予測分布を考えます。多変数正規分布のベイズ公式より

\begin{align}
p(\boldsymbol{x}^{(t)} | \boldsymbol X_{t-1}) = \mathcal{N}(\boldsymbol{x}^{(t)} | \boldsymbol{C}\boldsymbol{A}\boldsymbol \mu_{t-1}, \boldsymbol{R} + \boldsymbol{C}\boldsymbol Q_{t-1}\boldsymbol{C}^{T})
\end{align}

QAは状態変数のパラメメータなので状態変数のパラメータが影響して正規分布のモデルが変化しています。さらに平均は一つ前の平均、分散は一つ前の状態変数の分散に影響して決定しています。
下記の図のように正規分布が変化していくイメージです。


実装

状態空間モデルの実装はPyFluxというライブラリを使用します。
これによって下記に示すわずかなコードで実現が可能ですが再現しているのはローカルレベルモデルのため若干状態空間モデルと異なります。


#!/usr/bin/env python3

import numpy as np
from pandas import Series
import pyflux as pf
import matplotlib.pyplot as plt

# Reference
#   data : https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/master/csv/texmex/nidd.csv  # noqa


def train_autoregression(train_data, test_data):
    # train autoregression
    print(train_data)
    model = pf.LLEV(data=train_data, target='Nile')
    x = model.fit()
    x.summary()
    model.plot_fit(figsize=(15, 10))
    plt.figure(figsize=(15, 5))
    for i in range(10):
        model_index = np.array(model.index)
        plt.plot(model_index, model.simulation_smoother(
            model.latent_variables.get_z_values())[1][0:model_index.shape[0]])
    plt.show()


def main():
    wave_data = Series.from_csv('../data/raw/daily-minimum-temperatures-in-me.csv', header=0)  # noqa
    X = np.array(wave_data.values, dtype=float)
    train_data_size = 3000
    train_data, test_data = X[1:train_data_size], X[train_data_size:]
    train_autoregression(train_data, test_data)


if __name__ == '__main__':
    main()



下記に上から実データに対するモデルのフィット度合いの確認、予測されたモデルとそのモデルが95%の確率で取りうる範囲、データから観測されたノイズを確認できます。


異常スコアの計算

正規分布による異常度を計算すれば良いので下記のような式になります。

\begin{align}
a(\boldsymbol{x}^{(t)}) = (\boldsymbol{x}^{(t)} - \boldsymbol{C}\boldsymbol{A}\boldsymbol \mu_{t-1})\boldsymbol{\sigma}^-1_{t}(\boldsymbol{x}^{(t)} - \boldsymbol{C}\boldsymbol{A}\boldsymbol \mu_{t-1})
\end{align}

\begin{align}
\boldsymbol{\sigma_t} \equiv \boldsymbol{R}+\boldsymbol{C}\boldsymbol Q_{t-1}\boldsymbol{C}^{T}
\end{align}

閾値の設定

検証データで設定するか分位点による設定になります。
状態空間モデルも異常度は各点において計算されるので異常と判定する場合に何度、閾値を超えたかも設定する必要があります。
その代わり各点における異常度が計算できているため変化点の検知も可能になっています。

最後に

弊社では異常検知以外にも物体検出、3次元データ検索エンジンの開発をしています。これらに興味があるエンジニアがいらっしゃれば絶賛採用中なので是非、弊社へ応募してください。

参考

時系列解析_理論編
入門 機械学習による異常検知―Rによる実践ガイド
時系列データに対する異常検知
Autoregression Models for Time Series Forecasting With Python
Understanding Kalman Filters, Part 4: Optimal State Estimator Algorithm
Gaussian State Space Models

その他の記事

Other Articles

2019/11/06
GraphQLの入門書を翻訳しました

2019/09/20
Kabuku Connect 即時見積機能のバックエンド開発

2019/08/14
Maker Faire Tokyo 2019でARゲームを出展しました

2019/07/25
夏休みだョ!WebAssembly Proposal全員集合!!

2019/07/08
鵜呑みにしないで! —— 書籍『クリーンアーキテクチャ』所感 ≪null 篇≫

2019/07/03
W3C Workshop on Web Games参加レポート

2019/06/28
TypeScriptでObject.assign()に正しい型をつける

2019/06/25
カブクエンジニア開発合宿に行ってきました 2019夏

2019/06/21
Hola! KubeCon Europe 2019の参加レポート

2019/06/19
Clean Resume きれいな環境できれいな履歴書を作成する

2019/05/20
[Web フロントエンド] 状態更新ロジックをフレームワークから独立させる

2019/04/16
C++のenable_shared_from_thisを使う

2019/04/12
OpenAPI 3 ファーストな Web アプリケーション開発(Python で API 編)

2019/04/08
WebGLでレイマーチングを使ったCSGを実現する

2019/04/02
『エンジニア採用最前線』に感化されて2週間でエンジニア主導の求人票更新フローを構築した話

2019/03/29
その1 Jetson TX2でk3s(枯山水)を動かしてみた

2019/03/27
任意のブラウザ上でJestで書いたテストを実行する

2019/02/08
TypeScript で “radian” と “degree” を間違えないようにする

2019/02/05
Python3でGoogle Cloud ML Engineをローカルで動作する方法

2019/01/18
SIGGRAPH Asia 2018 参加レポート

2019/01/08
お正月だョ!ECMAScript Proposal全員集合!!

2019/01/08
カブクエンジニア開発合宿に行ってきました 2018秋

2018/12/25
OpenAPI 3 ファーストな Web アプリケーション開発(環境編)

2018/12/23
いまMLKitカスタムモデル(TF Lite)は使えるのか

2018/12/21
[IoT] Docker on JetsonでMQTTを使ってCloud IoT Coreと通信する

2018/12/11
TypeScriptで実現する型安全な多言語対応(Angularを例に)

2018/12/05
GASでCompute Engineの時間に応じた自動停止/起動ツールを作成する 〜GASで簡単に好きなGoogle APIを叩く方法〜

2018/12/02
single quotes な Black を vendoring して packaging

2018/11/14
3次元データに2次元データの深層学習の技術(Inception V3, ResNet)を適用

2018/11/04
Node Knockout 2018 に参戦しました

2018/10/24
SIGGRAPH 2018参加レポート-後編(VR/AR)

2018/10/11
Angular 4アプリケーションをAngular 6に移行する

2018/10/05
SIGGRAPH 2018参加レポート-特別編(VR@50)

2018/10/03
Three.jsでVRしたい

2018/10/02
SIGGRAPH 2018参加レポート-前編

2018/09/27
ズーム可能なSVGを実装する方法の解説

2018/09/25
Kerasを用いた複数入力モデル精度向上のためのTips

2018/09/21
競技プログラミングの勉強会を開催している話

2018/09/19
Ladder Netwoksによる半教師あり学習

2018/08/10
「Maker Faire Tokyo 2018」に出展しました

2018/08/02
Kerasを用いた複数時系列データを1つの深層学習モデルで学習させる方法

2018/07/26
Apollo GraphQLでWebサービスを開発してわかったこと

2018/07/19
【深層学習】時系列データに対する1次元畳み込み層の出力を可視化

2018/07/11
きたない requirements.txt から Pipenv への移行

2018/06/26
CSS Houdiniを味見する

2018/06/25
不確実性を考慮した時系列データ予測

2018/06/20
Google Colaboratory を自分のマシンで走らせる

2018/06/18
Go言語でWebAssembly

2018/06/15
カブクエンジニア開発合宿に行ってきました 2018春

2018/06/08
2018 年の tree shaking

2018/06/07
隠れマルコフモデル 入門

2018/05/30
DASKによる探索的データ分析(EDA)

2018/05/10
TensorFlowをソースからビルドする方法とその効果

2018/04/23
EGLとOpenGLを使用するコードのビルド方法〜libGLからlibOpenGLへ

2018/04/23
技術書典4にサークル参加してきました

2018/04/13
Python で Cura をバッチ実行するためには

2018/04/04
ARCoreで3Dプリント風エフェクトを実現する〜呪文による積層造形映像制作の舞台裏〜

2018/04/02
深層学習を用いた時系列データにおける異常検知

2018/04/01
音声ユーザーインターフェースを用いた新方式積層造形装置の提案

2018/03/31
Container builderでコンテナイメージをBuildしてSlackで結果を受け取る開発スタイルが捗る

2018/03/23
ngUpgrade を使って AngularJS から Angular に移行

2018/03/14
Three.jsのパフォーマンスTips

2018/02/14
C++17の新機能を試す〜その1「3次元版hypot」

2018/01/11
異常検知の基礎

2018/01/09
three.ar.jsを使ったスマホAR入門

2017/12/17
Python OpenAPIライブラリ bravado-core の発展的な使い方

2017/12/15
WebAssembly(wat)を手書きする

2017/12/14
AngularJS を Angular に移行: ng-annotate 相当の機能を TypeScrpt ファイルに適用

2017/12/08
Android Thingsで4足ロボットを作る ~ Android ThingsとPCA9685でサーボ制御)

2017/12/06
Raspberry PIとDialogflow & Google Cloud Platformを利用した、3Dプリンターボット(仮)の開発 (概要編)

2017/11/20
カブクエンジニア開発合宿に行ってきました 2017秋

2017/10/19
Android Thingsを使って3Dプリント戦車を作ろう ① ハードウェア準備編

2017/10/13
第2回 魁!! GPUクラスタ on GKE ~PodからGPUを使う編~

2017/10/05
第1回 魁!! GPUクラスタ on GKE ~GPUクラスタ構築編~

2017/09/13
「Maker Faire Tokyo 2017」に出展しました。

2017/09/11
PyConJP2017に参加しました

2017/09/08
bravado-coreによるOpenAPIを利用したPythonアプリケーション開発

2017/08/23
OpenAPIのご紹介

2017/08/18
EuroPython2017で2名登壇しました。

2017/07/26
3DプリンターでLチカ

2017/07/03
Three.js r86で何が変わったのか

2017/06/21
3次元データへの深層学習の適用

2017/06/01
カブクエンジニア開発合宿に行ってきました 2017春

2017/05/08
Three.js r85で何が変わったのか

2017/04/10
GCPのGPUインスタンスでレンダリングを高速化

2017/02/07
Three.js r84で何が変わったのか

2017/01/27
Google App EngineのFlexible EnvironmentにTmpfsを導入する

2016/12/21
Three.js r83で何が変わったのか

2016/12/02
Three.jsでのクリッピング平面の利用

2016/11/08
Three.js r82で何が変わったのか

2016/12/17
SIGGRAPH 2016 レポート

2016/11/02
カブクエンジニア開発合宿に行ってきました 2016秋

2016/10/28
PyConJP2016 行きました

2016/10/17
EuroPython2016で登壇しました

2016/10/13
Angular 2.0.0ファイナルへのアップグレード

2016/10/04
Three.js r81で何が変わったのか

2016/09/14
カブクのエンジニアインターンシッププログラムについての詩

2016/09/05
カブクのエンジニアインターンとして3ヶ月でやった事 〜高橋知成の場合〜

2016/08/30
Three.js r80で何が変わったのか

2016/07/15
Three.js r79で何が変わったのか

2016/06/02
Vulkanを試してみた

2016/05/20
MakerGoの作り方

2016/05/08
TensorFlow on DockerでGPUを使えるようにする方法

2016/04/27
Blenderの3DデータをMinecraftに送りこむ

2016/04/20
Tensorflowを使ったDeep LearningにおけるGPU性能調査

→
←

関連職種

Recruit

フロントエンドエンジニア(TypeScript)

業務内容

当ポジションは弊社Webサービスのフロントエンド機能設計及び実装を担当します。 サービス毎の開発チームで2週間スプリントのスクラム開発を実施しています。 週次で開発チームミーティングを実施し、実装設計の相談や工数見積もりを行います。 全ての開発コードはレビューと自動テストによって品質を保っています。 また、リファクタリングやフレームワークのバージョンアップも開発フローに組込み、技術的負債を放置しない開発を目指しています。

インターン(Webエンジニア)

業務内容

業務から独立した、調査・研究系のタスクをおまかせしています。コードレビュー、 社内での報告会、 ブログ記事執筆を通して着実にスキルアップしていただくことを目指しています。 (希望があれば、プロダクトの開発業務もおまかせします。)

→
←

お客様のご要望に「Kabuku」はお応えいたします。
ぜひお気軽にご相談ください。

お電話でも受け付けております
03-6380-2750
営業時間:09:30~18:00
※土日祝は除く