EeePCの軌跡: 2023年7月

« 2023年6月 | トップページ | 2023年8月 »

2023年7月30日 (日)

ユニバーサルスタジオジャパン(USJ)へ行ってきました

単なる旅行記ですが、地味に活躍したガジェットの紹介もあり。
表題の通り、USJへ行ってきました。

会社の福利厚生でもらえる旅行ポイントなるものがあって、それが結構溜まっていることもあり、少し豪華？な旅でした。
なお、乗ったのは7時11分発のこだま783号の

グリーン車です。
こだまのグリーン車なんて、利用するやついるんかいな？と思っておりましたが、はい、思いっきり利用しました。
結構ガラガラですが、ちょくちょく人がいました。
名古屋から新大阪なので、4駅ほど。こだまでも、さほどかかりません。

前回、乗ったのが普通車(東京出張時)でしたが、この広さを知ってしまうと、快適過ぎて癖になります。

で、朝食はこちら。朝早く出たため、経路上のコンビニがどこも開いてなかったため、新幹線構内のキヨスクで購入。

この途上、ちょっと珍しいエピソードがありまして、米原駅到着直前というかもう駅のホームにどっぷり入っている状態で、ドアも開かず停車。そのまま、3分ほどが経過。
何があったのかはよく分かりませんが、安全点検のため停車しているとのアナウンスが入ります。
そういえば以前、名古屋駅構内で車体の異常を感じて緊急点検をしたら、車台に亀裂が入ってた事件がありましたね。こうして何らかの予兆を感知したら、安全点検を実施しているみたいです。

そういえばですが、この間、東京出張した時には流れていた車内チャイム「Ambitious Japan」が、

「会いに行こう」に変わってました。

悪い曲ではありませんが、長年、親しんできた「Ambitious Japan」だけに、ちょっと寂しい気分でしたね。

てことがありましたが、予定より3分ほど遅れて新大阪に到着しました。

途中をすっぱかしますが、なんやかんやあってUSJに到着。
ご覧の通り、すごいいい天気でした。

ちなみに、例のNikonの40倍カメラ「COOLPIX B500」を持っていったため、

最大望遠で、御覧のような写真も撮影可能です。
案外、この望遠がちょくちょく活躍しました。

例えば、USJのスーパー・ニンテンドー・ワールドですが、あれって今、入場制限がかかってるんですよね。
で、到着して入場券のQRコードをアプリで読み込ませ、入場資格の抽選をやるようなのですが、私が着いた時点ではすでに18時以降しか空いてなくて、行くのを断念しました。

それじゃあ、あまりにも口惜しいということで、

この40倍ズームカメラで、その一部を写真に収めてきました。
なお、これを撮影した場所はスパイダーマンとミニオンズのエリアの中間ほどのところにある池のほとり。
ほぼ対角上のエリアにある建物を、この通り捉えることができました。

他にも、メルズ・ドライブインというところで食事をしていたのですが、いきなりパレードが始まり、それをこの40倍カメラで撮影。
結構店の奥にいたので、スマホじゃこんな距離からの撮影は無理でしたね。

ちなみに、その時食べたのはこのBBQベーコンチーズバーガーセット。このセットで、お値段は1,800円。うーん、観光地価格ですね。

高いだけあって、ハンバーガー本体はかなりボリューミーでした。

さて、肝心の園内ですが、なんていうか、ユニバーサルスタジオらしさが爆発・・・

と思いきや、日本製コンテンツが結構幅を利かせてました。
なお、ここはスパイダーマンのエリアの真ん前にある、スーパーマリオ・パワーアップ・サマーのショー会場。
スーパー・ニンテンドー・ワールドには入れませんでしたが、ところどころこんな具合にマリオがちりばめられています。

なお、この場所には2時より行われるショーにこれだけの人が集まり、

大量の水をばらまいてました。結構、濡れます。濡れてもすぐに乾くものの、顔にかかると汗が一緒に流れるのか、かなり目に沁みます。タオルは必須。

他にも、ちょっと気になったのはこれ。
ジュラシックパークなエリアに、なぜか「ワンピース」です。このカフェの真ん中にある恐竜の骨っぽいものが飾られたエリアの壁画には、ワンピースが描かれてました(写真、右下です)。
うまく溶け込んではいるものの、こう言っては何ですが、世界観が雑多な感じです。
まあ、子供らは大喜びでしたが。

ところでこのUSJ、アメ車の旧車があちこちに置かれてました。ある場所にも唐突に、イエローキャブが置かれてます。

が、この車の窓ガラスには、こんな注意書きが。
博物館レベルって・・・いや、そんな大事な車、こんなところに唐突に置かないで。トヨタ博物館にでも寄贈して。

他にも、この巨大なジョーズのモニュメント？があったのですが、

その傍らには、こんな車が。こっちは、特に張り紙もなし。
先ほど昼食を食べたメルズ・ドライブインというお店の周辺にも、旧型のアメ車がずらりと並んでました。ちょっと写真にちゃんと収めてませんでしたが。

実はいろいろとトラブル続きでして、まず妻が段差に足を取られて転んでひざを派手に擦りむいて、園内の診療所で応急処置の後、外にあるマツキヨに消毒液を買いに行ったり、私自身も暑さにやられて急に気分が悪くなったりと、なんだかいつもとは違う感じでした。行きの新幹線の3分遅れは、まさにトラブルの予兆だったのでしょうか？

