Pythonで機械学習の練習

Pythonでロスアラモス大学が配布しているSHMTools（構造ヘルスモニタリングの機械学習デモ）をやってみた
開発環境は以下の記事通り
renga.hatenadiary.com

www.lanl.gov

とりあえず、データのセットアップと特徴量探索の途中まで

# coding: Shift_JIS
import numpy as np
import pandas as pd
import scipy.io

## データセットの構成
#
# ロスアラモス大学の配布ページからダウンロードする
# SHMTools Software
# http://www.lanl.gov/projects/national-security-education-center/engineering/software/shm-data-sets-and-software.php
#
# \Examples\ExampleData\data3SS.mat
#

data3SS = scipy.io.loadmat("data3SS.mat")

for ch in range(1,5):
    data = np.squeeze(data3SS["dataset"][:,ch,:])

    np.save('data_ch'+str(ch), data)


State = data3SS["states"].T.astype(np.int)
Cond = np.empty((0,1),int)
for i in range(0, State.size):
    if State[i] >=10:
        Cond=np.append(Cond,np.array([['damaged']]),axis=0)
    else:
        Cond=np.append(Cond,np.array([['undamaged']]),axis=0)


df = pd.DataFrame({'State': State.T[0],'Cond':Cond.T[0]})
df.to_csv('StateCond.csv')

# coding: Shift_JIS
import numpy as np
#import scipy.io
from scipy import signal
import matplotlib.pyplot as plt
from statsmodels.tsa import ar_model


## データの読み込み
data = np.load('data_ch4.npy')

dataGood = data[:,1-1]
dataDamaged = data[:,170-1]

Fs = 320
time = np.arange(0, dataGood.size) / Fs


# 時系列データの表示
plt.figure()
plt.subplot(2,1,1);
plt.plot(time,dataGood)
plt.xlabel('time [s]'), plt.ylabel('acc [m/s^2]'), plt.title('undamaged')
plt.axis('tight'), plt.ylim([-2,2])

plt.subplot(2,1,2);
plt.plot(time,dataDamaged);
plt.xlabel('time [s]'), plt.ylabel('acc [m/s^2]'), plt.title('damaged')
plt.axis('tight'), plt.ylim([-2,2])

plt.tight_layout()
plt.show()


## ヒストグラムの表示
plt.figure()
plt.subplot(2, 1, 1)
plt.hist(dataGood, 50)
plt.xlim([-3,3])
plt.title('histgram : undamaged')

plt.subplot(2, 1, 2)
plt.hist(dataDamaged, 50)
plt.xlim([-3,3])
plt.title('histgram : damaged')

plt.tight_layout()
plt.show()


## 周波数スペクトルの表示
freq1, P1 = signal.welch(dataGood, Fs)
freq2, P2 = signal.welch(dataDamaged, Fs)

plt.figure()
plt.subplot(2, 1, 1)
plt.plot(freq1, 10*np.log10(P1), linewidth=2, label="nseg=n/4")
plt.xlabel("Frequency[Hz]"), plt.ylabel("Power/frequency[dB/Hz]")
plt.title('PSD : undamaged')

plt.subplot(2, 1, 2)
plt.plot(freq2, 10*np.log10(P2), linewidth=2, label="nseg=n/4")
plt.xlabel("Frequency[Hz]"), plt.ylabel("Power/frequency[dB/Hz]")
plt.title('PSD : undamaged')

plt.tight_layout()
plt.show()


## 自己回帰モデルの係数取得と表示
arOrder = 10;

# TODO : モジュール化
N = data.shape[1]   # dataの列数取得

coef = np.zeros((N, arOrder+1))

for id in range(0,N):
    model = ar_model.AR(data[:,id])
    results = model.fit(maxlag = arOrder)
    coef[id,:] = results.params

Features = coef

# 可視化
plt.figure()
for id in range(0,N):
    if id>90:
        plt.plot(Features[id],'b')
    else:
        plt.plot(Features[id],'r')
plt.show()