この記事、何?

前回の続き。
主成分分析をやってみた。

前回の記事

これ。
■Python手遊び(機械学習・外部検証)
https://qiita.com/siinai/items/0a5b89543c9a7bbd243b

主成分分析?

ちょっとだけかじった内容をOutput

多次元空間内で「分散を最大化する」ための行列演算。
その目的は次元圧縮だったりノイズの除去だったり。
具体的には、計算時間の削除。速く答えがもらえる。ラッキー♪
つまり、その分、GridSearch等で試行錯誤できる回数を増やせる。ハッピー♪

方針

機械学習として目的は高い精度のモデルが欲しい。
なので、短い時間で高い精度のモデルが得られれば良しとしようと。
・・・だた、今回、前回の投稿との絡みで「時間」ではなく、「次元数」で比較することにした。
(正直、PLSでは時間比較が難しかったの...orz)

できあがり

# 60-01.py
import psycopg2
import numpy as np
from sklearn.cross_decomposition import PLSRegression
from sklearn.decomposition import PCA
from sklearn.model_selection import GridSearchCV
from datetime import datetime

conn = psycopg2.connect('postgresql://postgres:manager@localhost:5432/postgres')
cur = conn.cursor()

logs = []

# ------------------------------------------------------
# 関数群
# ------------------------------------------------------

# 最後にまとめてログを出力
def output_log():
    for log in logs:
        print(log)

# ログとして出力をStore
def write_log(log):
    print(log)
    logs.append(log)

# 現在時刻の出力
def print_time(comment):
    date = datetime.now()
    tmp = str(date.minute).zfill(2) + ':' + str(date.second).zfill(2) + '.' + str(date.microsecond).zfill(6) + ' : ' + comment
    write_log(tmp)

# 項目のいずれかにnull, nanがあるレコードをにフラグを立てて件数を取得
def update_flg():
    print_time('start: update_flg()')
    sqls = []

    #最初にfalseに設定
    sqls.append('update t_sample set c195 = 0 where c195 != 0')
    #sqls.append('update t_sample set c195 = 0')

    #フラグを更新
    for i in range(1, 115):
        col = 'c' + str(i).zfill(3)
        sql = 'update t_sample set c195 = 1 where c195 = 0 and '
        #sql += ' (' + col + ' in (' + "'nan'" + ', 0) or ' + col + ' is null)'
        sql += ' (' + col + ' in (' + "'nan'" + ') or ' + col + ' is null)'
        sqls.append(sql)

    #実行
    for sql in sqls:
        #print(sql)
        cur.execute(sql)
        conn.commit()

    #フラグのたっていないレコードの件数を把握
    sql = 'select count(sampleid) from t_sample'
    sql += ' where c195 = 0'
    cur.execute(sql)
    result = cur.fetchall()
    return int(result[0][0])

# 項目の値が0であるレコードの件数を把握し、0が一定数以下の項目をリスト化
def get_parameters(record_count):
    print_time('start: get_parameters()')
    condition_rate = 0.9
    parameters = []
    for i in range(115):
        col_id = i+1
        sql = 'select count(sampleid) from t_sample'
        sql += ' where c195 = 0 and c' + str(col_id).zfill(3) + ' = 0'
        cur.execute(sql)
        result = cur.fetchall()
        column_count = int(result[0][0])
        write_log('id=' + str(col_id) + ': ' + str(column_count) + '/' + str(record_count))
        if (column_count / record_count) < condition_rate:
            parameters.append(col_id)
    return parameters

# データを分割
def split_data(record_count):
    sql = 'update t_sample set'
    sql += ' c196 = 0, c197 = 1'
    sql += ' where c195 = 0'
    cur.execute(sql)
    sql = 'update t_sample set'
    sql += ' c196 = 1, c197 = 0'
    sql += ' where c195 = 0'
    sql += '   and sampleid in (select sampleid from t_sample order by random() limit ' + str(record_count*0.8) + ')'
    cur.execute(sql)
    conn.commit()

# 処理時間がかかるので計算済みの値をRetrun
def get_dummy_list():
    print_time('start: get_dummy_list()')
    dummy_list = []
    dummy_list.extend([1, 2, 3, 4, 5, 6, 7, 8, 9])
    dummy_list.extend([11, 12, 13, 14, 15, 16, 17, 18, 19, 20])
    dummy_list.extend([21, 22, 23, 24, 25, 26, 27, 28, 29, 30])
    #C033は極端に大きな値があるので除外...
    #dummy_list.extend([31, 32, 33, 34, 35, 36, 37, 38, 39, 40])
    dummy_list.extend([31, 32, 34, 35, 36, 37, 38, 39, 40])
    dummy_list.extend([41, 42, 43, 44, 45, 46, 47, 48, 49, 50])
    dummy_list.extend([51, 52, 53, 55, 56, 57, 58, 59])
    dummy_list.extend([61, 62, 63, 64, 65, 66, 67, 68, 69, 70])
    dummy_list.extend([72, 74, 75, 76, 78, 79, 80])
    dummy_list.extend([81, 82, 83, 84, 85, 87, 94, 95, 96, 97, 98, 99, 100])
    dummy_list.extend([101, 102, 103, 104, 105, 106, 107, 108, 109])
    dummy_list.extend([111, 112, 113, 114, 115])
    return dummy_list