それにしても、ユニバーサルスタジオといえば、私にとっては「バック・トゥ・ザ・フィーチャー」なんですけど、音楽は流れたものの、それに関するエリアもグッズもありませんでした。えっ、なんでバック・トゥ・ザ・フィーチャーがないの？

などなど、いろいろあったUSJを後にします。

それにしても、名古屋と東京はかなり慣れましたが、大阪は全然慣れませんね。特に大阪と新大阪の間を行き来する電車は、相変わらず分かりづらい。

なお、帰りの新幹線はひかりで、同じくグリーン車でした。

戦利品です。

色々買ってます。カールは6個ぐらい買いました。

そしてこの強烈なApple Watch日焼け。この日は19000歩近く歩きました。疲れましたねぇ・・・

そんな我々家族を支えてくれたのは、これらのガジェット。
Nikon COOLPIX B500、2台のモバイルバッテリーとケーブル、そしてXPeria ACE III(なぜかくしゃくしゃのレシートと一緒に映ってます)。
XPeriaは、メインスマホであるiPhone 12の電池消耗を減らすため、マップ表示や電車の乗り換え案内などを担わせました。おかげで私のiPhone 12は最後までバッテリーが持続。
一方で家族のiPhoneは軒並み途中でバッテリーを消耗し、これで充電してました。

始めて行ったUSJでしたが、もう行くことはあるのかどうか。

ＵＳＪを劇的に変えた、たった１つの考え方　成功を引き寄せるマーケティング入門 (角川書店単行本)

固定リンク | 0 | コメント (2)

2023年7月27日 (木)

在宅ワークで大容量バッテリー「JVC-BN-RB10」を使ったら思わぬ落とし穴が

表題の通りですが、先日買った1002Whの大容量バッテリー「JVC-BN-RB10」を使いました。

椅子で作業していたら首が痛くなったので、Yogiboで仕事しようと思い立ちました。

会社のノートPCですが、これにリモート会議用のスピーカーマイクを引っ付けてます。
なお、手が腱鞘炎で痛いため、マウスではなく敢えてトラックパッドのみで使います。

ノートPC自体を使うことは、以前の記事「JVC BN-RB10CでメインPCを動かしてみた: EeePCの軌跡」でもやってるので、さほど目新しいネタではない・・・
かと思いきや思わぬ落とし穴があったので、ご報告。

このバッテリーにはUSB-C端子があるので、当初はここから給電しておりました。
が、1時間ほどリモート会議をやったら、ノートPC側のバッテリーが減少していることが発覚。

あれ、でかいバッテリーにつないでいるんだが、何で減るの？と思いつつ、端子をよく見ると、

えっ、このUSB-C端子、18Wまでしか対応してないの！？
そうです、大容量バッテリーだから、てっきり40～60W出るものだと思ってました。
が、このUSB-C端子を使う限りでは、うちの10000mAh(37Wh)のバッテリーが出せる出力と何ら変わらないことが判明。

会社のPCは18Wでも動かないことはないんですが、リモート会議時の出力だと18Wでは足りないため、本体バッテリーを消耗してしまうようです。あらら・・・

てことで、カッコ悪いですが、御覧の通りAC電源の方から取ることにしました。

これにつないだ途端、一気に出力は50～70Wほどまで上昇。

と思いきや、時々1Wまで下がります。これは本体の仕様でしょうかね。
ともかく、消耗したバッテリーをあっという間に充電してくれました。

どうやらこのバッテリーのUSB-C端子は、スマホ/タブレットのみを想定している模様です。ノートPCまでは考えていない。というか、ノートPCはACを使え、ということみたいです。

大きいバッテリーだからといって、USB-C端子が必ずしも高出力というわけではない、というお話でした。

JVCケンウッドポータブル電源 BN-RB10-CA 充電池容量 278,400mAh/1,002Wh

固定リンク | 0 | コメント (0)

2023年7月26日 (水)

MT法で音声の異常検知をやってみた

久しぶりに「機械学習」的なネタです。今回、MT法というのを使いました。

「MT法とは何ぞや？」という方は、以下を参照。

MT法で異常検知（python） - Qiita

こっちの動画もおすすめです。

要するにですね、正規分布に従っているデータからの外れっぷりを示す「マハラノビス距離」を使って、そのデータの異常度合を測るという方法で、以下のように使います。

(１) まず、正常データのみをそろえて、そのデータの平均値、分散を計算する。
(２) その平均、分散を使って、未知のデータのマハラノビス距離を計算する。

すると、異常なデータほど、このマハラノビス距離ってやつが大きくなるため、その数値を見て異常発生キターッ！と分かるというわけです。

この手法の特徴としては、正常データの分布を元に予測するため、大量の「正常データ」からごくわずかな「異常データ」の出現を測ることができます。
現実問題として、正常がたくさんで異常が滅多に出ないという製品が多いですから、そういう現場では重宝されます。

なので、まともそうな業者の保険修理件数があれば、ビッ〇モーターの保険修理の異常っぷりも検出できた・・・
いや、それは異常検知以前の問題ですね。

さて、今回はこのMT法を使って「音」の異常を検知してみます。

