在室判定

Occupancy Detection データセットによる在室判定

データセットは日付と時刻(毎分)、温度(摂氏)、相対湿度(%)、明るさ(Lux)、CO2濃度(ppm)、重量絶対湿度(kgwater-vapor/kg-air)、在室情報(0 or 1)から構成されます。このチュートリアルではニューラルネットワークのモデルを構築し温度、相対湿度、明るさ、CO2濃度、重量絶対湿度から在室か不在かを判定します。
在室判定はスマートビルディングに応用でき、正確な判定は建物の制御システムの節電につながります[1]。

Datasetのリファレンスは以下になります。

Accurate occupancy detection of an office room from light, temperature, relative humidity and CO2 measurements using statistical learning models. Luis M. Candanedo, Véronique Feldheim. Energy and Buildings. Volume 112, 15 January 2016, Pages 28-39.

Required Libaries

  • matplotlib 2.0.2
  • numpy 1.13.1
  • scikit-learn 0.18.2
In [1]:
from __future__ import division, print_function
import datetime
import numpy as np
import pandas as pd

import renom as rm
from renom.optimizer import Adam

from sklearn.metrics import confusion_matrix, classification_report, roc_curve, accuracy_score
from sklearn.preprocessing import StandardScaler, MinMaxScaler

import matplotlib.pyplot as plt

from renom.cuda import set_cuda_active

# If you would like to use GPU, set True, otherwise you should be set to False.
set_cuda_active(False)

Load Data

データセットをhttps://github.com/LuisM78/Occupancy-detection-data からダウンロードし、作業フォルダに入れます。
‘datatraining.txt’をトレーニングデータとして、‘dataset2.txt’をテストデータとして読み込みます。
In [2]:
df_train = pd.read_csv('./datatraining.txt')
df_test = pd.read_csv('./datatest2.txt')

可視化のための前処理

可視化のためにデータセットの前処理をします。
全てのセンサデータを同じグラフ上に描画するために値を0から1に正規化します。
文字列の日付・時刻データをdatetime型に変換します。加えて、状態が「在」になっている範囲を調べます。
In [3]:
temperature = np.array(df_train['Temperature'])
humidity = np.array(df_train['Humidity'])
light = np.array(df_train['Light'])
CO2 = np.array(df_train['CO2'])
humidity_ratio = np.array(df_train['HumidityRatio'])
occupancy = np.array(df_train['Occupancy'])

mms = MinMaxScaler()

# Normalize data.
temperature_norm = mms.fit_transform(temperature.reshape(-1, 1))
humidity_norm = mms.fit_transform(humidity.reshape(-1, 1))
light_norm = mms.fit_transform(light.reshape(-1, 1))
CO2_norm = mms.fit_transform(CO2.reshape(-1, 1))
humidity_ratio_norm =mms.fit_transform(humidity_ratio.reshape(-1, 1))

# Survey the range that status is 'Occupied'.
date = []
occupied_rage = []
isOccupied = False
occupied_start = []
occupied_end = []

for i, d in enumerate(df_train['date']):
    dt = datetime.datetime.strptime(d, "%Y-%m-%d %H:%M:%S")
    date.append(dt)
    if isOccupied == False and occupancy[i] == 1:
        occupied_start.append(dt)
        isOccupied = True
    elif isOccupied == True and occupancy[i] == 0:
        occupied_end.append(date[i-1])
        isOccupied = False

if isOccupied == True:
    occupied_end.append(date[-1])
    isOccupied = False

データの可視化

x軸を日付・時刻、y軸を正規化されたセンサ値としてデータをそれぞれグラフ上に示します。
背景が灰色になっている部分は状態が「在室」になっている範囲です。
このグラフから状態が「在室」であるとき、全てのセンサ値が上昇していることがわかります。
いくつかの範囲では、「不在」であるにもかかわらず湿度や温度が上昇しています。
In [4]:
plt.figure(figsize=(20, 10))
for s, e in zip(occupied_start, occupied_end):
    plt.axvspan(s, e, facecolor='black', alpha=0.2)
plt.plot(date, temperature_norm, label='Temperature')
plt.plot(date, humidity_norm, label='Relative Humidity')
plt.plot(date, light_norm, label='Light')
plt.plot(date, CO2_norm, label='CO2')
plt.plot(date, humidity_ratio_norm, label='Humidity Ratio')
plt.xlim(min(date), max(date))
plt.xlabel('Date')
plt.legend()
plt.show()
../../../_images/notebooks_clustering_occupancy-detection_notebook_9_0.png

学習のための前処理

時刻以外のデータをndarrayに変換し、標準化します。

In [5]:
x_train, y_train = np.array(df_train.iloc[:, 1:6]), np.array(df_train.iloc[:, 6:])
x_test, y_test = np.array(df_test.iloc[:, 1:6]), np.array(df_test.iloc[:, 6:])

sc = StandardScaler()

sc.fit(x_train)
x_train_std = sc.transform(x_train)
x_test_std = sc.transform(x_test)

モデルの定義