# 項目一覧からnumpyのarrayを取得
def get_array_from_parameters(parameters, flg):
    print_time('start: get_array_from_parameters()')
    #対象項目からselect句
    sql = 'select sampleid'
    for parameter in parameters:
        sql += ', c' + str(parameter).zfill(3)
    #from句
    sql += ' from t_sample'
    #where句
    sql += ' where c195 = 0'
    #教師データとテストデータを選択
    if flg == 0:
        sql += '   and c196 = 1'
    elif flg == 1:
        sql += '   and c197 = 1'
    #コーディング時の速度向上のために件数を制限
    #sql += ' where c195 = 0 and sampleid between 1 and 3000'
    #order by句
    sql += ' order by sampleid'
    cur.execute(sql)
    result = cur.fetchall()
    #numpyのarrayへ変換
    return np.array(result)

# データを分離して取得
def get_Xy(narray):
    print_time('start: get_Xy()')
    X = np.hstack([narray[:, 1:8], narray[:, 9:]])
    y = narray[:, 8:9]
    return X, y

# PLSモデルを取得
def get_pls_model(X, y, n):
    print_time('start: get_pls_model()')
    pls = PLSRegression(n_components=n)
    pls.fit(X, y)
    return pls

# PCAモデルを取得
def get_pca_model(X, n):
    print_time('start: get_pca_model()')
    pca = PCA(n_components=n)
    pca.fit(X)
    return pca

# ------------------------------------------------------
# 本体
# ------------------------------------------------------

# 項目のいずれかにnull, nanがあるレコードをにフラグを立てて件数を取得
# record_count = update_flg()
# record_count = 22762

# 項目の値が0であるレコードの件数を把握し、0が一定数以下の項目をリスト化
# parameters = get_parameters(record_count)

# データを分割
# split_data(record_count)

# 処理時間がかかるので計算済みの値をRetrun
parameters = get_dummy_list()

# 項目一覧からnumpyのarrayを取得
narray0 = get_array_from_parameters(parameters, 0)
narray1 = get_array_from_parameters(parameters, 1)

cur.close()
conn.close()

# データを分離して取得
X0, y0 = get_Xy(narray0)
X1, y1 = get_Xy(narray1)

# 分離したデータの形を確認
write_log(X0.shape)
write_log(y0.shape)
write_log(X1.shape)
write_log(y1.shape)

# 主成分分析の結果を取得
axes = [2, 5, 10]
pcas = []
for axis in axes:
    pcas.append(get_pca_model(X0, axis))

# PLSモデルを取得
print_time('#PLSモデルを取得')
plss = []
xc0s = []
xc1s = []
for i in range(3):
    xc0s.append(pcas[i].transform(X0))
    xc1s.append(pcas[i].transform(X1))
for i in range(3):
    plss.append(get_pls_model(xc0s[i], y0, axes[i]))

# 各パラメータでのPLSモデルのスコア
write_log('--- X0 score ---')
for i in range(3):
    write_log(plss[i].score(xc0s[i], y0))
write_log('--- X1 score ---')
for i in range(3):
    write_log(plss[i].score(xc1s[i], y1))
write_log(' ')

print_time('おしまい')

output_log()

結果の比較

# 前回の素直にPLSした結果
--- X0 score ---
0.9725608197451778
--- X1 score ---
0.9648904257102765

--- X0 score ---
0.9883775547452807
--- X1 score ---
0.9893330893478123

--- X0 score ---
0.9965083747897437
--- X1 score ---
0.9926695838964517

--- X0 score ---
0.9996288154173829
--- X1 score ---
0.9708325903922626

# 今回、PCAしてPLSした結果
--- X0 score ---
0.8524159223950235
0.9968350863374471
0.9997679790212506
--- X1 score ---
0.851138560747049
0.9961703642694612
0.999509698338353

評価/検討

・・・微妙。
目的を達したような、疑問点が増えたような。

目的を達した?

主成分5軸の時点で外部検証のスコアが0.9961703642694612。
まあ・・・MolWeightの予測だし、極めて高くてもまあ、そんなものといえばそんなもの、かも?

それに対して、PCA抜きでPLSした場合、10項目使っても0.9926695838964517。
都合よく(?)考えると、半分まで次元圧縮して、上回るスコアを出した、といえる、かも?

疑問点

ん・・・疑問点だらけ。今後つぶしてゆこう。

まず、単純なPLSの外部検証のスコアが10因子から20因子に増えたにもかかわらず大幅に落ちている。

0.9926695838964517 → 0.9708325903922626

前回投稿後に気づいたような極端に大きな値とか、そういうものが邪魔しているのかも。
とはいえ、それでこうやってスパイク気味に落ちるのがちょっとイミフ・・・

感想

前回・今回の投稿の目的は外部検証や主成分分析のコードサンプルをアップすること自体でもあるのでいったんアップします。
ただ、いろいろ問題があるので改めて上げ直そう。
っていうか、PCAの効能を見てみたいと思って何でもいいから選んだ予測モデルアルゴリズムがPLSって・・・考えなさすぎでしょ。。。

そういう意味ではプランニング自体が問題だった気がしますわ。
まあ、それに気づいたということも一歩前進ということで。。。