使った音は、うちの充電式シェーバーで髭を剃る時の音、です。

この何の変哲もない髭剃り機の何の異常を検知するのか？

この充電式のシェーバー、充電したてのころは元気がいいんですが、二週間近くなるとだんだんと電池が切れかかってきたのか、音が明らかに弱ってきます。といっても、毎日聞いていると気づかないのですが、充電した直後に明らかに音が変わるので、すぐに分かります。

つまり、徐々にではあるんですが日々、音が弱くなってきているんだろうな、と。

その衰えっぷりを、このマハラノビス距離ってやつで表現できないかと考えた次第。

てことで、今回用意したのは、充電した直後の１日目のデータ(hige01.wav)と、充電切れかけヤバヤバの12日目のデータ(hige12.wav)を使いました。どちらもだいたい50秒程度のデータ。

で、どうやってこの異常を検知しようか？

方針を決めるためにはまず、データを可視化してみましょう。

まずは充電ほやほやの元気のいいシェーバーの音。
周波数は0～5000Hzに絞ってます。
だいたい1500～3000Hzに音の中心が集まっているのが分かります。

一方こちらは充電切れ寸前のシェーバーの音データ。

明らかに音圧が減ってますが、ところどころ2000～3000Hzの音がガクッと下がってますね。

ひげを剃ってる間、負荷によってモーターの回転数が下がるみたいで、それがこのまだらな模様の原因となっているものと推測されます。

てことで、なんとなく分かったのは、「1500～3000Hz辺りの音データから特徴を取り出してやれば、シェーバーの充電時期が検知できるんじゃないか？」ということです。

ということで、以下のような方針を立てました。

(1) 1500～3000Hzの音を、250Hz刻みで平均化して出力
(2) その際、1秒ごとに音を分類できるように、IDをつけておく

上のスペクトログラムを各プログラムと合わせて、この(1)、(2)を実行するプログラムを以下に載せておきます。

import os
import sys
import numpy as np
import librosa
import matplotlib.pyplot as plt
import librosa.display
import pandas as pd 
#音声ファイル読み込み
args = sys.argv
wav_filename = "hige01.wav"

data, rate = librosa.load(wav_filename, sr=None, mono=False)
# フレーム長
fft_size = 2048
# フレームシフト長 
hop_length = int(fft_size * 0.2)
# 時間の刻み幅 (s)
step_time = hop_length / rate
# 周波数の刻み幅 (Hz)
step_freq = rate / fft_size
# 短時間フーリエ変換実行
amplitude = np.abs(librosa.core.stft(data, n_fft=fft_size, hop_length=hop_length))
# 振幅をデシベル単位に変換
log_power = librosa.core.amplitude_to_db(amplitude)
fig_filename = os.path.splitext(os.path.basename(wav_filename))[0]
fft_spect = log_power.T
print(fft_spect.shape)
time_ax = []
freq_ax = []
for i in range(fft_spect.shape[0]):
    time_ax.append(step_time * i)
for j in range(fft_spect.shape[1]):
    freq_ax.append(step_freq * j)
time_ax = np.array(time_ax)
freq_ax = np.array(freq_ax)
df_spect = pd.DataFrame(fft_spect, columns=freq_ax)
df_time = pd.Series(time_ax,name="Time[s]")
df_spect_time = pd.concat([df_time,df_spect],axis=1)
df_spect_time.to_csv(fig_filename + '_pd.csv')
# グラフ表示
librosa.display.specshow(log_power, sr=rate, hop_length=hop_length, x_axis='time', y_axis='hz', cmap='jet')
plt.ylim(0,4000)
plt.xlim(0,10)
plt.colorbar(format='%+2.0f dB')  
plt.title('スペクトログラム' ,fontname="MS Gothic")
plt.savefig("10s_4000Hz_" + fig_filename)
#plt.show()
# データ加工
# 1500～3000Hzを250Hzづつ足して平均化し出力する
start_freq = 1500
interval_freq = 250
wave_feat = []

for k in range(fft_spect.shape[0]):
    tmp_wave = 0
    for l in range(int(interval_freq / step_freq)):
        tmp_wave += fft_spect[k, int(start_freq / step_freq) + l]
        tmp_wave = tmp_wave / (500 / step_freq)
    wave_feat.append(tmp_wave)

start_freq1 = 1750
wave_feat1 = []

for k in range(fft_spect.shape[0]):
    tmp_wave = 0
    for l in range(int(interval_freq / step_freq)):
        tmp_wave += fft_spect[k, int(start_freq1 / step_freq) + l]
        tmp_wave = tmp_wave / (500 / step_freq)
    wave_feat1.append(tmp_wave)

start_freq2 = 2000
wave_feat2 = []

for k in range(fft_spect.shape[0]):
    tmp_wave = 0
    for l in range(int(interval_freq / step_freq)):
        tmp_wave += fft_spect[k, int(start_freq2 / step_freq) + l]
        tmp_wave = tmp_wave / (500 / step_freq)
    wave_feat2.append(tmp_wave)

start_freq3 = 2250
wave_feat3 = []