出力値を0から1の範囲にするために最後のレイヤーにシグモイド関数を加えます。
閾値を決定し、出力値を0か1に変換します。この変換により、状態が「在」か「不在」かを決定します。
In [6]:
model = rm.Sequential([
    rm.Dense(6),
    rm.Relu(),
    rm.Dense(1),
    rm.Sigmoid()
])

最適化手法の決定

最適化手法としてadamを用います。学習率は0.001とします。

In [7]:
optimizer = Adam(lr=0.001)

学習ループ

numpy.random.permutationを用いてランダムにインデックスを作成し,バッチデータを作成します。このチュートリアルではバッチサイズを128、エポック数を100にし、平均二乗誤差を誤差関数として使用します。

In [8]:
# parameters
EPOCH = 100 # Number of epochs
BATCH =128 # Mini-batch size

# Learning curves
learning_curve = []
test_curve = []

# Training loop
for i in range(1, 1+EPOCH):

    N = x_train_std.shape[0] # Number of records in training data
    perm = np.random.permutation(N)
    train_loss = 0

    for j in range(N//BATCH):
        # Make mini-batch
        index = perm[j*BATCH:(j+1)*BATCH]
        train_batch_x = x_train_std[index]
        train_batch_y = y_train[index]

        # Forward propagation
        with model.train():
            z = model(train_batch_x)
            loss = rm.mean_squared_error(z, train_batch_y)

        # Backpropagation
        grad = loss.grad()

        # Update
        grad.update(optimizer)

        train_loss += loss.as_ndarray()

    # calculate mean squared error for training data
    train_loss = train_loss / (N // BATCH)
    learning_curve.append(train_loss)

    # calculate mean squared error for testidation data
    y_test_pred = model(x_test_std)
    test_loss = rm.mean_squared_error(y_test_pred, y_test).as_ndarray()
    test_curve.append(test_loss)

    # print training progress
    if i % 10 == 0:
        print("Epoch %d - loss: %f - test_loss: %f" % (i, train_loss, test_loss))

print('Finished!')
Epoch 10 - loss: 0.011556 - test_loss: 0.029600
Epoch 20 - loss: 0.007123 - test_loss: 0.020418
Epoch 30 - loss: 0.006330 - test_loss: 0.018179
Epoch 40 - loss: 0.006049 - test_loss: 0.016938
Epoch 50 - loss: 0.005893 - test_loss: 0.015705
Epoch 60 - loss: 0.005739 - test_loss: 0.014743
Epoch 70 - loss: 0.005738 - test_loss: 0.013876
Epoch 80 - loss: 0.005751 - test_loss: 0.013107
Epoch 90 - loss: 0.005648 - test_loss: 0.012373
Epoch 100 - loss: 0.005553 - test_loss: 0.011171
Finished!

学習曲線

In [9]:
plt.figure(figsize=(10, 4))
plt.plot(learning_curve, label='loss')
plt.plot(test_curve, label='test_loss', alpha=0.6)
plt.title('Learning curve')
plt.xlabel("Epoch")
plt.ylabel("MSE")
plt.legend()
plt.grid()
../../../_images/notebooks_clustering_occupancy-detection_notebook_19_0.png

モデルの評価

閾値以上の出力値を1に、未満の出力値を0に変換します。このチュートリアルでは0.5を閾値とします。scilit-learnを使い、正解率(accuracy)、混同行列、適合率(precision)、再現率(recall)、F値を計算します。
混同行列は次のように示されます。
[TN FP] [FN TP]
In [10]:
y_pred = model(x_test_std)
y_pred = np.asarray(y_pred)

threshold = 0.5
y_pred_binary = np.add(y_pred, threshold)
y_pred_binary = np.floor(y_pred_binary)
y_pred_binary = y_pred_binary.astype(np.int64)

print("accuracy : {}".format(accuracy_score(y_test, y_pred_binary)))
print(confusion_matrix(y_test, y_pred_binary))
print(classification_report(y_test, y_pred_binary))
accuracy : 0.9707752255947498
[[7426  277]
 [   8 2041]]
             precision    recall  f1-score   support

          0       1.00      0.96      0.98      7703
          1       0.88      1.00      0.93      2049

avg / total       0.97      0.97      0.97      9752

ROC曲線

ROC(receiver operating characteristics)曲線は分類器の性能を評価する方法のひとつです。
閾値を変化させながらTP、FPをプロットしていくことでその曲線は描かれます。
(0,0)は分類器が必ず0を出力する場合の点で、
(1,1)は分類器が必ず1を出力する場合の点です。
このグラフを見ることで、TPとFPのトレードオフを考慮して閾値を決定することができます。
In [11]:
fpr, tpr, thresholds = roc_curve(y_test, y_pred)

plt.figure(figsize=(10, 4))
plt.plot(fpr, tpr)
plt.xlabel("False Positive Rate")
plt.ylabel("True Positive Rate")
plt.title("ROC Curve")
plt.grid()
../../../_images/notebooks_clustering_occupancy-detection_notebook_23_0.png

References

[1]
Kemal Tutuncu, Ozcan Cataltas, Murat Koklu, OCCUPANCY DETECTION THROUGH LIGHT, TEMPERATURE, HUMIDITY AND CO2 SENSORS USING ANN, Proceedings of ISER 45th International Conference, 2016