for k in range(fft_spect.shape[0]):
    tmp_wave = 0
    for l in range(int(interval_freq / step_freq)):
        tmp_wave += fft_spect[k, int(start_freq2 / step_freq) + l]
        tmp_wave = tmp_wave / (500 / step_freq)
    wave_feat3.append(tmp_wave)

start_freq4 = 2500
wave_feat4 = []

for k in range(fft_spect.shape[0]):
    tmp_wave = 0
    for l in range(int(interval_freq / step_freq)):
        tmp_wave += fft_spect[k, int(start_freq4 / step_freq) + l]
        tmp_wave = tmp_wave / (500 / step_freq)
    wave_feat4.append(tmp_wave)

start_freq5 = 2750
wave_feat5 = []

for k in range(fft_spect.shape[0]):
    tmp_wave = 0
    for l in range(int(interval_freq / step_freq)):
        tmp_wave += fft_spect[k, int(start_freq5 / step_freq) + l]
        tmp_wave = tmp_wave / (500 / step_freq)
    wave_feat5.append(tmp_wave)

time_id_step = 100 # IDをつけるためのID刻み幅
time_id = []

for i in range(time_ax.shape[0]):
    time_id.append(i // time_id_step)
wave_feat_all = np.array([time_ax, wave_feat,wave_feat1,
　　　　　　　　　wave_feat2,wave_feat3,wave_feat4,wave_feat5]).T

df_wave_feat = pd.DataFrame(wave_feat_all, columns=["time",start_freq,start_freq1,
                            start_freq2,start_freq3,start_freq4,start_freq5])
df_id = pd.Series(time_id,name="id")
df_wave_time = pd.concat([df_id,df_wave_feat],axis=1)
# IDで出力範囲指定する
df_wave_time = df_wave_time.query('0< id < 54')
df_wave_time.to_csv(fig_filename + '_feat.csv')
#
df_wave_feat.plot(subplots=True, sharex=True, figsize=(6,12))
plt.savefig("1500-3000Hz_wave_feat_" + fig_filename)

そういえばですが、予めnumpy、matplotlib、pandas、scikit-learn、scipy、seabornなどをpipコマンドでインストールしておいてください。
そしてこのプログラムでは「librosa」というライブラリを使ってます。
上のプログラム、ちょっとやっつけな部分がありますがご容赦ください。このプログラム名を「sound_spect.py」としておきます。
で、これをまず実行。

> python sound_spect.py

これを、10行目の「wave_filename = 」の後ろを「hige01.wav」「hige12.wav」と変えてそれぞれ実行すると、「hige01_feat.csv」「hige12_feat.csv」が出力されます。
なお、ファイルごとに録音時間がバラバラなので、このプログラムでは1～53秒までの音を抜き出すようにしております。

こんな感じのファイルが作成されます。なお、列名が1500となっているのは、1500Hzから250Hz分の平均値、という意味です。また2列目のIDは、1～53まで変わります。

さて、この時系列データの特徴のとり方というのもいろいろな方法があるんですが、ここではオーソドックスに、

・ 1秒ごとに(つまり同じIDごとに)、1500～2750Hz帯それぞれの平均、最大、最小、標準偏差を計算する

という方法を取ります。

で、それをやってくれるプログラムが、以下。

import pandas as pd
import numpy as np
def process_timeseries_data(input_file, output_file):
    # CSVファイルを読み込む
    df = pd.read_csv(input_file)
    # idごとにグループ化し、time、No.1、No.2、No.3の平均、最大、最小、標準偏差を計算
    grouped = df.groupby('id').agg({
        'time': 'first',   # 各グループ内の最初のtimeを取得
        '1500': ['mean', 'max', 'min', 'std'],  # No.1の平均、最大、最小、標準偏差を計算
        '1750': ['mean', 'max', 'min', 'std'],  # No.2の平均、最大、最小、標準偏差を計算
        '2000': ['mean', 'max', 'min', 'std'],   # No.3の平均、最大、最小、標準偏差を計算
        '2250': ['mean', 'max', 'min', 'std'],   # No.4の平均、最大、最小、標準偏差を計算
        '2500': ['mean', 'max', 'min', 'std'],   # No.5の平均、最大、最小、標準偏差を計算
        '2750': ['mean', 'max', 'min', 'std'],   # No.6の平均、最大、最小、標準偏差を計算
    })
    # 列名を整理
    grouped.columns = ['_'.join(col) for col in grouped.columns]
    # 出力ファイルに書き込む
    grouped.to_csv(output_file, index=True)
if __name__ == "__main__":
    input_file = "hige01_feat.csv"    # 入力ファイル名を適切に指定してください
    output_file = "output_hige01.csv"  # 出力ファイル名を適切に指定してください
    process_timeseries_data(input_file, output_file)

実はこのコード、考えるのがめんどくさくなって全部ChatGPTに作らせました。なんと、一発で動きましたよ。賢いでちゅね～、ChatGPT。
このプログラムは「spect_feat.py」と名付けます。
これを実行します。

> python spect_feat.py

なお、読み込むファイル名は下から３行目の「input_file =」の後ろ、出力ファイル名はその下の「output_file =」を変えます。

で、これを使って2つの音源ファイルから作った「hige01_feat.csv」「hige12_feat.csv」を読み込ませ、「output_hige01.csv」「output_hige12.csv」を得ます。

なお、出力されたファイルの中身はこんな感じ。ID、時間に続いて、1500Hz帯の平均、最大、最小、標準偏差・・・という感じに並んでます。

それぞれの周波数帯に対して、4つの値。これが5つあって、かつ時間が53秒分あるため、全部で20×53のデータになります。
元が1025×6000とかいう巨大なデータなので、これでもかなり圧縮されております。

なお、「output_hige01.csv」のデータを「正常データ」、「output_hige12.csv」を「異常データ」とみなして、これをMT法によってより分けられるかを検証していきます。

まずは、この2つのデータをExcel等を使い「id」「time」の2列を消し、上から01、12の順に結合しておきます。ファイル名は、ここでは「hige01_12.csv」としておきます。

このファイルは、1行目に列名、次の53行は「正常データ」、その後ろの53行は「異常データ」ということになってます。

その上で、次のプログラムを実行します。

# 直交表から任意の因子数についての要因効果図を作る
# 使う直交表はL12～L64まで
# パッケージインポート
import osimport re
import collections
import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt
import math
import sklearn
from sklearn import preprocessing
import scipy
from scipy.spatial import distance
import functools
import seaborn as sns
import sys
import tqdm as tq
from tqdm import tqdm
import gzip
import glob
import datetime as dt
import gc
import time
sns.set()

# =========== ユーザー定義 ==================
# データ読み込み
data0 = pd.read_csv('hige01_02.csv', encoding='cp932') # 51番目からは異常判定のあるデータ
# 正常データの行数 53行ある場合は、0始まりの52行目、という意味で「52」と入れる
s_data = 52
# ==========================================

# データを表示(表、グラフ)
print(data0.head())
data0.plot(figsize=(10, 6))
plt.show()

# マハラノビス距離計算前に、まずデータの正規化しておく(平均値０、標準偏差１のデータにする)
ss = preprocessing.StandardScaler()
data1 = ss.fit_transform(data0)
data = pd.DataFrame(data=data1, columns=data0.columns.values)
# 以下、要因効果図作成のため使用する直交表の一覧(2水準系)を定義
# L12 → 因子数が11以下
orthogonal_array_L12 = [
    [1,1,1,1,1,1,1,1,1,1,1],
    [2,1,2,1,1,1,2,2,2,1,2],
    [2,2,1,2,1,1,1,2,2,2,1],
    [1,2,2,1,2,1,1,1,2,2,2],
    [2,1,2,2,1,2,1,1,1,2,2],
    [2,2,1,2,2,1,2,1,1,1,2],
    [2,2,2,1,2,2,1,2,1,1,1],
    [1,2,2,2,1,2,2,1,2,1,1],
    [1,1,2,2,2,1,2,2,1,2,1],
    [1,1,1,2,2,2,1,2,2,1,2],
    [2,1,1,1,2,2,2,1,2,2,1],
    [1,2,1,1,1,2,2,2,1,2,2]
    ]
# L32 → 因子数が21以下
orthogonal_array_L32 = [
    [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1],
    [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2],
    [1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2],
    [1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1],
    [1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2],
    [1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1],
    [1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1],
    [1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,1,1,1,1,2,2,2,2],
    [1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2],
    [1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1],
    [1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,1,1,2,2,2,2,1,1,2,2,1,1],
    [1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,1,1,2,2,1,1,2,2],
    [1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1],
    [1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2],
    [1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,1,1,2,2],
    [1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,1,1,2,2,2,2,1,1],
    [2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2],
    [2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1],
    [2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,1,2,1,2,2,1,2,1,2,1,2,1],
    [2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,1,2,1,2,1,2,1,2],
    [2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1],
    [2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2],
    [2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,1,2,1,2],
    [2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,1,2,1,2,2,1,2,1],
    [2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1],
    [2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2],
    [2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2,1,2,1,1,2,2,1,1,2],
    [2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2,1],
    [2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2],
    [2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1],
    [2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1],
    [2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,1,2,2,1,2,1,1,2]
    ]
# L64 → 因子数が63以下
orthogonal_array_L64 = [
    [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1],
    [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,
2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2],
    [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2],
    [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1],
    [1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,1,1,
1,1,1,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2],
    [1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,
2,2,2,2,2,2,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1],    [1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1],
    [1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,2,2,
2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2],
    [1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,
1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2],
    [1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,2,2,
2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1],
    [1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,1,1,
1,1,2,2,2,2,1,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1],
    [1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,
2,2,1,1,1,1,2,2,2,2,1,1,1,1,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2],
    [1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,1,1,
1,1,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1],
    [1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,2,2,
2,2,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,2,2,2,2],
    [1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,1,1,1,1,2,2,2,2,1,1,
1,1,2,2,2,2,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,1,1,1,1,2,2,2,2],
    [1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,1,1,1,1,2,2,2,2,2,2,
2,2,1,1,1,1,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1],
    [1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,
2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2],
    [1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,2,2,
1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1],
    [1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,1,1,
2,2,1,1,2,2,1,1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1],
    [1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,
1,1,2,2,1,1,2,2,1,1,2,2,1,1,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2],
    [1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,1,1,
2,2,1,1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,1,1,2,2,2,2,1,1,2,2,1,1],
    [1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,2,2,
1,1,2,2,1,1,1,1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,1,1,2,2],
    [1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,1,1,2,2,1,1,2,2,1,1,
2,2,1,1,2,2,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,1,1,2,2,1,1,2,2],
    [1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,1,1,2,2,1,1,2,2,2,2,
1,1,2,2,1,1,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,2,2,1,1,2,2,1,1],
    [1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,
2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1],
    [1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,2,2,
1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2],
    [1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,1,1,
2,2,2,2,1,1,1,1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2],
    [1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,
1,1,1,1,2,2,2,2,1,1,1,1,2,2,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1],
    [1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,1,1,
2,2,2,2,1,1,2,2,1,1,1,1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,1,1,2,2],
    [1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,2,2,
1,1,1,1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,1,1,2,2,2,2,1,1],
    [1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,1,1,2,2,2,2,1,1,1,1,
2,2,2,2,1,1,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,1,1,2,2,2,2,1,1],
    [1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,1,1,2,2,2,2,1,1,2,2,
1,1,1,1,2,2,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,2,2,1,1,1,1,2,2],
    [2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,
1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2],
    [2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,2,1,
2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1],
    [2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,1,2,
1,2,1,2,1,2,1,2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1],
    [2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,
2,1,2,1,2,1,2,1,2,1,2,1,2,1,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2],
    [2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,1,2,
1,2,1,2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,1,2,1,2,2,1,2,1,2,1,2,1],
    [2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,2,1,
2,1,2,1,2,1,1,2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,1,2,1,2],
    [2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,1,2,1,2,1,2,1,2,1,2,
1,2,1,2,1,2,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,1,2,1,2,1,2,1,2],
    [2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,1,2,1,2,1,2,1,2,2,1,
2,1,2,1,2,1,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,2,1,2,1,2,1,2,1],
    [2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,
1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1],
    [2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,2,1,
2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2],
    [2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,1,2,
1,2,2,1,2,1,1,2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2],
    [2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,
2,1,1,2,1,2,2,1,2,1,1,2,1,2,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1],
    [2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,1,2,
1,2,2,1,2,1,2,1,2,1,1,2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,1,2,1,2],
    [2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,2,1,
2,1,1,2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,1,2,1,2,2,1,2,1],
    [2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,1,2,1,2,2,1,2,1,1,2,
1,2,2,1,2,1,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,1,2,1,2,2,1,2,1],
    [2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,1,2,1,2,2,1,2,1,2,1,
2,1,1,2,1,2,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,2,1,2,1,1,2,1,2],
    [2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,
2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1],
    [2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,2,1,
1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2],
    [2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,1,2,
2,1,1,2,2,1,1,2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2],
   [2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,
1,2,2,1,1,2,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1],
    [2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,1,2,
2,1,1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2,1,2,1,1,2,2,1,1,2],
    [2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,2,1,
1,2,2,1,1,2,1,2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2,1],
    [2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2,1,1,2,
2,1,1,2,2,1,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2,1],
    [2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2,1,2,1,
1,2,2,1,1,2,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,2,1,1,2,2,1,1,2],
    [2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,
2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2],
    [2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,2,1,
1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1],
    [2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,1,2,
2,1,2,1,1,2,1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1],
    [2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,
1,2,1,2,2,1,2,1,1,2,1,2,2,1,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2],
    [2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,1,2,
2,1,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1],
    [2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,2,1,
1,2,1,2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2,1,2,1,1,2],
    [2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,1,2,2,1,2,1,1,2,1,2,
2,1,2,1,1,2,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,1,2,2,1,2,1,1,2],
    [2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,1,2,2,1,2,1,1,2,2,1,
1,2,1,2,2,1,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1]
]

# MDの関数を定義
def get_mean_cov(feat):

    mean = np.mean(feat.astype(float).values, axis=0)
    cov = np.cov(feat.astype(float).values.T)
    return mean, cov
def get_score(feat, mean, cov):
    result = []
    # 分散共分散行列の逆行列を計算
    cov_i = np.linalg.pinv(cov)
    feat = feat.astype(float).values
    for i in range(len(feat)):
        result.append(distance.mahalanobis(feat[i], mean, cov_i))
    return result
# 正常と異常データを分ける
normal_feat = data.iloc[:s_data,:]
anomaly_feat = data.iloc[s_data+1:,:]
print(normal_feat.shape, anomaly_feat.shape)
# 平均、分散を求める
mean, cov = get_mean_cov(normal_feat)
# マハラノビス距離の計算(正常、異常それぞれ)
normal_score = get_score(normal_feat, mean, cov)
anomaly_score = get_score(anomaly_feat, mean, cov)
# 正常、異常のマハラノビス距離をプロット
plt.plot(normal_score)
plt.plot(anomaly_score)
plt.show()

md_ot = []

num = data.shape[1] # データから因子数を取得

# 因子数ごとに使う直交表を選ぶ (最大因子数は127まで)
if num<12:
    orthogonal_array = orthogonal_array_L12
elif num<31:
    orthogonal_array = orthogonal_array_L32
elif num<63:
    orthogonal_array = orthogonal_array_L64
else:
    print("因子数は最大63まで！")
    sys.exit()
k = 0
for c in orthogonal_array:
    b = []
    for i in range(num):
        if c[i]==2:
            b.append(i) # 直交表で「2」の因子を取り除く為のインデックスの配列を作成
    data_orth = data.drop(columns=data.columns[b]) # 直交表で”2”の項目を取り除く
    normal_feat = data_orth.iloc[:s_data,:]     # 正常データ
    anomaly_feat = data_orth.iloc[s_data+1:,:]  # 異常データ
    # 直交表によって選択された因子のみでマハラノビス計算用の平均、分散を計算
    mean, cov = get_mean_cov(normal_feat)
    
    # 上の平均、分散を用いて、すべての異常値についてのマハラノビス距離を計算
    anomaly_score = get_score(anomaly_feat, mean, cov)
    # 異常値のマハラノビス距離の平均値を結合
    md_ot.append(sum(anomaly_score)/len(anomaly_score))
sum_md = np.zeros((2,num))
count_md_p = np.zeros(num)
count_md_n = np.zeros(num)
sum_msd = np.zeros((2,num))
# マハラノビス距離計算後の処理：因子ありなしの差を可視化する
k=0
for c in orthogonal_array:
    for i in range(num):
        if c[i]==1:
            sum_md[0,i] += md_ot[k]
            sum_msd[0,i] += 1 / (md_ot[k])**2
            count_md_p[i] +=1
            
        else:
            sum_md[1,i] += md_ot[k]
            sum_msd[1,i] += 1 / (md_ot[k])**2
            count_md_n[i] +=1
    k +=1
# それぞれの因子のあり(posi)、なし(nega)のMD平均値を格納する
posi_mean_md = np.zeros(num)
nega_mean_md = np.zeros(num)
for i in range(num):
    posi_mean_md[i] = sum_md[0,i] / count_md_p[i]
    nega_mean_md[i] = sum_md[1,i] / count_md_n[i]
# 各々の因子有無のマハラノビス距離
print(posi_mean_md)
print(nega_mean_md)
# マハラノビス距離についての要因効果図を頑張って書く
# Y軸の範囲を固定
ymin = min([np.min(posi_mean_md), np.min(nega_mean_md)])*0.95
ymax = max([np.max(posi_mean_md), np.max(nega_mean_md)])*1.05
t = np.linspace(0, 1, 2)
y =  np.zeros((num,2))
fig , axes = plt.subplots(1 , num , figsize=(10,5))
for i in range(num):
    y[i] = [posi_mean_md[i], nega_mean_md[i]]
    axes[i].set_title(data.columns.values[i], fontsize=7)
    axes[i].set_xlim(-0.3, 1.3)
    axes[i].set_ylim(ymin, ymax)
    axes[i].plot(t, y[i], color="black", marker='.',markersize=10, label=str(i+1))
    if i>0:
        axes[i].tick_params(labelleft=False, labelright=False,
labeltop=False, length=0, labelsize=7)
    else:
        axes[i].tick_params(length=0, labelsize=7)
plt.show()

これが、今回の肝となるコードになります。プログラム名を「MT_method.py」としておきます。
これを、例によって実行します。

＞ python MT_method.py

このプログラムがやっていることは、以下の3つ。

(1) まず最初の53行の「正常データ」を使って、平均、分散を計算する
(2) 続いて、この「正常データ」とその下53行の「異常データ」それぞれのマハラノビス距離を算出しグラフ化
(3）たくさんの因子を使って計算したのですが、その因子の寄与度を求めるために、正常データを使ってその因子の有り無しの平均、分散を計算し、異常値のマハラノビス距離を計算、その距離値がどれくらい変動するかを「0：あり」「1：なし」で比較する

ということをやってます。

特にこの（3）が肝で、プログラムの途中に見える1と2だけの強烈な配列（直交表）は、この寄与度を算出するために入れております。
このプログラム自体は汎用性があって、最大63因子までのMT法と、その寄与度分析が可能となっております。
先日の東京出張は、まさにこれを作るために行ったようなものです。

では、結果を見ていきましょう。まずは(2)のマハラノビス距離の比較です。

こんなのが出ました。青線が正常、オレンジ色が異常データのものです。
若干ラップしてますが、大体5のところに線を引けば、概ね分けられます。

一般にMT法では、マハラノビス距離が4で分ける、というのが定番ですけど、正常データで4を超えているところが多過ぎて・・・あまりいいデータじゃないってことですかね。

で、このマハラノビス距離の寄与が大きい因子を見るため、（3）を実行した結果が以下。

ここでは「要因効果図」と呼びますが、この図では後ろから4つ目、5つ目が大きく変動してますね。これはそれぞれ、「2500_std」「2750_mean」、つまり「2500～2750Hzの標準偏差」と「2750～3000Hzの平均値」がマハラノビス距離に対して効いていることを示しています。

試しにこの二つを縦横軸にしてプロットしてみると、御覧の通り、この2軸だけでまあまあ正常/異常が分離できてます。
他の因子を縦横にしても、この青とオレンジが結構入り乱れてました。

最初の可視化の時点で、3000Hz付近の音の揺れの違いがあるとみていたので、これはまあ納得の結果、といったところでしょうか。

という具合に、音の異常を身近なもので実験してみました。

もっとも、たまたま分離できただけかもしれません。本来ならもっとたくさんのデータを使ったほうが精度が高いと思われます。

正直言いますと、1日目、12日目だけでなく、2日目、3日目・・・のデータもあるんですが、実はどれも似たようなマハラノビス距離の分離っぷりになってしまいました。1日目から見れば、どいつもこいつも「異常」ということに。
なので、シェーバーの音自体というよりも、録音環境（マイクの位置、髭の剃り方の違い、など）の方が効いているのかもしれません。

また、時系列データの特徴量抽出にはよく「tsfresh」というライブラリが使われます。今回は時系列データの特徴量としては基本の平均、最大、最小、標準偏差だけを使いましたが、tsfreshでは一つの時系列データから779種類の特徴量を作り出してくれます。そちらを使えば、もっといい特徴量が見つかったかも。

まあ、そういうわけで今回の話は、ご参考ということで。

よくわかるMTシステム―品質工学によるパターン認識の新技術

固定リンク | 0 | コメント (0)

2023年7月24日 (月)

iPad miniの激安ケースを購入

久しぶりに、iPad miniの話題です。

といっても、ケースを購入しました、という話だけですが。

こういうのを買いました。お値段は大体千円。

なお、半透明な黒です。いっそ本体色の活かせる透明の方が良かったのですが、それゆえに安かったのだと思います。実際、同じ会社の製品でも、この黒が一番安いです。

おかげさまで、せっかくのパープルが台無しに・・・というほどに酷いわけでもなく、ちょっと浅黒いかなぁ程度の黒です。

ご覧の通り、ピッタリです。

この手のケースのお約束で、こうやってスタンド形態として使用可能

かと思いきや、すぐに滑るので、全然支えになりません。

カバーとして使うだけだと割り切る必要がありそうです。

もっとも、そのつもりで購入しているので、私は問題ありません。

しばらくカバーなしで使ってましたが、これを購入した理由は、ホコリよけのため。

というのも、最近あまりiPad miniを使っておりません。

ここ最近、五十肩になりかけの二の腕と腱鞘炎を抱えた右手首の痛みのコラボレーションが酷くて、iPad miniが持てなくなりました。

３月に買ったキーボードで、文書打ちマシンとして活躍してはいたんですが

iPad mini用のサンワサプライ製スタンド付きのキーボードを買った: EeePCの軌跡

このキーボード、腱鞘炎の右手首と非常に相性が悪くて、結局今は使っておりません。

ということで、電源を落とした状態で置いていることが多いため、ケースを買っておこうと思った次第。

なお、この間まで使っていたケースは、糸がほつれてきたので捨ててしまいました。

早いとこ、この腱鞘炎も直して、復帰させたいものですねぇ。

JEDirect iPad mini 1 2 3 ケース三つ折スタンドオートスリープ機能 (ブラック)

固定リンク | 0 | コメント (0)

2023年7月16日 (日)

4TBのUSB HDDをバックアップ用に購入

先日、メインPCを2TBのSSDに換装しましたが、そのバックアップを取るべくメディアがありませんでした。

てことで、そのバックアップ用に、最近行われていたAmazon Primeデーで安くなったUSB HDDを購入。

BUFFALO製USB HDDです。容量は4TB。2.5インチのUSB供給の電力で動くやつです。

今、2TBのは持ってるんですが、SSD自体が2TBであること、またSSD外にも保存しているデータがあるため、2TBではちょっと足りません。そんな事情で、4TBのHDDを購入しました。

出してみると、思った以上に分厚い。

2TB HDDと比べると一目瞭然です。どうやら2.5インチHDDは、2TBを超えると分厚くなるようです。

この2TB HDDにあるファイル転送も兼ねて、早速つなぎます。

ファイルの転送は、PCに出てきたHDDのアイコンにドラッグ＆ドロップするだけで行けるんですが、SSDをシステム丸ごとバックアップする方法はちょっと工夫が要ります。

この辺の手順は、以下にも書かれてます。

システムイメージを作成しバックアップする方法 | ドスパラサポートFAQ よくあるご質問｜お客様の｢困った｣や｢知りたい｣にお応えします。

スタートメニューの「検索」で「コントロールパネル」と打ち込み、「システムとセキュリティ」‐「バックアップと復元(Windows 7)」というやつを選びます。

その中にある「システムイメージの作成」というのをクリック。

バックアップ元、そしてバックアップ先を聞かれます。バックアップ先にはこの購入したHDDを選択。

あとは終わるのを待つだけ。

・・・だったんですが、ちょっと問題が発生。
夜中に見てみたら、なんとWindowsアップデートによる再起動がかかって、ログイン画面に戻ってました。
このため、バックアップがきちんと終了したかどうかが分からなくなったため、念のためもう一度これを実行する羽目に。

勝手に再起動しないで欲しいなぁ。なんでこうWindowsってユーザーの感覚を度外視するんでしょう？

さすがに一晩かけたら、ちゃんとバックアップが終了しました。

で、もう一つ、Windowsが起動できなくなった時に備えて、回復ドライブというのを作っておきます。
上のドスパラの手順だと「システム修復ディスク」を作ることになってますが、DVD/CDのメディアが必要なため、ちょっと不便です。このため、回復ドライブを使ってUSBメディアから起動できるやつを作っておきます。

これもスタートメニューの検索から「回復ドライブ」と打ち込むと、メニューが出てくるので、USBメモリーを挿して