Gnuplotでマンデルブロ集合をなんだか書きたくなったのでPythonを利用して書いてみた。まだまだアラが目立つのでブラッシュアップは必要。決定的なのはGnuplotのset pm3dを使っても等高線がグラデーションで表現されないという弱点あり。 mandelbrot_ver1.py #!/usr/bin/env python import numpy as np def columns_save(data, filename='tmp.dat'): ''' PyGnuplotのs関数は等高線を表示させるには力不足なので、別途空白行を出力できるsave関数を作成してデータ出力する。 s関数のように data = [xs,ys] ではなく data= [[x1,y1],[x2,y2] ...]というようなデータ構造に変更。 ''' file = open(filename, 'w') columns = len(data) rows = len(data[0]) for i in range(columns): if data[i] == []: file.write('\n') continue for j in range(rows): file.write(str(data[i][j])) file.write(' ') file.write('\n') if i % 1000 == 0 : file.flush() # write once after every 1000 entries file.close() # write the rest import PyGnuplot as gp gp.colss = columns_save class Mandelbrot: #マンデルブロ集合を生成するクラス。 scale = 2.0 increments = 0.01 def __init__(self): self._xbase = np.arange(-self.scale,self.scale,self.increments) self._ybase = np.arange(-self.scale,self.scale,self.increments) self._cbases = [] self.mandelbrot_set = [] def __gen_cbases(self): #マンデルブロ集合の土台となる複素領域を生成する関数。 for x in self._xbase: self._cbases.append([]) for y in self._ybase: self._cbases[-1].append(complex(x,y)) return def mandelbrot(self,c): #マンデルブロ集合の力学系の計算をする関数。返り値は反復回数。反復回数の数に応じた等高線がマンデルブロ集合のあの絵になる。 z = 0j for n in range(1,2000): z = z**2 + c if np.abs(z) > 2.0: return np.log(n) return 0 def gen_mandelbrot(self): #マンデルブロ集合を生成する関数。 self.__gen_cbases() for cs in self._cbases: for c in cs: self.mandelbrot_set.append([c.real, c.imag, self.mandelbrot(c)]) self.mandelbrot_set.append([]) return if __name__ == "__main__": m = Mandelbrot() m.gen_mandelbrot() gp.colss(m.mandelbrot_set) gp.c(""" #gnuplotの設定読み込み。 set size square; set nosurface; set contour base; set view map; """) gp.c('splot "tmp.dat" using 1:2:3 with lines') マンデルブロ集合をそのままクラスと見立てているので長くなっている。 生成される画像 今後の課題 マンデルブロ集合となっているが解像度が悪い。 ズームすると再起的に同じ図形が出てくると言うところをどう表現すればいいのかわからない。 なんかコードが冗長 等高線が汚いのはcbbaseがメッシュ状になっていないからか？ 上記課題をなんとかしていきたい。

何がしたいか paizaのスキルチェックはWeb上で完結しておりますが、ローカルのIDEやテキストエディタでコードを書きたい場合もあると思います。そこで、独自のinput()を作ってファイルから入力できるようにします。 フォルダ構成 フォルダ構成は以下を想定しております。main.pyが実行するファイルで、utils.pyにinput()の実装を書きます。 paiza/ in_1.txt main.py utils.py 処理内容 実際の処理は以下の通りとなります。とりあえず固定値でin_1.txtを渡すようにしてますが、複数ファイルを扱いたい場合はこのあたりを修正していただければ良いと思います。 utils.py import builtins import linecache line_num = 1 def input() -> str: global line_num target_line = linecache.getline('in_1.txt', line_num) line_num += 1 return target_line.strip() builtins.input = input 使い方 既存のinput()をオーバーライドしているため、基本的にはpaizaで書いたコードを使用していただくことができます。ただしimport utilsだけは忘れないようにしてください。 import utils num = int(input()) (以下略) 次のステップ この処理だけですとアウトプットが目視での検証になってしまうため、そちらも自動化する予定です。

はじめに ラズパイに接続されたカメラ映像のWebストリーミングアプリケーションを作成しました。 また、取得した画像フレームから眠気検知を同時に行っています。 本記事のソースコードは以下に置いています。 https://github.com/AkitoArai709/RaspiSmartCamera-WebStreaming 環境 Raspberry Pi 4 Model B Camera Module V2 For raspberry pi Python : 3.7.3 flask : 1.0.2 opencv : 4.5.1.48 dlib : 19.22.0 imutils : 0.5.4 scipy : 1.6.2 greenlet : 1.0.0 システム構成 　ラズパイに接続されたカメラ映像から画像フレームを取得して、OpenCVとDlibを使用して顔検出と眠気検知を行い、眠気検知の結果を描画した画像フレームをFlaskを使ってローカルネットワーク上でWebストリーミングを行います。 動作デモ 　以下のコマンドでアプリを実行します。 python3 src/main.py 　アプリ実行後にhttp://xx.xx.xx.xx/streamにアクセスするとWebストリーミングが開始されます。目の大きさによって眠気の検知を行い、眠気を検知すると「Look sleepy!」と表示します。 (※「xx.xx.xx.xx」はラズパイのIPアドレス) ファイル構成 . ├─ models │ ├─ dlib ─ shape_predictor_68_face_landmarks.dat │ └─ opencv ─ haarcascade_frontalface_alt2.xml └─ src ├─ static ─ style.css ├─ templates │ ├─ index.html │ └─ stream.html ├─ baseCamera.py ├─ buffer.py ├─ camera.py ├─ cameraEvent.py ├─ detectionSleepiness.py ├─ main.py └─ webStreamingApp.py FlaskによるWebストリーミング 　FlaskとはPython上で動作するウェブアプリケーションのフレームワークで、必要最低限の機能した搭載されていませんが、他フレームワークに比べて動作が軽く簡単に実装が可能です。 　main.pyではwebStreamingApp.pyで実装されたWebアプリを読み込んで実行しています。 　Webアプリでは/、/stream、/video_feedのアドレスを定義し、/でトップページ、/streamでストリーミングのページを表示します。/video_feedではCameraクラスから画像フレームの取得と表示を行い、ストリーミングページに埋め込んでいます。 main.py #!/usr/bin/env python from webStreamingApp import webStreaming def main(): webStreamingApp = webStreaming() webStreamingApp.run(host='0.0.0.0', threaded=True, debug=True) if __name__ == '__main__': main() webStreamingApp.py import os from camera import Camera from flask import Flask, render_template, Response def webStreaming(): return app app = Flask(__name__) def gen(camera): while True: frame = camera.get_frame() yield (b'--frame\r\n' b'Content-Type: image/jpeg\r\n\r\n' + frame + b'\r\n') @app.route('/') def index(): return render_template('index.html') @app.route('/stream') def stream(): return render_template('stream.html') @app.route('/video_feed') def video_feed(): return Response(gen(Camera()), mimetype='multipart/x-mixed-replace; boundary=frame') @app.context_processor def add_staticfile(): def staticfile_cp(fname): path = os.path.join(app.root_path, 'static', fname) mtime = str(int(os.stat(path).st_mtime)) return '/static/' + fname + '?v=' + str(mtime) return dict(staticfile=staticfile_cp) カメラから画像フレームの取得 　camera.py、baseCamera.py、cameraEvent.pyにカメラフレームの処理を実装しています。 Cameraクラス 　OpenCVによりカメラから画像フレームを取得します。OpenCVによるカメラへのアクセスは排他処理になっており、ブラウザごとにカメラアクセスを行うと、複数ブラウザからアクセスした場合にストリーミングが行われているブラウザ以外は固まってしまいます。そこで、画像フレーム取得処理をstaticmethodで実装することで、各インスタンスで同じ処理を参照するようにしています。 camera.py import cv2 from baseCamera import BaseCamera from detectionSleepiness import DetectionSleepiness class Camera(BaseCamera): tick = 0 fpsColor = (0, 255, 0) infApp = DetectionSleepiness() def __init__(self): super().__init__() @staticmethod def frames(): camera = cv2.VideoCapture(0) if not camera.isOpened(): raise RuntimeError('Could not start camera.') while True: # read current frame _, frame = camera.read() frame = Camera.infApp.getDetectResultFrame(frame) yield cv2.imencode('.jpg', Camera.__drawingFps(frame))[1].tobytes() @staticmethod def __drawingFps(frame): fps = 0 if Camera.tick != 0: fps = cv2.getTickFrequency() / (cv2.getTickCount() - Camera.tick) Camera.tick = cv2.getTickCount() return cv2.putText(frame, "FPS:{} ".format(int(fps)), (520, 30), cv2.FONT_HERSHEY_DUPLEX, 1, Camera.fpsColor, 1, cv2.LINE_AA) BaseCameraクラス 　基底クラスにてバックグランドによる画像フレームの取得処理を実行してクラス変数を保持することで、各インスタンスから共通の画像フレームを参照します。また、コンストラクタで起動するスレッドもクラス変数で管理して、処理自体も共通化としています。 baseCamera.py import copy import time import threading from abc import abstractmethod from cameraEvent import CameraEvent class BaseCamera(object): # background thread that reads frames from camera thread = None # current frame is stored here by background thread frame = None # time of last client access to the camera last_access = 0 event = CameraEvent() def __init__(self): if BaseCamera.thread is None: BaseCamera.last_access = time.time() # start background frame thread BaseCamera.thread = threading.Thread(target=self._thread) BaseCamera.thread.start() # wait until frames are available while self.get_frame() is None: time.sleep(0) def get_frame(self): BaseCamera.last_access = time.time() # wait for a signal from the camera thread BaseCamera.event.wait() BaseCamera.event.clear() return BaseCamera.frame @staticmethod @abstractmethod def frames(): raise RuntimeError('Must be implemented by subclasses.') @classmethod def _thread(cls): print('Starting camera thread.') frames_iterator = cls.frames() for frame in frames_iterator: BaseCamera.frame = frame BaseCamera.event.set() # send signal to clients # if there hasn't been any clients asking for frames in # the last 10 seconds then stop the thread if time.time() - BaseCamera.last_access > 10: frames_iterator.close() print('Stopping camera thread due to inactivity.') break BaseCamera.thread = None 　_threadでは、Cameraクラスで定義したframesメソッドを呼び出して画像フレームの取得とクラス変数への保持を行います。framesメソッドはyieldで画像フレームを返しているため、forループで全て処理するまで動作し続けます。クラス変数への保持後にCameraEventのSetメソッドを呼び出して、get_frameメソッド内のwaitの待機を解放します。このイベント処理によって、カメラフレーム取得処理が終わってから、値を返すようにしています。 (※イベント処理についてはCameraEventクラスにて記述しています。) 　保持した画像フレームはget_frameメソッドにて外部モジュールから取得されます。取得時にアクセス時間をlast_accessに保持して一定期間、外部モジュールからの画像フレーム取得がなかった場合にはカメラ処理を停止するようにしています。 CameraEventクラス 　このクラスでは、各ブラウザから呼び出されたフレーム表示スレッドの管理と、イベント処理による同期処理を提供します。 cameraEvent.py import threading import time from greenlet import getcurrent as get_ident class CameraEvent(object): def __init__(self): self.events = {} def wait(self): ident = get_ident() if ident not in self.events: # this is a new client # add an entry for it in the self.events dict # each entry has two elements, a threading.Event() and a timestamp self.events[ident] = [threading.Event(), time.time()] return self.events[ident][0].wait() def set(self): now = time.time() remove = [] for ident, event in self.events.items(): if not event[0].isSet(): # if this client's event is not set, then set it # also update the last set timestamp to now event[0].set() event[1] = now else: # if the client's event is already set, it means the client # did not process a previous frame # if the event stays set for more than 5 seconds, then assume # the client is gone and remove it if now - event[1] > 5: remove.append(ident) for ident in remove: del self.events[ident] def clear(self): self.events[get_ident()][0].clear() 　waitメソッドでは呼び出されたスレッドのID毎にthreading.Event()取得してself.eventsに格納し、Eventのwaitメソッドを呼び出します。このwaitはEventのsetメソッドが呼び出されるまで処理を待機させます。 　setでは待機状態のEventについて、Eventのsetメソッドを呼び出して待機状態の解放を行います。また、非待機状態のイベントが一定時間経過した場合に、管理していたEventの削除を行います。 眠気検知 　DetectionSleepinessクラスにて画像フレームから顔検出と目の大きさによる眠気検知を行います。検知にはCV系処理ライブラリのOpenCVとDlibを使用します。それぞれライブラリで使用する処理モデルについては、modelsフォルダに配置しています。 detectionSleepiness.py import cv2 import dlib from buffer import Buffer from imutils import face_utils from scipy.spatial import distance class DetectionSleepiness: def __init__(self): # Learning result model file path self.faceCascadePath = "./models/opencv/haarcascade_frontalface_alt2.xml" self.faceLandmarksPath = "./models/dlib/shape_predictor_68_face_landmarks.dat" # Learning model self.faceCascade = cv2.CascadeClassifier(self.faceCascadePath) self.faceLandmarksCascade = dlib.shape_predictor(self.faceLandmarksPath) # Drawing color self.faceColor = (255, 255, 255) self.msgColor = (0, 0, 255) # Minimum buffer size required for detection sleepiness self.bufferSize = 50 self.requiredBufferSize = 30 self.SleepinessEARThreshold = 0.58 # EAR buffer # Using for detection sleepiness self.EARbuffer = Buffer(self.bufferSize) def getDetectResultFrame(self, frame): frame, _ = self.__detection(frame, True) return frame def isSleepy(self, frame): _, ret = self.__detection(frame, False) return ret def __detection(self, frame, isDrawing): isSleepy = None # detect person face rect = self.faceCascade.detectMultiScale(frame, scaleFactor=1.11, minNeighbors=3, minSize=(200, 200)) if len(rect) > 0: # resize to face size # convert frame to dlib rectangle resizedFace = self.__resizeFace(frame, rect) faceDlibRectangle = dlib.rectangle(0, 0, resizedFace.shape[1], resizedFace.shape[0]) # caltulation EAR # detect sleepiness left_EAR, right_EAR = self.__getEARs(resizedFace, faceDlibRectangle) isSleepy = self.__detectSleepiness(left_EAR, right_EAR) # drawing result if isDrawing: # drawing a square around the face x, y, w, h = rect[0,:] cv2.rectangle(frame, (x, y), (x+w, y+h), self.faceColor) # drawing left & right EAR(eyes aspect ratio) cv2.putText(frame,"leftEAR:{}".format(round(left_EAR,2)), (10,30), cv2.FONT_HERSHEY_DUPLEX, 1, self.msgColor, 1, 1) cv2.putText(frame,"rightEAR:{}".format(round(right_EAR,2)), (220,30), cv2.FONT_HERSHEY_DUPLEX, 1, self.msgColor, 1, 1) # drawing sleepiness result if isSleepy: cv2.putText(frame,"Look sleepy!", (10,70), cv2.FONT_HERSHEY_DUPLEX, 1, self.msgColor, 1, 1) else: # extract the contents of the buffer if it is not detected self.EARbuffer.pop() return frame, isSleepy def __detectSleepiness(self, left_EAR, right_EAR): ret = True self.EARbuffer.push(left_EAR + right_EAR) if self.EARbuffer.size() >= self.requiredBufferSize and \ self.EARbuffer.getAvg() > self.SleepinessEARThreshold: ret = False return ret def __getEARs(self, frame, face): rect = self.faceLandmarksCascade(frame, face) rect = face_utils.shape_to_np(rect) left_EAR = self.__calcEAR(rect[42:48]) right_EAR = self.__calcEAR(rect[36:42]) return left_EAR, right_EAR def __calcEAR(self, eye): A = distance.euclidean(eye[1], eye[5]) B = distance.euclidean(eye[2], eye[4]) C = distance.euclidean(eye[0], eye[3]) eye_ear = (A + B) / (2.0 * C) return round(eye_ear, 3) def __resizeFace(self, frame, range): # Since the face detection range is small, increase the range. x, y, w, h = range[0,:] w -= 10 y += 10 h += 10 w += 10 face = frame[y :(y + h), x :(x + w)] scale = 480 / h return cv2.resize(face, dsize=None, fx=scale, fy=scale) 　眠気検知には以下の論文を参考に目の高さと幅との間のアスペクト比(ERA)を計算して、直近30個以上のEAR値の平均値から眠気を判定を行います。 引用元：Real-Time Eye Blink Detection using Facial Landmarks 　DetectionSleepinessはAPIとして、getDetectResultFrameメソッドとisSleepyメソッドを提供しています。getDetectResultFrameメソッドでは引数から受け取った画像フレームに対して眠気判定を行った結果を描画した画像フレームを返し、isSleepyメソッドでは眠気判定結果のみを返します。 (※isSleepyメソッドは使用していませんが、別のプロジェクトで使用することを想定して用意しています。) 　眠気検知のためまず、カメラから取得した画像フレームからOpenCVライブラリのself.faceCascade.detectMultiScaleから顔の検出を行います。検出結果から、顔の部分のみを切り抜きてリサイズした画像を使用してDlibライブラリのself.faceLandmarksCascadeから顔のランドマークを検出して37～48の目の部分のみを利用してEARを計算します。リサイズすることで、判定する顔の大きさを統一してカメラからの距離による判定のブレを削減しています。 　顔のランドマークとは、瞳孔、鼻の先など、顔の輪郭のマッピングデータのことを指します。以下の68点の37～48が目の部分になります。 　眠気の検知として以下のパラメータを定義しています。bufferSizeはEAR値を格納するバッファサイズ、requiredBufferSizeは眠気判定の計算に必要なEAR値の最低サイズ、SleepinessEARThresholdは眠気判定の閾値となっています。自身の目の大きさによってSleepinessEARThresholdのパラメータを調整して下さい。 (※現状は結構判定が厳しめなパラメータとなっています。) self.bufferSize = 50 self.requiredBufferSize = 30 self.SleepinessEARThreshold = 0.58 最後に 　ラズパイでリアルタイムな画像処理による眠気検知とWebストリーミングアプリケーションを作成することが出来ました。今回作成した眠気検知モジュールを使用して別のアプリケーションを作成しようと思っています。 (※何を作成するかはまだ未定。。。) 参考サイト 以下の情報を参考にさせていただきました。 FlaskとOpenCVでカメラ画像をストリーミングして複数ブラウザでアクセスする 眠気を判定！目のまばたき検知をDlibとOpenCVを組み合わせて数十行で作る Real-Time Eye Blink Detection using Facial Landmarks

こんにちは。 今回は異常検知の古典的手法であるホテリング理論を Pythonで実践してみようと思います。 まず、ホテリング理論について解説します。 次に、解説した内容を踏まえてPythonによる コーディングを行います。 また、今回の記事を書く上で 以下の書籍の第2章をかなり 参考にさせていただきました。 入門 機械学習による異常検知―Rによる実践ガイド ホテリング理論とは ホテリング理論は一言でいえば、 「観測データが独立に同一の正規分布に従っているという仮定のもとで 外れ値を検出する手法」 です。 以下、数式を用いて説明します。 各対象$i, i=1,\ldots,N$からデータ$x^{(i)}$が観測されているとします。 今回は説明を簡潔にするため、各対象から観測されるデータは一次元データとします。 観測されているデータをまとめて$D=\{x^{(1)},x^{(2)},\ldots, x^{(N)}\}$ と表します。 ここで、$D$の中には異常な値をとるデータは含まれていないか、 含まれていてもその影響は無視できるものと仮定します。 さらに、$D=\{x^{(1)},x^{(2)},\ldots, x^{(N)}\}$の各観測値が 独立に同一の正規分布： N(\mu, \sigma^2)\equiv \frac{1}{\sqrt{2\pi\sigma}} {\rm exp}\{-\frac{1}{2\sigma^2}\left(x-\mu\right)^2\} に従うとします。 このとき、広く知られているように 正規分布のパラメータ$\mu, \sigma^2$の 最尤推定量は以下の式により求められます。 {\hat \mu}=\frac{1}{N}\sum_{i=1}^{N}x^{(i)}, \\ {\hat \sigma^2}=\frac{1}{N}\sum_{i=1}^{N}\left(x^{(i)}-{\hat \mu}\right)^2. 仮定より、${\hat \mu}, {\hat \sigma^2}$を用いることで 正常時の確率分布を求めることができます。 次に異常検知特有の概念である、異常度を考えます。 異常度とは「異常の度合いを関数を用いて数値化したもの」です。 ここで、異常度を計算するための関数が満たしていて ほしい性質として、以下の2つが挙げられます。 正常時に出現しやすい観測値は異常度が低い 正常時に出現しにくい観測値は異常度が高い これらの性質を満たすことから、観測値$x$の一つの自然な 異常度$a(x)$として a(x)= -{\rm ln}p(x) がよく採用されるそうです。 ここで、$p(x)$は正常時の確率分布、${\rm ln}$は自然対数です。 今回も自然対数によって異常度を定義すれば、 \begin{align} &-{\rm ln}\left( \frac{1}{\sqrt{2\pi{\hat\sigma}}} {\rm exp}\{-\frac{1}{2{\hat \sigma}^2}\left(x-{\hat \mu}\right)^2\} \right) \\ &=\frac{1}{2}{\rm ln}\left(2\pi{\hat \sigma}^2\right) + \frac{1}{2{\hat \sigma}^2}\left( x-{\hat \mu} \right)^2 \end{align} と表せます。 実際には上式の右辺の第1項は$x$に依存しないので 計算する必要がなく、第2項の形を整えるために2をかけた a(x)\equiv \frac{\left(x-{\hat \mu}\right)^2}{{\hat \sigma}^2} を異常度として用いることができます。 $a(x)$を解釈してみましょう。 分子の$\left(x-{\hat \mu}\right)^2$は、$x$とデータの中心${\hat \mu}$のズレの大きさ に対応しています。 一方、分母はデータの標本分散なので、$a(x)$は 「$x$がデータの中心からどれだけ離れているかを データのばらつきで割ったもの」と言えます。 ${\hat \sigma}^2$で割ることで、 ばらつきが大きいデータに対しては、 $x$がデータの中心から離れていても異常度が 大きくなり過ぎないように調整しているわけですね。 さて、ここまでの工程をまとめると以下のように表せます。 STEP1：観測されているデータ$D$から正常時の確率分布を推定 STEP2：推定した確率分布を用いて異常度を計算するための関数を導出 よって、STEP2までで、 新たに観測されたデータ$x'$に対して、 そのデータの異常度$a(x')$を計算できるようになりました。 ここまでくればもう一息です。 最後のSTEPは「閾値（しきいち）の設定」です。 閾値とは、「異常と正常の境目にあたる値」です。 これを客観的に決めることで、新たに観測されたデータが 異常かどうかを判定できるようになります。 閾値を決めるために次の結果を用います。 $a(x')$の定数倍は、自由度$(1,N-1)$のF分布に従う： \frac{N-1}{N+1}a(x') \sim F(1,N-1) 特に$N \gg 1$のとき、$a(x')$は自由度$1$、スケール因子$1$のカイ二乗分布に従う： a(x') \sim \chi^2(1,1). この結果は、これまで置いていた仮定に加え、 新たに観測された$x'$は$D=\{x^{(1)},x^{(2)},\ldots, x^{(N)}\}$と 同じ分布に従う、という仮定の下で成り立ちます。 基本的に$N \gg 1$は成り立つため、 $a(x')$はカイ二乗分布に従っていると考えることができます。 よって、 確率値$\alpha$（ex. $0.01$や$0.05$）を用いて、 カイ二乗分布の表から異常度の閾値$a_{\rm th}$を決めること ができます。 ここまでがホテリング理論の一連の流れです。 Pythonを用いたホテリング理論の実装 Pythonで実装してみましょう。 ここでは、正常データとして $\mu=170, \sigma^2=5$の正規分布に従うような 人工データを100個生成します。 成人男性の身長のデータをイメージしています。 import numpy as np np.random.seed(10) D = np.random.normal( loc = 170, scale = 5, size = 100, ) 次に、人工データを用いてパラメータを推定します。 mean_D = sum(a)/len(a) sigma_D =1/(len(a)) * sum((a - mean_a) ** 2) print(mean_D) # 標本平均 print(sigma_D) # 標本分散 # 170.39708331468447 # 23.379133757408734 推定した結果、標本平均は170.40、標本分散は23.38でした。 異常度を計算し、カイ二乗分布による５%水準の閾値を用いて 新たに得られたデータが異常かどうかを判断する関数を 作成します。 def judge (x): ''' 新たに得られたデータが異常な値かどうか判断する。 Args: ｘ: 新たに得られたデータ Returns: str: True なら「異常な値です」、False　なら「正常な値です」 ''' from scipy.stats import chi2 ax = ((int(x) - mean_D) / np.sqrt(sigma_D)) ** 2 th = chi2.ppf(q = 0.95, df = 1) if ax > th: return "異常な値です" else: return "正常な値です" これで、新たに得られたデータが異常かどうか判定できるようになりました。 試しに$x'_1=180, x'_2=162$が異常かどうか判定してみました。 print(judge (180)) print(judge(162)) # 異常な値です # 正常な値です 180cmの人は異常、162cmの人は正常と判定されました。 まとめ 今回は一次元データの場合のホテリング理論を説明しました。 次回は多次元データの場合のホテリング理論を説明したいと思います。 最後まで読んでいただき、ありがとうございました。

こんにちは。 今回は異常検知の古典的手法であるホテリング理論を Pythonで実践してみようと思います。 まず、ホテリング理論について解説します。 次に、解説した内容を踏まえてPythonによる コーディングを行います。 また、今回の記事を書く上で 以下の書籍の第2章をかなり 参考にさせていただきました。 入門 機械学習による異常検知―Rによる実践ガイド ホテリング理論とは ホテリング理論は一言でいえば、 「観測データが独立に同一の正規分布に従っているという仮定のもとで 外れ値を検出する手法」 です。 以下、数式を用いて説明します。 各対象$i, i=1,\ldots,N$からデータ$x^{(i)}$が観測されているとします。 今回は説明を簡潔にするため、各対象から観測されるデータは1次元データとします。 観測されているデータをまとめて$D=\{x^{(1)},x^{(2)},\ldots, x^{(N)}\}$ と表します。 ここで、$D$の中には異常な値をとるデータは含まれていないか、 含まれていてもその影響は無視できるものと仮定します。 さらに、$D=\{x^{(1)},x^{(2)},\ldots, x^{(N)}\}$の各観測値が 独立に同一の正規分布： N(\mu, \sigma^2)\equiv \frac{1}{\sqrt{2\pi\sigma}} {\rm exp}\{-\frac{1}{2\sigma^2}\left(x-\mu\right)^2\} に従うとします。 このとき、広く知られているように 正規分布のパラメータ$\mu, \sigma^2$の 最尤推定量は以下の式により求められます。 {\hat \mu}=\frac{1}{N}\sum_{i=1}^{N}x^{(i)}, \\ {\hat \sigma^2}=\frac{1}{N}\sum_{i=1}^{N}\left(x^{(i)}-{\hat \mu}\right)^2. 仮定より、${\hat \mu}, {\hat \sigma^2}$を用いることで 正常時の確率分布を求めることができます。 次に異常検知特有の概念である、異常度を考えます。 異常度とは「異常の度合いを関数を用いて数値化したもの」です。 ここで、異常度を計算するための関数が満たしていて ほしい性質として、以下の2つが挙げられます。 正常時に出現しやすい観測値は異常度が低い 正常時に出現しにくい観測値は異常度が高い これらの性質を満たすことから、観測値$x$の一つの自然な 異常度$a(x)$として a(x)= -{\rm ln}p(x) がよく採用されるそうです。 ここで、$p(x)$は正常時の確率分布、${\rm ln}$は自然対数です。 今回も自然対数によって異常度を定義すれば、 \begin{align} &-{\rm ln}\left( \frac{1}{\sqrt{2\pi{\hat\sigma}}} {\rm exp}\{-\frac{1}{2{\hat \sigma}^2}\left(x-{\hat \mu}\right)^2\} \right) \\ &=\frac{1}{2}{\rm ln}\left(2\pi{\hat \sigma}^2\right) + \frac{1}{2{\hat \sigma}^2}\left( x-{\hat \mu} \right)^2 \end{align} と表せます。 実際には上式の右辺の第1項は$x$に依存しないので 計算する必要がなく、第2項の形を整えるために2をかけた a(x)\equiv \frac{\left(x-{\hat \mu}\right)^2}{{\hat \sigma}^2} を異常度として用いることができます。 $a(x)$を解釈してみましょう。 分子の$\left(x-{\hat \mu}\right)^2$は、$x$とデータの中心${\hat \mu}$のズレの大きさ に対応しています。 一方、分母はデータの標本分散なので、$a(x)$は 「$x$がデータの中心からどれだけ離れているかを データのばらつきで割ったもの」と言えます。 ${\hat \sigma}^2$で割ることで、 ばらつきが大きいデータに対しては、 $x$がデータの中心から離れていても異常度が 大きくなり過ぎないように調整しているわけですね。 さて、ここまでの工程をまとめると以下のように表せます。 STEP1：観測されているデータ$D$から正常時の確率分布を推定 STEP2：推定した確率分布を用いて異常度を計算するための関数を導出 よって、STEP2までで、 新たに観測されたデータ$x'$に対して、 そのデータの異常度$a(x')$を計算できるようになりました。 ここまでくればもう一息です。 最後のSTEPは「閾値（しきいち）の設定」です。 閾値とは、「異常と正常の境目にあたる値」です。 これを客観的に決めることで、新たに観測されたデータが 異常かどうかを判定できるようになります。 閾値を決めるために次の結果を用います。 $a(x')$の定数倍は、自由度$(1,N-1)$のF分布に従う： \frac{N-1}{N+1}a(x') \sim F(1,N-1) 特に$N \gg 1$のとき、$a(x')$は自由度$1$、スケール因子$1$のカイ二乗分布に従う： a(x') \sim \chi^2(1,1). この結果は、これまで置いていた仮定に加え、 新たに観測された$x'$は$D=\{x^{(1)},x^{(2)},\ldots, x^{(N)}\}$と 同じ分布に従う、という仮定の下で成り立ちます。 基本的に$N \gg 1$は成り立つため、 $a(x')$はカイ二乗分布に従っていると考えることができます。 よって、 確率値$\alpha$（ex. $0.01$や$0.05$）を用いて、 カイ二乗分布の表から異常度の閾値$a_{\rm th}$を決めること ができます。 ここまでがホテリング理論の一連の流れです。 Pythonを用いたホテリング理論の実装 Pythonで実装してみましょう。 ここでは、正常データとして $\mu=170, \sigma^2=5$の正規分布に従うような 人工データを100個生成します。 成人男性の身長のデータをイメージしています。 import numpy as np np.random.seed(10) D = np.random.normal( loc = 170, scale = 5, size = 100, ) 次に、人工データを用いてパラメータを推定します。 mean_D = sum(D) / len(D) sigma_D = 1 / (len(D)) * sum((D - mean_D) ** 2) print(mean_D) # 標本平均 print(sigma_D) # 標本分散 # 170.39708331468447 # 23.379133757408734 推定した結果、標本平均は170.40、標本分散は23.38でした。 異常度を計算し、カイ二乗分布による５%水準の閾値を用いて 新たに得られたデータが異常かどうかを判断する関数を 作成します。 def judge(x): ''' 新たに得られたデータが異常な値かどうか判断する。 Args: ｘ: 新たに得られたデータ Returns: str: Trueなら「異常な値です」、Falseなら「正常な値です」 ''' from scipy.stats import chi2 ax = ((int(x) - mean_D) / np.sqrt(sigma_D)) ** 2 th = chi2.ppf(q = 0.95, df = 1) if ax > th: return "異常な値です" else: return "正常な値です" これで、新たに得られたデータが異常かどうか判定できるようになりました。 試しに$x'_1=180, x'_2=162$が異常かどうか判定してみました。 print(judge(180)) print(judge(162)) # 異常な値です # 正常な値です 180cmの人は異常、162cmの人は正常と判定されました。 まとめ 今回は1次元データの場合のホテリング理論を説明しました。 次回は多次元データの場合のホテリング理論を説明したいと思います。 最後まで読んでいただき、ありがとうございました。

Pythonウェブページのデータを取得 Pythonにライブラリrequests、解析ライブラリBeautifulSoupを介して情報を取得します、ここではYahoo（https://www.yahoo.com/）を例にします、ホームページTrending Nowの内容を取得します。 1）まずはライブラリrequests、BeautifulSoupをインポート 2）requests.get(url)を使ってウェブページのすべてのコンテンツを取得 java requests.get(url) 3）取得したウェブページのコンテンツをBeautifulSoupタイプに変換 java bs = BeautifulSoup(result, 'html.parser') 4）解析Yahoo（https://www.yahoo.com/）のTrending Nowエリア、探し出す他の内容との違い、分析し発見するclass属性は唯一です。この属性を使って値を取得します。 find_all(class_='XXX')方法を使って値を取得します。 ※classがpythonのキーワード、だからclass_を使用します。 その他ラベル、属性は直接使用できます。 例： python find_all('a') find_all(id='link2') find_all(["a", "span"]) など。 イメージ： ソースコード # ライブラリ import requests # 解析ライブラリ from bs4 import BeautifulSoup # ウェブページurl r = requests.get("https://www.yahoo.com/") # タイプ result = r.text # BeautifulSoupタイプに変換 bs = BeautifulSoup(result, 'html.parser') print("データ：") # class属性を使って値を取得します data1 = bs.find_all(class_='Mstart(-2px) C($c-fuji-grape-jelly):h C($c-fuji-inkwell) Fz(12px) Fw(600)') # ループ出力 for i in data1: print(i.text) 実行結果 C:\>"C:\Program Files\Python37\python.exe" D:/WeChat_Robot/CatchData/catchdata.py データ： Carey Mulligan Hideki Matsuyama Kid Cudi Prince Phillip Biglietti da visita Antivirus Kaspersky Forni elettrici Camicie uomo Mascherine online Spesa online Process finished with exit code 0

最近ようやくCallBackに手を付ける気になりましたので 手始めにTensorflowのHistoryを改造する方法について調べました。 この記事では次の環境を使用しております。 Python 3.7.7 tensorflow 2.4.1 標準Historyの問題点 Tensorflowではfitの戻り値としてHistoryオブジェクトが返ってきます 出力される値は、Batch毎の出力を平均した値です LossやAccならまだ使えると思いますが F1Scoreでは全く用をなしません 以下の記事を参照しました。 https://qiita.com/koshian2/items/81abfc0a75ea99f726b9 HistoryCallBack 標準が使えないならば自作してみましょう。 まずはDocument巡りです。 TensorflowのDocument https://www.tensorflow.org/api_docs/python/tf/keras/callbacks/History 生成されたHistoryの使い方については書いてありますが CallBackとしてHistoryを加工する情報はありませんでした。 KerasのDocument https://keras.io/ja/callbacks/ 　検索キー 　　例:損失を記録する 有益な情報が記載されていました。 やってみましょう。 コード import numpy as np,tensorflow as tf class myHistory(tf.keras.callbacks.History): def __init__(self,**kwargs): super().__init__(**kwargs) def on_epoch_end(self, epoch, logs={}): print("on_epoch_end",logs) # super().on_epoch_end(epoch,logs) super().on_epoch_end(epoch, {"このDictが":"Historyに出力される"}) import numpy as np,tensorflow as tf x_train = np.ones(shape=(1,1)) y_train = np.ones(shape=(1,1)) inputs = tf.keras.layers.Input(shape=(1,)) pred = tf.keras.layers.Dense(1,activation="softmax")(inputs) model = tf.keras.models.Model(inputs=inputs, outputs=pred) model.compile(loss='mse',optimizer='sgd') modelhist = model.fit(x_train,y_train, callbacks=[myHistory()]) print("history out",modelhist.history) 出力結果 1/1 [==============================] - 0s 125ms/step - loss: 0.0000e+00 on_epoch_end {'loss': 0.0} history out {'このDictが': ['Historyに出力される']} 結論 tf.keras.callbacks.Historyを継承して on_epoch_end の logs を加工し super().on_epoch_end に渡せば良い

Djangoでプロジェクトをつくってデプロイで失敗し、何時間も時間を無駄にしてしまった方は多いと思います。私もその一人で、毎回、試行錯誤を繰り返していますが、うまくいったケースを参考までにメモとして残しておきます。なお、これでうまくいかないケースやセキュリティ上大丈夫かなどの問題も想定されるので、記事の内容の適用については、あくまで自己責任でお願いします。 1. settings.pyを本番環境用に修正する settings.py（重要な箇所のみコードを抽出しています）はデフォルトというかローカル開発環境では、以下のようになっていると思います。ローカル開発環境で動作させるだけであれば、アプリを追加した場合にアプリ名をINSTALLED_APPSに追加する程度の修正のみだと思います。（今後ローカル環境と本番環境で使い分けるため、ローカル環境用のものは、settings_local.pyのようにファイル名を変えておくことをおすすめします。） settings_local.py from pathlib import Path BASE_DIR = Path(__file__).resolve().parent.parent DEBUG = True ALLOWED_HOSTS = [] INSTALLED_APPS = [ 'django.contrib.admin', 'django.contrib.auth', 'django.contrib.contenttypes', 'django.contrib.sessions', 'django.contrib.messages', 'django.contrib.staticfiles', ] MIDDLEWARE = [ 'django.middleware.security.SecurityMiddleware', 'django.contrib.sessions.middleware.SessionMiddleware', 'django.middleware.common.CommonMiddleware', 'django.middleware.csrf.CsrfViewMiddleware', 'django.contrib.auth.middleware.AuthenticationMiddleware', 'django.contrib.messages.middleware.MessageMiddleware', 'django.middleware.clickjacking.XFrameOptionsMiddleware', #　ローカル環境でも、アプリを追加（manage.py startapp ＜アプリ名＞）した場合、 # ここにアプリ名を登録していると思いますが、それはデプロイ後もそのまま維持します。 'test_app', ] ROOT_URLCONF = 'apps.urls' TEMPLATES = [ { 'BACKEND': 'django.template.backends.django.DjangoTemplates', 'DIRS': [], 'APP_DIRS': True, 'OPTIONS': { 'context_processors': [ 'django.template.context_processors.debug', 'django.template.context_processors.request', 'django.contrib.auth.context_processors.auth', 'django.contrib.messages.context_processors.messages', ], }, }, ] DATABASES = { 'default': { 'ENGINE': 'django.db.backends.sqlite3', 'NAME': BASE_DIR / 'db.sqlite3', } } STATIC_URL = '/static/' Herokuにデプロイして動作させるためには、上記のローカル開発環境用の設定ファイルに何箇所か修正を入れる必要があります。以下が、実際に私がherokuへのデプロイに成功したファイルです。 おそらくポイントは大きく３つ。１つ目は、ホワイトノイズを使えるようにすること、２つ目は、テンプレートの場所の指定、３つ目は、静的ファイルの場所の視点、です。 このあたりをミスすることでエラーが発生しているように思いますので、ご自身の環境にあわせてチェックをされることをおすすめします。 settings.py from pathlib import Path import os import whitenoise #修正 BASE_DIR = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) #追加 SETTINGS_PATH = os.path.dirname(os.path.dirname(__file__)) #ここは、最初はTrueでやってみて、エラーが出ないか確認してからFalseに変えるのも一案。 #（Falseにするとエラーがあったときに原因の究明ができない） DEBUG = False #デプロイしたホスト名を記載する ALLOWED_HOSTS = ['abc.fff.com'] INSTALLED_APPS = [ 'django.contrib.admin', 'django.contrib.auth', 'django.contrib.contenttypes', 'django.contrib.sessions', 'django.contrib.messages', 'django.contrib.staticfiles', #アプリ名はそのまま。もし動かない場合には、'test_app.apps.Test_appConfig'と #する必要があるかもしれません。 'test_app', ] MIDDLEWARE = [ 'django.middleware.security.SecurityMiddleware', 'django.contrib.sessions.middleware.SessionMiddleware', 'django.middleware.common.CommonMiddleware', 'django.middleware.csrf.CsrfViewMiddleware', 'django.contrib.auth.middleware.AuthenticationMiddleware', 'django.contrib.messages.middleware.MessageMiddleware', 'django.middleware.clickjacking.XFrameOptionsMiddleware', #WhiteNoiseを追加。SecurytyMiddlewareがその前に入っているのことを確認しましょう。 #settings.pyの頭にimport whitenoiseとしてモジュールをインポートしておく。 'whitenoise.middleware.WhiteNoiseMiddleware', ] TEMPLATES = [ { 'BACKEND': 'django.template.backends.django.DjangoTemplates', #ここがポイント。テンプレートディレクトリをこのような形で明示してやることが必要。 'DIRS': [os.path.join(BASE_DIR, 'templates')], 'APP_DIRS': True, 'OPTIONS': { 'context_processors': [ 'django.template.context_processors.debug', 'django.template.context_processors.request', 'django.contrib.auth.context_processors.auth', 'django.contrib.messages.context_processors.messages', ], }, }, ] # データベースは特にいじらなくても、問題なく動きました。HerokuはPostgreSQLに対応している # ので、そちらに移植するという選択肢もあるかもしれない。 DATABASES = { 'default': { 'ENGINE': 'django.db.backends.sqlite3', 'NAME': os.path.join(BASE_DIR, 'db.sqlite3'), } } #ここからがポイント。静的ファイルの扱いは常にトラブルになりがちなので慎重に設定しましょう。 #ここは特に変更する必要なし。 STATIC_URL = '/static/' #これを新たに追加する必要があります。静的ファイルのルートを指定します。 STATIC_ROOT = os.path.join(BASE_DIR, 'staticfiles') # STATICFILES_DIRSの設定をしているケースもあるようですが、試行錯誤した結果、個人的には # 不要ではないかと思います。以下はなくても動きました。 # STATICFILES_DIRS = ( # os.path.join(BASE_DIR, 'static'), #) #これも追加。Whitenoiseを使うために必要らしい。 STATICFILES_STORAGE = 'whitenoise.storage.CompressedManifestStaticFilesStorage' 2. Heroku上の環境構築に必要な設定ファイルの準備 ３つのファイルを準備します。runtime.txtはPythonのインストール用（これは不要との説もあります。）、Procfileは、Heroku上でmanage.py runserverを実行して動かすために必要なファイル（説明が間違っているかも）、requirements.txtは、モジュールのインストールに必要なファイルです。参考までに、私が使用しているそれぞれのファイルを掲載しておきます。runtime.txtに記載するPythonのバージョンは、Herokuで対応しているものが限定的なので、Herokuのサイトでバージョンを確認された方が良いと思います。 runtime.txt python-3.8.10 Procfile web: gunicorn quiz.wsgi --log-file - requirements.txt dj-database-url==0.5.0 Django==2.2.6 django-heroku==0.3.1 gunicorn psycopg2==2.8.4 pytz==2019.3 sqlparse==0.3.0 whitenoise==4.1.4 3. Herokuにデプロイ まず、PythonのプロジェクトをHerokuにデプロイする際に、はじめての方は公式サイトのチュートリアルにざっと目を通してなんとなく感覚をつかむことをオススメします。 https://devcenter.heroku.com/ja/articles/getting-started-with-python Herokuのアカウント作成方法は各種サイトに掲載されていると思いますので、ここでは省略したいと思います。また、GitHubとリンクする方法もありますが、とりあえず作ったものをアップするだけの方法を紹介したいと思います。バージョン管理が必要な場合には、GitHubを使ってデプロイをしてみてください。 ＜（参考）Heroku CLIのインストール（MACの場合）＞ $ brew install heroku/brew/heroku ＜（参考）アカウント作成後のログイン＞ $ heroku login ＜（参考）プロジェクト作成　※任意のプロジェクト名が自動生成されます。 Herokuのウェブサイトから作成しても良いかも。＞ $ heroku create　 ＜アプリ名を表示＞ $ heroku apps ＜Gitと特定のアプリを連携＞ $ heroku git:remote -a アプリ名 $ git init $ git add . $ git commit -m "initial commit" $ git push heroku master ＜アプリを動かすコマンド＞ $ heroku ps:scale web=1 ＜アプリを開くコマンド＞ $ heroku open ＜アプリを止めるコマンド＞ $ heroku ps:scale web=0 ＜（参考）マイグレーションと特権管理者の設定＞ $ heroku run python manage.py migrate $ heroku run python manage.py createsuperuser （参考）Collectstaticについて 自動的にHerokuではCollectstatic（静的ファイルを自動的に収集）作業をやってくれるのですが、どうもエラーが出るという記事をよく見るように思います。私はエラーが出なかったのですが、ローカル開発環境下でcollectstaticはやっておいて、Heroku上では行わないという手法もあるようです。その場合の無効化コマンドは以下の通り。 $ heroku config:set DISABLE_COLLECTSTATIC=1

始めに Pythonを利用して、csvファイルをそのままpostgresqlにインポートするスクリプトを作成した。sqlaclchemyとpandasを利用することで、テーブルの作成の手間がなくなったり、DBの違いも吸収できるようになった。 postgresqlの準備 実験用のpsotgresqlは、docker-composeを利用して立ち上げる。 version: '3' services: postgres: image: postgres:10 container_name: postgresql ports: - "5432:5432" environment: POSTGRES_USER: postgresql_user POSTGRES_PASSWORD: postgresql_password volumes: - ./db:/var/lib/postgresql/data 起動コマンドを以下に示す。 docker-compose up -d Docker Compose is now in the Docker CLI, try `docker compose up` Creating csv2db_postgres_1 ... done csvをpostgresqlにインポート 必要ライブラリのインストール 利用するライブラリはpandasとsqlalchemyの2つのみ。各ライブラリのバージョンは以下の通り。 SQLAlchemy==1.4.15 pandas==1.2.4 以下のコマンドでインストールを実施する。 pip install pandas sqlalchemy スクリプトの作成 csvファイルをDBにインポートするPythonスクリプトは以下の通り。 main.py import argparse import pandas as pd from sqlalchemy import create_engine if __name__ == "__main__": parser = argparse.ArgumentParser() parser.add_argument("table_name", help="テーブル名") parser.add_argument("csv_path", help="読み込むCSVファイルのパス") args = parser.parse_args() # 引数の設定 connection_config = { 'user': 'postgresql_user', 'password': 'postgresql_password', 'host': 'localhost', 'port': '5432', 'database': 'postgres', 'database_type': 'postgresql' # 接続するDBのタイプ(mysqlの時: mysql) } engine = create_engine('{database_type}://{user}:{password}@{host}:{port}/{database}'.format(**connection_config), encoding='utf-8') df = pd.read_csv(args.csv_path) # PostgreSQLに書き込む df.to_sql(args.table_name, con=engine, if_exists='replace', index=False ) table_nameにDBのテーブル名を指定し、csv_pathに読み込むCSVファイルのパスを指定する。今回はkaggleのタイタニックの学習データをテーブルに保存してみる。 以下のURLからcsvをダウンロードして、dataフォルダ配下に配置する。 https://www.kaggle.com/c/titanic/data 実行コマンドを以下に示す。 python main.py taitanic_train data/train.csv テーブルにログインして保存されていることを確認。 docker exec -ti postgresql /bin/bash # Dockerのコンテナにログイン root@3ba4a2adf245:/# psql -d postgres -U postgresql_user psql (10.16 (Debian 10.16-1.pgdg90+1)) Type "help" for help. postgres=# postgres=# \z Access privileges Schema | Name | Type | Access privileges | Column privileges | Policies --------+----------------+-------+-------------------+-------------------+---------- public | taitanic_train | table | | | (1 row) postgres=# select * from taitanic_train; # Select文でテーブルの中身を表示。 PassengerId | Survived | Pclass | Name | Sex | Age | SibSp | Parch | Ticket | Fare | Cabin | Embarked -------------+----------+--------+------------------------------------------------------------------------------------+--------+------+-------+-------+--------------------+----------+ -----------------+---------- 1 | 0 | 3 | Braund, Mr. Owen Harris | male | 22 | 1 | 0 | A/5 21171 | 7.25 | | S 2 | 1 | 1 | Cumings, Mrs. John Bradley (Florence Briggs Thayer) | female | 38 | 1 | 0 | PC 17599 | 71.2833 | C85 | C 3 | 1 | 3 | Heikkinen, Miss. Laina | female | 26 | 0 | 0 | STON/O2. 3101282 | 7.925 | | S 4 | 1 | 1 | Futrelle, Mrs. Jacques Heath (Lily May Peel) | female | 35 | 1 | 0 | 113803 | 53.1 | C123 | S 5 | 0 | 3 | Allen, Mr. William Henry | male | 35 | 0 | 0 | 373450 | 8.05 | | S ===== 以下省略 ==== 指定したテーブル名でデータをインポートできた！ 保存されているデータの型を確認してみると、文字列はtext、数字はbigintやdouble precisionが適応されており、最低限のマッピングはなされているよう。もちろんNullやCollactionは未定義の状態になるので、これらの定義を付与したい場合は手動で定義する必要がある。 postgres=# \d taitanic_train Table "public.taitanic_train" Column | Type | Collation | Nullable | Default -------------+------------------+-----------+----------+--------- PassengerId | bigint | | | Survived | bigint | | | Pclass | bigint | | | Name | text | | | Sex | text | | | Age | double precision | | | SibSp | bigint | | | Parch | bigint | | | Ticket | text | | | Fare | double precision | | | Cabin | text | | | Embarked | text | | | 終わりに 今まで地味にめんどくさかったDBへcsvファイルのインポートが簡単に出来るようになった。細かい設定ができないのでシステム開発でDBにデータを挿入する際は利用できないが、データ分析のために1度DBに入れたい、といった場合には使えそう。 GitHub - sqlalchemy/sqlalchemy: The Database Toolkit for Python

What is this? (これは何？) kaggle の House Priceに取り組みます。 狙い kaggleの初心者向け記事はたくさんありますが、あんまり初心者向けじゃない気がしています。（私のレベルが低すぎて...） なので、より低レベルな内容に絞った記事（メモ）を作ろうとしました。 What is the skill of writer?（書いた人の技術力はどんなもの？） 記事を書いた人の能力は以下の通り 状況が近い人は、何か参考になるものがあるかもしれません 学校：文系大卒 お仕事：WEBアプリケーションエンジニア歴満4年 その他のエンジニア経歴はなし 言語：そもそも、Pythonを仕事ではあんまり使ってない 数学：Ⅲ/Ｃを習ってない 年代が違う人のために説明すると、微分方程式や行列は習っていない、ということです（年齢がばれる...） 理系の大学を出た人には何の参考にもならないのではないかと思っています... First step! (はじめの一歩) 0. What do you do?（何する？） はじめの一歩は、ほぼ全く前処理せずに以下の分析フローを一通り流してみました。 1-1. データ取得（読み込み） 1-2. 学習データと検証用データの分割 1-3. データ整形 分析可能にするための最低限のもののみ 1-4. モデルを学習させる 手を抜くためにXGBoostを使う 1-5. 予測してみる 学習させたXGBoostモデルで検証用データに対する予測を行い、成績を確認 1-6. submit（提出） この記事では、とりあえずここまでやりました。 前処理は、本当に何もしてないです。それでも、ものすごいつまづきました。 それでも諦めない(｀･ω･´)根性 1. Source cord 1-0. import modules 今回使ったモジュールです import numpy as np # linear algebra import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv) from sklearn.model_selection import train_test_split from sklearn.metrics import mean_squared_error from xgboost import XGBClassifier, XGBRegressor 1-1. get data （データ取得（読み込み）） import os for dirname, _, filenames in os.walk('/kaggle/input'): for filename in filenames: print(os.path.join(dirname, filename)) # kaggleなので、csvファイルがこのpathになってます train = pd.read_csv('/kaggle/input/house-prices-advanced-regression-techniques/train.csv') test = pd.read_csv('/kaggle/input/house-prices-advanced-regression-techniques/test.csv') train.head() # 処理結果 /kaggle/input/house-prices-advanced-regression-techniques/sample_submission.csv /kaggle/input/house-prices-advanced-regression-techniques/data_description.txt /kaggle/input/house-prices-advanced-regression-techniques/train.csv /kaggle/input/house-prices-advanced-regression-techniques/test.csv Id MSSubClass MSZoning LotFrontage LotArea Street Alley LotShape LandContour Utilities ... PoolArea PoolQC Fence MiscFeature MiscVal MoSold YrSold SaleType SaleCondition SalePrice 0 1 60 RL 65.0 8450 Pave NaN Reg Lvl AllPub ... 0 NaN NaN NaN 0 2 2008 WD Normal 208500 1 2 20 RL 80.0 9600 Pave NaN Reg Lvl AllPub ... 0 NaN NaN NaN 0 5 2007 WD Normal 181500 2 3 60 RL 68.0 11250 Pave NaN IR1 Lvl AllPub ... 0 NaN NaN NaN 0 9 2008 WD Normal 223500 3 4 70 RL 60.0 9550 Pave NaN IR1 Lvl AllPub ... 0 NaN NaN NaN 0 2 2006 WD Abnorml 140000 4 5 60 RL 84.0 14260 Pave NaN IR1 Lvl AllPub ... 0 NaN NaN NaN 0 12 2008 WD Normal 250000 5 rows × 81 columns カラムはどうなってるのかな train.columns # 結果：こんなにたくさんカラムが... Index(['Id', 'MSSubClass', 'MSZoning', 'LotFrontage', 'LotArea', 'Street', .... (略) .... 'SaleCondition', 'SalePrice'], dtype='object') 欠損値も見ておこう train.isnull().sum()[train.isnull().sum() > 0].sort_values(ascending=False) PoolQC 1453 MiscFeature 1406 Alley 1369 Fence 1179 .... (略) 欠損値の確認は、この記事を参考にしました 1-2. 学習データと検証用データの分割 まず、予測したい値（いわゆる目的変数）を分離しておきます。 House Priceの場合は、SalePriceが予測したい値みたいですね。 proto_x = train.drop('SalePrice', axis=1) proto_y = train['SalePrice'] protoというprefixをつけたのは、お試しだからです。 深い意味はないです。 次に、学習に使うデータと、検証(test)に使うデータを分割します x_train, x_test, y_train, y_test = \ train_test_split(proto_x, proto_y, test_size=0.20, random_state=0) 1-3. データ整形 さーて、試しにこれで学習してみようかなぁああ！！(^ワ^*) model = XGBRegressor(n_estimators=20, random_state=0) model.fit(x_train, y_train) 結果...さすがに気が早かったか.... --------------------------------------------------------------------------- ValueError Traceback (most recent call last) <ipython-input-36-1c4e8355ea5e> in <module> 1 model = XGBRegressor(n_estimators=20, random_state=0) ----> 2 model.fit(x_train, y_train) .... (略) .... ValueError: DataFrame.dtypes for data must be int, float or bool. Did not expect the data types in fields MSZoning, Street, Alley, .... (略) .... PoolQC, Fence, MiscFeature, SaleType, SaleCondition model.fit()に渡すDataFrame内の値は、int, float, bool のうちどれかでなければいけないらしい。 この条件に反するfield名は、 MSZoning, Street, ...., SaleConditionということのようだ。 全く前処理しない場合の成績を見てみたかったが、これらのフィールドについては、消すか変換してあげないといけない。 「前処理なし」により近いのは「消す」方かな！(^ワ^*)keshichae!! # 消えてしまえ！ not_expected_type_column_names = [ 'MSZoning', 'Street', 'Alley', 'LotShape', 'LandContour', 'Utilities', 'LotConfig', 'LandSlope', 'Neighborhood', 'Condition1', 'Condition2', 'BldgType', 'HouseStyle', 'RoofStyle', 'RoofMatl', 'Exterior1st', 'Exterior2nd', 'MasVnrType', 'ExterQual', 'ExterCond', 'Foundation', 'BsmtQual', 'BsmtCond', 'BsmtExposure', 'BsmtFinType1', 'BsmtFinType2', 'Heating', 'HeatingQC', 'CentralAir', 'Electrical', 'KitchenQual', 'Functional', 'FireplaceQu', 'GarageType', 'GarageFinish', 'GarageQual', 'GarageCond', 'PavedDrive', 'PoolQC', 'Fence', 'MiscFeature', 'SaleType', 'SaleCondition' ] # _0 というpostfixは、これからたくさんチャレンジするつもりなのでつけてみた x_train_0 = x_train.drop(not_expected_type_column_names, axis=1) x_test_0 = x_test.drop(not_expected_type_column_names, axis=1) x_train_0.head() # 結果: 81列あったデータが37列まで減った。スッキリ！（） Id MSSubClass LotFrontage LotArea OverallQual OverallCond YearBuilt YearRemodAdd MasVnrArea BsmtFinSF1 ... GarageArea WoodDeckSF OpenPorchSF EnclosedPorch 3SsnPorch ScreenPorch PoolArea MiscVal MoSold YrSold 618 619 20 90.0 11694 9 5 2007 2007 452.0 48 ... 774 0 108 0 0 260 0 0 7 2007 870 871 20 60.0 6600 5 5 1962 1962 0.0 0 ... 308 0 0 0 0 0 0 0 8 2009 92 93 30 80.0 13360 5 7 1921 2006 0.0 713 ... 432 0 0 44 0 0 0 0 8 2009 817 818 20 NaN 13265 8 5 2002 2002 148.0 1218 ... 857 150 59 0 0 0 0 0 7 2008 302 303 20 118.0 13704 7 5 2001 2002 150.0 0 ... 843 468 81 0 0 0 0 0 1 2006 5 rows × 37 columns 1-4. モデルを学習させる model = XGBRegressor(n_estimators=20, random_state=0) model.fit(x_train_0, y_train) # 結果: なんか出たぽよ～ XGBRegressor(base_score=0.5, booster=None, colsample_bylevel=1, colsample_bynode=1, colsample_bytree=1, gamma=0, gpu_id=-1, importance_type='gain', interaction_constraints=None, learning_rate=0.300000012, max_delta_step=0, max_depth=6, min_child_weight=1, missing=nan, monotone_constraints=None, n_estimators=20, n_jobs=0, num_parallel_tree=1, objective='reg:squarederror', random_state=0, reg_alpha=0, reg_lambda=1, scale_pos_weight=1, subsample=1, tree_method=None, validate_parameters=False, verbosity=None) 1-5. 予測してみる # 予測結果の出力 # 学習に使ったデータに対しての予測 predict_result_for_tr_data = model.predict(x_train_0) # 検証用データに対しての予測 predict_result = model.predict(x_test_0) 結果は以下の通り。 House Priceという課題では、RMSE(Root Mean Squared Error)という指標で成績を評価するらしいので、RMSEを計算する。 # 学習に使ったデータに対してモデルを適用した場合の誤差 rmse_of_train = np.sqrt(mean_squared_error(predict_result_for_tr_data, y_train)) rmse_of_train # 結果 >> 9283.227477186625 # 検証用データに対してモデルを適用した場合の誤差 rmse_of_test = np.sqrt(mean_squared_error(predict_result, y_test)) rmse_of_test # 結果 >> 37010.13476747102 test も train もそうですが、RMSE値が小さい値になるほど、より正解に近い予測ができているといえるようです。 なので、上記「9283.22..., 37010.13...」よりも小さい値を目指す感じですね(｀･ω･´) 1-6 submit 一応提出してみます。 kaggleの点数を上げたいわけではないのですが、どのくらい実力が上がっているのか、後で参考になると思うので、最低ラインの場合の実力を確認。 ここは特に解説しませんが、コードは載せておきます。 # trainデータに対して実施したのと同じ手順で、1-3 ~ 1-4を実施します # 1-3. データ整形 not_expected_type_column_names = [ 'MSZoning', 'Street', 'Alley', 'LotShape', 'LandContour', 'Utilities', 'LotConfig', 'LandSlope', 'Neighborhood', 'Condition1', 'Condition2', 'BldgType', 'HouseStyle', 'RoofStyle', 'RoofMatl', 'Exterior1st', 'Exterior2nd', 'MasVnrType', 'ExterQual', 'ExterCond', 'Foundation', 'BsmtQual', 'BsmtCond', 'BsmtExposure', 'BsmtFinType1', 'BsmtFinType2', 'Heating', 'HeatingQC', 'CentralAir', 'Electrical', 'KitchenQual', 'Functional', 'FireplaceQu', 'GarageType', 'GarageFinish', 'GarageQual', 'GarageCond', 'PavedDrive', 'PoolQC', 'Fence', 'MiscFeature', 'SaleType', 'SaleCondition' ] test_proto = test.drop(not_expected_type_column_names, axis=1) # 1-4. モデル学習 # モデルはさっき学習させたもの(model)を使います！ # 1-5. 予測 predict_proto = model.predict(test_proto) submission = pd.DataFrame({'Id': test['Id'], 'SalePrice': predict_proto}) submission.to_csv('submission_new.csv', index=False) なんか、何もしてないのに上位68%に入った...XGBoostしゅごい...(*´﹃｀*) 2. Result （結果と所感） これが正真正銘の最低ライン。完全に手抜きした（完全に何もしなかった）場合の結果。 ここからどこまで改善できるか。 個人的には、今回試した6ステップを、何も見ずにそらでできるようになるのが、まず最初のステップと捉えています。 結果を出すためには細かい手順も必要なのですが、それらは結局課題によって異なるので、まずはなるべく幅広い課題に使えるものから試しつつ覚えていきたいと思います。 Next step 次回以降は、 カテゴリ変数を数値化 カテゴリ変数をXGBoostで使えるようにする pd.get_dummies()で簡単にできるみたいです 欠損値の穴埋め あたりをやってみて、成績が多少変わるかどうか試してみようと思います。 ハイパーパラメータの探索や、他のモデルのお試しなどはもっと先になると思いますが、やってみたいと思います。 ド素人なので、ひとまずはこの課題で上位50パーセント以内に入ることを目指します ( *˙︶˙*)وｸﾞｯ! Reference（参考にしたもの） この本は、様々な手法を広く浅く辞書的に掲載した本になっている気がします。 コードを書きながら1歩1歩学んでいく系の本ではないです。 なので、前から順に読んでいっても、全然書けるようにならないです(´；ω；｀)ﾌﾞﾜｯ ただ、初心者向けなのか、私レベルの人でもとても理解しやすかったです！ さいごに.... もうちょっと、もうちょっと情報量の少ない記事にしたい... できれば3~5スクロールで最後まで見える記事がいい.... あ、喋るの好きだからムリかな....(*´﹃｀*)

はじめに カテゴリ変数を数値化する方法にはいろいろな種類があります。 一番有名な方法としては、One-Hotエンコーディングやラベルエンコーディングなどがあります。 勉強会を聞いているときにFeature Hashing(Hashing, Hashed Featureなど)という手法がありましたのでそちらについて、まとめていきたいと思います。 また、それ以外のカテゴリエンコーディングはこちらの記事を確認いただけると、とても分かりやすくまとめられております。 本記事でもこちらの記事を参考にFeature Hashingをやっていきたいと思います。 Feature Hashingの手法について 手法についてはこちらの記事に乗っていましたので、利用しました。 例えば以下のようなデータがあり、Feature Hashingをしようとします。 趣味 特技 Aさん 散歩 マラソン Bさん 読書 勉強 One-Hotエンコーディングをすると以下のようになります。 Aさん = \begin{pmatrix} 趣味:散歩 \\ 趣味:読書 \\ 特技:マラソン \\ 特技:勉強 \\ \end{pmatrix} = \begin{pmatrix} 1 \\ 0 \\ 1 \\ 0 \\ \end{pmatrix} Bさん = \begin{pmatrix} 趣味:散歩 \\ 趣味:読書 \\ 特技:マラソン \\ 特技:勉強 \\ \end{pmatrix} = \begin{pmatrix} 0 \\ 1 \\ 0 \\ 1 \\ \end{pmatrix} Feature Hashingは以下のアルゴリズムで行われます。 Aさんをカテゴリ変数にしたいと思います。 1. 初期化 Feature Hashingしたあとの特徴量の次元数をmとする。 今回はFeature Hashingを利用して、4次元の特徴量(m=4)にする。 特徴1 特徴2 特徴3 特徴4 Aさん 0 0 0 0 2. ハッシュ値を計算する $h=hash("趣味:散歩")$を計算する import hashlib dat = '散歩' # hash(int)を求める hs = hash(dat) # 出力 # 4893028045998397989 3. modを求める 値をハッシュ化し、整数に変換し、m(列数)で割った余りを求める。 今回はmod(hs, 4)で求められます。 x = divmod(hs, 4) print(x) # 出力 # (1223257011499599497, 1) 4. modの余りに対応した列に+1する 今回の例では、1余りとなっているので特徴量2を+1します。 (特徴量は1-4ですが、余りは0-3になるので0が特徴量1に対応します) 特徴1 特徴2 特徴3 特徴4 Aさん 0 1 0 0 1から4までをFeature Hashingするカラムすべてに対して行う。 # Aさんの特技:マラソンについても計算 import hashlib dat = 'マラソン' hs = hash(dat) print(hs) # 出力 # 2068416709068501704 x = divmod(hs, 4) print(x) # 出力 # (517104177267125426, 0) # 特徴量1を+1する 変換前 趣味 特技 Aさん 散歩 マラソン 変換後 特徴1 特徴2 特徴3 特徴4 Aさん 1 1 0 0 このようになりました。 注意しておくべきこと Feature Hashingはカテゴリ変数を好きな次元数mに変換することができます。 One-Hotエンコーディングのように次元数が多くなってしまうのを避けることができます。 また、One-Hotエンコーディングは疎なベクトルになりますが、mを小さくすることで密なベクトルにすることができます。 しかし、問題もあるかと思います。 単純にカテゴリ変数の文字コードをハッシュ化して整数値にしていることです。 文字コードが似ているために、似たベクトルになってしまい機械的には特徴同士にあるはずのない関係が生まれる可能性があります。 使いどころは限られてしまうとは思いますが、カテゴリ変数の次元が多すぎる場合の手段の一つにはなるかと思います。 Pythonでの実装 Pythonで利用する場合は、Category Encodersというライブラリを使うと簡単にできます。 コードはCategory Encodersのすゝめ【AI道場「Kaggle」への道 by 日経 xTECH ビジネスAI① Advent Calendar 2019 10日目】を利用させていただきます。 こちらにはほかのカテゴリエンコーディングの例がわかりやすく説明されています。 今回はGoogle Colabratoryで動かします。 まずはデータを用意します。 !pip install featuretools import featuretools as ft import numpy as np import pandas as pd data = ft.demo.load_mock_customer() df = data['sessions'] df = df[['customer_id', 'device']].head(12) df[7:8]['device'] = 'laptop' # Target Encodingの検証用にカラムを追加 np.random.seed(seed=3) df['target'] = np.random.randint(1, 10, size=len(df)) deviseという列を利用して、Feature Hashingをします。 実際にFeature Hashingするコードは以下となります。 また、次元数を指定しない場合はm=8(デフォルト値)となります。 !pip install category_encoders import category_encoders as ce cate_col = 'device' he = ce.HashingEncoder(cols=cate_col) he_df = he.fit_transform(df[cate_col]) pd.concat([df[cate_col], he_df], axis=1) 公式ドキュメントにはn_componentsという引数で指定ができます。 n_components=4と指定した場合は以下のようになります。 he = ce.HashingEncoder(cols=cate_col, drop_invariant=True, n_components=4) he_df = he.fit_transform(df[cate_col]) pd.concat([df[cate_col], he_df], axis=1) 4次元になっています。 また、行列をみるとすべてが0の列があることが確認できます。 例えば、col_7はすべて0の列です。 そのような列を削除することもできます。 he = ce.HashingEncoder(cols=cate_col, drop_invariant=True) he_df = he.fit_transform(df[cate_col]) pd.concat([df[cate_col], he_df], axis=1) さいごに 今回は、ほかの記事をまとめたものになりました。 Feature Hashingを調べていて、どうも内容がバラバラしておりわかりづらかったためまとめています。 使いどころは難しいですが、手段の一つとして知っておくのはアリだと感じました。 参考文献

はじめに カテゴリ変数を数値化する方法にはいろいろな種類があります。 一番有名な方法としては、One-Hotエンコーディングやラベルエンコーディングなどがあります。 勉強会を聞いているときにFeature Hashing(Hashing, Hashed Featureなど)という手法がありましたのでそちらについて、まとめていきたいと思います。 また、それ以外のカテゴリエンコーディングはこちらの記事を確認いただけると、とても分かりやすくまとめられております。 本記事でもこちらの記事を参考にFeature Hashingをやっていきたいと思います。 Feature Hashingの手法について 手法についてはこちらの記事に乗っていましたので、利用しました。 例えば以下のようなデータがあり、Feature Hashingをしようとします。 趣味 特技 Aさん 散歩 マラソン Bさん 読書 勉強 One-Hotエンコーディングをすると以下のようになります。 Aさん = \begin{pmatrix} 趣味:散歩 \\ 趣味:読書 \\ 特技:マラソン \\ 特技:勉強 \\ \end{pmatrix} = \begin{pmatrix} 1 \\ 0 \\ 1 \\ 0 \\ \end{pmatrix} Bさん = \begin{pmatrix} 趣味:散歩 \\ 趣味:読書 \\ 特技:マラソン \\ 特技:勉強 \\ \end{pmatrix} = \begin{pmatrix} 0 \\ 1 \\ 0 \\ 1 \\ \end{pmatrix} Feature Hashingは以下のアルゴリズムで行われます。 Aさんをカテゴリ変数にしたいと思います。 1. 初期化 Feature Hashingしたあとの特徴量の次元数をmとする。 今回はFeature Hashingを利用して、4次元の特徴量(m=4)にする。 特徴1 特徴2 特徴3 特徴4 Aさん 0 0 0 0 2. ハッシュ値を計算する $h=hash("趣味:散歩")$を計算する import hashlib dat = '散歩' # hash(int)を求める hs = hash(dat) # 出力 # 4893028045998397989 3. modを求める 値をハッシュ化し、整数に変換し、m(列数)で割った余りを求める。 今回はmod(hs, 4)で求められます。 x = divmod(hs, 4) print(x) # 出力 # (1223257011499599497, 1) 4. modの余りに対応した列に+1する 今回の例では、1余りとなっているので特徴量2を+1します。 (特徴量は1-4ですが、余りは0-3になるので0が特徴量1に対応します) 特徴1 特徴2 特徴3 特徴4 Aさん 0 1 0 0 1から4までをFeature Hashingするカラムすべてに対して行う。 # Aさんの特技:マラソンについても計算 import hashlib dat = 'マラソン' hs = hash(dat) print(hs) # 出力 # 2068416709068501704 x = divmod(hs, 4) print(x) # 出力 # (517104177267125426, 0) # 特徴量1を+1する 変換前 趣味 特技 Aさん 散歩 マラソン 変換後 特徴1 特徴2 特徴3 特徴4 Aさん 1 1 0 0 このようになりました。 注意しておくべきこと Feature Hashingはカテゴリ変数を好きな次元数mに変換することができます。 One-Hotエンコーディングのように次元数が多くなってしまうのを避けることができます。 また、One-Hotエンコーディングは疎なベクトルになりますが、mを小さくすることで密なベクトルにすることができます。 しかし、問題もあるかと思います。 単純にカテゴリ変数の文字コードをハッシュ化して整数値にしていることです。 文字コードが似ているために、似たベクトルになってしまい機械的には特徴同士にあるはずのない関係が生まれる可能性があります。 使いどころは限られてしまうとは思いますが、カテゴリ変数の次元が多すぎる場合の手段の一つにはなるかと思います。 Pythonでの実装 Pythonで利用する場合は、Category Encodersというライブラリを使うと簡単にできます。 コードはCategory Encodersのすゝめ【AI道場「Kaggle」への道 by 日経 xTECH ビジネスAI① Advent Calendar 2019 10日目】を利用させていただきます。 こちらにはほかのカテゴリエンコーディングの例がわかりやすく説明されています。 今回はGoogle Colabratoryで動かします。 まずはデータを用意します。 !pip install featuretools import featuretools as ft import numpy as np import pandas as pd data = ft.demo.load_mock_customer() df = data['sessions'] df = df[['customer_id', 'device']].head(12) df[7:8]['device'] = 'laptop' # Target Encodingの検証用にカラムを追加 np.random.seed(seed=3) df['target'] = np.random.randint(1, 10, size=len(df)) deviseという列を利用して、Feature Hashingをします。 実際にFeature Hashingするコードは以下となります。 また、次元数を指定しない場合はm=8(デフォルト値)となります。 !pip install category_encoders import category_encoders as ce cate_col = 'device' he = ce.HashingEncoder(cols=cate_col) he_df = he.fit_transform(df[cate_col]) pd.concat([df[cate_col], he_df], axis=1) 公式ドキュメントにはn_componentsという引数で指定ができます。 n_components=4と指定した場合は以下のようになります。 he = ce.HashingEncoder(cols=cate_col, drop_invariant=True, n_components=4) he_df = he.fit_transform(df[cate_col]) pd.concat([df[cate_col], he_df], axis=1) 4次元になっています。 また、行列をみるとすべてが0の列があることが確認できます。 例えば、col_7はすべて0の列です。 そのような列を削除することもできます。 he = ce.HashingEncoder(cols=cate_col, drop_invariant=True) he_df = he.fit_transform(df[cate_col]) pd.concat([df[cate_col], he_df], axis=1) さいごに 今回は、ほかの記事をまとめたものになりました。 Feature Hashingを調べていて、どうも内容がバラバラしておりわかりづらかったためまとめています。 使いどころは難しいですが、手段の一つとして知っておくのはアリだと感じました。 参考文献

概要 深層学習を用いた、単一画像における超解像手法であるSRCNNの実装したので、それのまとめの記事です。 Python + Tensorflow(Keras)で実装を行いました。 今回は、2倍拡大の超解像にチャレンジしました。 以前、3つに分けて記事を上げていたのですが、コードが汚かった＋1つにまとめたかったので再投稿しています。ご了承ください。 今回紹介するコードはGithubにも載せています。 1. 超解像のおさらい 超解像について簡単に説明をします。 超解像とは解像度の低い画像に対して、解像度を向上させる技術のことです。 ここでいう解像度が低いとは、画素数が少なかったり、高周波成分（輪郭などの鮮鋭な部分を表す成分）がないような画像のことです。 以下の図で例を示します。(図は[論文]より引用) (a)は原画像、(b)は画素数の少ない画像を見やすいように原画像と同じ大きさにした画像、(c)は高周波成分を含まない画像の例です。 (b)と(c)は、荒かったりぼやけていたりしていると思います。 このような状態を解像度が低い画像といいます。 そして、超解像はこのような解像度が低い画像に処理を行い、(a)のような精細な画像を出力することを目的としています。 2. SRCNNの超解像アルゴリズム SRCNNは初めて超解像分野に深層学習を導入したモデルです。 モデルとしては単純な構造にはなっていて、Convolutionを3層組み合わせたモデルとなっています。 SRCNNのアルゴリズムの概要図は以下の通りです。(図は論文から引用) 今回実装したSRCNNは、事前にbicubicで拡大処理を行います。 3. 実装したアルゴリズム 今回、実装したSRCNNのモデルは以下のように組みました。(コードの一部を抽出) def SRCNN(input_channels = 1): #input single image input_shape = Input((None, None, input_channels)) #convolution conv2d_0 = Conv2D(filters = 64, kernel_size = (9, 9), padding = "same", activation = "relu")(input_shape) conv2d_1 = Conv2D(filters = 32, kernel_size = (1, 1), padding = "same", activation = "relu")(conv2d_0) conv2d_2 = Conv2D(filters = input_channels, kernel_size = (5, 5), padding = "same")(conv2d_1) model = Model(inputs = input_shape, outputs = [conv2d_2]) model.summary() return model 前述の通り、3層のConvolution層を組み合わせています。 4. 使用したデータセット 今回は、データセットにDIV2K datasetを使用しました。 このデータセットは、単一画像のデータセットで、学習用が800種、検証用とテスト用が100種類ずつのデータセットです。 今回は、学習用データと検証用データを使用しました。 パスの構造はこんな感じです。 train_sharp - 0001.png - 0002.png - ... - 0800.png val_sharp - 0801.png - 0802.png - ... - 0900.png このデータをbicubicで縮小したりしてデータセットを生成しました。 5. 画像評価指標PSNR 今回は、画像評価指標としてPSNRを使用しました。 PSNR とは Peak Signal-to-Noise Ratio(ピーク信号対雑音比) の略で、単位はデジベル (dB) で表せます。 PSNR は信号の理論ピーク値と誤差の2乗平均を用いて評価しており、8bit画像の場合、255(最大濃淡値)を誤差の標準偏差で割った値です。 今回は、8bit画像を使用しましたが、計算量を減らすため、全画素値を255で割って使用しました。 そのため、最小濃淡値が0で最大濃淡値が1です。 dB値が高いほど拡大した画像が元画像に近いことを表します。 PSNRの式は以下のとおりです。 PSNR = 10\log_{10} \frac{1^2 * w * h}{\sum_{x=0}^{w-1}\sum_{y=0}^{h-1}(p_1(x,y) - p_2(x,y))^2 } なお、$w$は画像の幅、$h$は画像の高さを表しており、$p_1$は元画像、$p_2$はPSNRを計測する画像を示しています。 6. コードの使用方法 このコード使用方法は、自分が執筆した別の実装記事とほとんど同じです。 ① 学習データ生成 まず、Githubからコードを一式ダウンロードして、カレントディレクトリにします。 Windowsのコマンドでいうとこんな感じ。 C:~/keras_SRCNN> 次に、main.pyから生成するデータセットのサイズ・大きさ・切り取る枚数、ファイルのパスなどを指定します。 main.pyの12~21行目です。 使うPCのメモリ数などに応じで、画像サイズや学習データ数の調整が必要です。 main.py parser.add_argument('--train_height', type = int, default = 33, help = "Train data size(height)") parser.add_argument('--train_width', type = int, default = 33, help = "Train data size(width)") parser.add_argument('--test_height', type = int, default = 700, help = "Test data size(height)") parser.add_argument('--test_width', type = int, default = 700, help = "Test data size(width)") parser.add_argument('--train_dataset_num', type = int, default = 10000, help = "Number of train datasets to generate") parser.add_argument('--test_dataset_num', type = int, default = 5, help = "Number of test datasets to generate") parser.add_argument('--train_cut_num', type = int, default = 10, help = "Number of train data to be generated from a single image") parser.add_argument('--test_cut_num', type = int, default = 1, help = "Number of test data to be generated from a single image") parser.add_argument('--train_path', type = str, default = "../../dataset/DIV2K_train_HR", help = "The path containing the train image") parser.add_argument('--test_path', type = str, default = "../../dataset/DIV2K_valid_HR", help = "The path containing the test image") 指定したら、コマンドでデータセットの生成をします。 C:~/keras_SRCNN>python main.py --mode train_datacreate これで、train_data_list.npzというファイルのデータセットが生成されます。 ついでにテストデータも同じようにコマンドで生成します。コマンドはこれです。 C:~/keras_SRCNN>python main.py --mode test_datacreate ② 学習 次に学習を行います。 設定するパラメータの箇所は、epoch数と学習率、今回のモデルの層の数です。 まずは、main.pyの22~25行目 main.py parser.add_argument('--learning_rate', type = float, default = 1e-4, help = "Learning_rate") parser.add_argument('--BATCH_SIZE', type = int, default = 32, help = "Training batch size") parser.add_argument('--EPOCHS', type = int, default = 500, help = "Number of epochs to train for") 後は、学習のパラメータをあれこれ好きな値に設定します。73~86行目です。 main.py train_model = model.SRCNN() optimizers = tf.keras.optimizers.Adam(lr = args.learning_rate) train_model.compile(loss = "mean_squared_error", optimizer = optimizers, metrics = [psnr]) train_model.fit(train_x, train_y, epochs = args.EPOCHS, verbose = 2, batch_size = args.BATCH_SIZE) train_model.save("SRCNN_model.h5") optimizerはmomentum、損失関数はmean_squared_errorを使用しています。 学習はデータ生成と同じようにコマンドで行います。 C:~/keras_SRCNN>python main.py --mode train_model これで、学習が終わるとモデルが出力されます。 ③ 評価 最後にモデルを使用してテストデータで評価を行います。 これも同様にコマンドで行いますが、事前に①でテストデータも生成しておいてください。 C:~/keras_SRCNN>python main.py --mode evaluate このコマンドで、画像を出力してくれます。 7. 結果 出力した画像はこのようになりました。 なお、今回は輝度値のみで学習を行っているため、カラー画像には対応していません。 元画像 低解像度画像 PSNR:30.83 生成画像 PSNR:31.43 論文ほどの数値ではありませんが、補間法よりはPSNRの値は大きくなりました。 論文の数値に近づけるには、学習回数などのパラメータ設定・学習データ数の増加などを行えばよさそうです。 最後に元画像・低解像度画像・生成画像の一部を並べて表示してみます。 数値の差があまりないので目視ではあまり変化は見えないですね... 8. コードの全容 前述の通り、Githubに載せています。 pythonのファイルは主に3つあります。 各ファイルの役割は以下の通りです。 data_create.py : データ生成に関するコード。 model.py : 超解像のアルゴリズムに関するコード。 main.py : 主に使用するコード。 9. まとめ 今回は、以前作っていたSRCNNのコードをもう一度見直し、記事も見直しました。 記事が長くなってしまいましたが、最後まで読んでくださりありがとうございました。 参考文献 ・Image Super-Resolution Using Deep Convolutional Networks 　実装の参考にした論文。 ・画素数の壁を打ち破る 複数画像からの超解像技術 　超解像の説明のために使用。 ・DIV2K dataset 　今回使用したデータセット。

はじめに レコメンドエンジンの仕組みってよくわかりませんよね。仕組みを知ろうと思って説明を見てみても複雑な数式とか難しそうな言葉がたくさん出てきたり、そもそも日本語じゃなかったり... 僕もレコメンドエンジンについて何も知らない状態から色々と調べました。それでほんの少しだけレコメンドエンジンについてわかるようになったので、ここに記録を残しておこうと思いました。 この記事の対象者 レコメンドエンジンの名前くらいしかまだ知らない人 レコメンドエンジンについて少し興味がある人 レコメンドエンジンは「おすすめシステム」のこと YoutubeやSpotifyのように、何かコンテンツを提供するサービスで例えます。 コンテンツをユーザーに提供する際に、どのコンテンツが気に入られるのかというのはもちろんユーザーによって異なります。 そのため、ただソートされただけのコンテンツやただ人気度の高いコンテンツを提供するだけでは、それぞれのユーザーにとって自分の欲しいコンテンツが表示される機会が少なくなってしまいます。 そこで最近の大手サービスにはレコメンドエンジンが搭載されています。優秀なエンジニアたちは高性能なレコメンドエンジンをサービスに搭載することで、ユーザーひとりひとりがより興味を持つようなコンテンツを表示してサービスの魅力を向上させています。 レコメンドエンジンの開発自体は難易度が高く、作り方も様々です。今回の記事ではThompson Samplingという手法を使って簡易エンジンの開発をしてみます。 レコメンドの仕組み 典型的なレコメンドの仕組みとして、2つの例を挙げることができます。 協調フィルタリング コンテンツベースのフィルタリング もしくは、これらのハイブリッドの手法が用いられたりもします。 協調フィルタリング 協調フィルタリングは、「オススメをしたいユーザーに似たようなユーザーの好みのコンテンツをレコメンドする」仕組みです。 コンテンツベースのフィルタリング 一方コンテンツベースのフィルタリングは、「オススメをしたいユーザーの過去の行動履歴からそのユーザーの好きそうなコンテンツをレコメンドする」仕組みです。 両者の欠点 これら２つともよく出来た仕組みなのですが、両者共通の欠点として「ユーザーの情報が集まらなければレコメンドの効果があまり出せない」ことが挙げられます。 この欠点のことをコールドスタート問題と呼んだりもします。 このコールドスタート問題をなんとかするために、今から私たちは多腕バンディット問題について考えていきます。 多腕バンディット問題 多腕バンディット問題について考えるとき、よくスロットの例が挙げられます。 なので僕もスロットの例を使おうと思います。 多腕バンディット問題 ここに4つのスロットマシンがあります。このスロットマシンを100回まで引けるとき、どのようにすれば報酬を最大化できるでしょうか。 この問題を解く方法は複数提案されています。 ε-greedyという手法は比較的イメージがしやすいと思います。「決められた数だけスロットマシンを試してみて、もっとも結果の良かったマシンを残りの回数分ずっと引き続ける」というものです。よって最も結果の良かったマシンを引ける回数は100-4n回となります。 しかしこの解法にはジレンマがあります。最初に検証のために4n回スロットマシンを回すのですが、検証の回数が多いと結果の良いマシンを回す回数が減ってしまい、逆に検証の回数が少ないと結果の良いマシンを選ぶことができる確率が下がってしまいます。 このように難しい多腕バンディット問題ですが、今回は先程も名前が出てきたThompson Samplingという方法をPythonで実行していきます。 Thompson Sampling Thompson Samplingでは3つの段階に分けて処理が行われます。 それぞれのスロットマシンが当たる確率を確率分布から推定する 確率分布の値が最大であるスロットマシンを引く 結果から、スロットマシンの当たる確率を推定する このThompson Samplingは昔からある手法なのですが、当たる確率が高いことでも知られています。 実装 Thompson Samplingを用いて、今回はゲームを推薦するシステムを試しに実装してみます。 ※この記事を非常に参考にしています！ 前提 Python3系がインストールされているものとします。 ライブラリは必要に応じてinstallしてもらえれば大丈夫です。 動作確認用CSVの自動生成 id ユーザーid ゲーム名 ユーザーのアクションから得られたreward というCSVデータが仮の行動履歴として保存されるようにします。 rewardの列は、ユーザーがそのゲームに対して何かプラスのアクションをしたことにより記録される値とします。 makecsv.py import csv import random records = 10000 users = 250 print("Making %d records...", records) fieldnames = ['id', 'user_id', 'game_name', 'gain'] writer = csv.DictWriter(open("data.csv", "w"), fieldnames=fieldnames) user_id_list = range(1, users) game_names = ["Fortnite", "Valorant", "PUBG", "Clash Royale", "Apex Legends", "FIFA", "Shadowverse", "CSGO", "CoD"] gains = [1, 5, 7] writer.writerow(dict(zip(fieldnames, fieldnames))) for i in range(0, records): writer.writerow(dict([ ('id', i), ('user_id', random.choice(user_id_list)), ('game_name', random.choice(game_names)), ('gain', random.choice(gains)) ])) 250人のユーザー 9種類のゲーム 10000行生成 という設定のプログラムです。このプログラムを実行すると10000行の行動履歴CSVデータが生成されます。 Thompson Sampling用クラスの実装 次はレコメンドエンジン本体となるクラスを実装します。 まずは親クラスからです。 thompson_sampling_replayer.py import numpy as np from tqdm import tqdm class ReplaySimulator(object): ''' Class to provide base functionality for simulatoing the replayer method for online algorithms. ''' def __init__(self, n_visits, reward_history, item_col_name, visitor_col_name, reward_col_name, n_iterations=1, random_seed=1): np.random.seed(random_seed) self.reward_history = reward_history self.item_col_name = item_col_name self.visitor_col_name = visitor_col_name self.reward_col_name = reward_col_name # number of visits to replay/simulate self.n_visits = n_visits # number of runs to average over self.n_iterations = n_iterations # items under test self.items = self.reward_history[self.item_col_name].unique() self.n_items = len(self.items) # visitors in the historical reward_history (e.g., from gain df) self.visitors = self.reward_history[self.visitor_col_name].unique() self.n_visitors = len(self.visitors) def reset(self): # number of times each item has been sampled self.n_item_samples = np.zeros(self.n_items) # fraction of time each item has resulted in a reward self.n_item_rewards = np.zeros(self.n_items) def replay(self): results = [] for iteration in tqdm(range(0, self.n_iterations)): self.reset() total_rewards = 0 fraction_relevant = np.zeros(self.n_visits) for visit in range(0, self.n_visits): found_match = False while not found_match: # choose a random visitor visitor_idx = np.random.randint(self.n_visitors) visitor_id = self.visitors[visitor_idx] # select an item to offer the visitor item_idx = self.select_item() item_id = self.items[item_idx] reward = self.reward_history.query('{} == @item_id and {} == @visitor_id'.format(self.item_col_name, self.visitor_col_name))[self.reward_col_name] found_match = reward.shape[0] > 0 reward_value = reward.iloc[0] self.record_result(visit, item_idx, reward_value) ## record metrics total_rewards += reward_value fraction_relevant[visit] = total_rewards * 1. / (visit + 1) result = {} result['iteration'] = iteration result['visit'] = visit result['item_id'] = item_id result['visitor_id'] = visitor_id result['reward'] = reward_value result['total_reward'] = total_rewards result['fraction_relevant'] = total_rewards * 1. / (visit + 1) results.append(result) return results def select_item(self): game_names = ["Fortnite", "Valorant", "PUBG", "Clash Royale", "Apex Legends", "FIFA", "Shadowverse", "CSGO", "CoD"] return game_names[np.random.randint(self.n_items)] def record_result(self, visit, item_idx, reward): self.n_item_samples[item_idx] += 1 alpha = 1. / self.n_item_samples[item_idx] self.n_item_rewards[item_idx] += alpha * (reward - self.n_item_rewards[item_idx]) 今回不必要な処理もありますが、僕はこのクラスを親クラスとして使っています。 こちらはThompson Samplingをするためのクラスです。 thompson_sampling_replayer.py class ThompsonSamplingReplayer(ReplaySimulator): ''' Class to provide functionality for simulating the replayer method on a Thompson Sampling bandit algorithm. ''' def reset(self): self.alphas = np.ones(self.n_items) self.betas = np.ones(self.n_items) def select_item(self): samples = [np.random.beta(a, b) for a, b in zip(self.alphas, self.betas)] return np.argmax(samples) def record_result(self, visit, item_idx, reward): ## update value estimate if reward == 1: self.alphas[item_idx] += 1 else: self.betas[item_idx] += 1 先程の親クラスを継承し、Thompson Samplingの処理を入れてあります。 実行 ではこのクラスを実行してレコメンドをさせてみます。 from appのところのディレクトリ構成だけ自分の環境に合わせて変えてください。 recommender.py import numpy as np import pandas as pd from app import ThompsonSamplingReplayer gain_df = pd.read_csv('./output.csv') reward_threshold = 5 gain_df['reward'] = gain_df.eval('gain >= @reward_threshold').astype(int) n_visits = 250 n_iterations = 20 reward_history = gain_df item_col_name = 'game_name' visitor_col_name = 'user_id' reward_col_name = 'reward' thompson_results = ThompsonSamplingReplayer(n_visits, reward_history, item_col_name, visitor_col_name, reward_col_name, n_iterations=n_iterations).replay() thompson_results_df = pd.DataFrame(thompson_results) thompson_results_df.head() thompson_results_df.to_csv('result/thompson_sampling.csv') このファイルを実行すると、thompson_sampling.csvという名前のファイルが生成されています。 そしてそこにはこんな感じのデータが載っています。 ,iteration,visit,item_id,visitor_id,reward,total_reward,fraction_relevant 0,0,0,DOTA2,888,1,1,1.0 1,0,1,CoD,436,1,2,1.0 2,0,2,Pokemon,541,1,3,1.0 3,0,3,DOTA2,906,0,3,0.75 4,0,4,Clash Royale,1096,0,3,0.6 ループ数が少ないのでそこまで信憑性は高くないのですが、それぞれのユーザーに対するレコメンドの結果が出ていると思います。 さいごに レコメンドエンジンは概念も難しく理解しきれていないところもありますが、学習を進めていこうと思います。 見てくださった方ありがとうございました。 参考 https://developers.cyberagent.co.jp/blog/archives/25099/ Icons made by Nikita Golubev from www.flaticon.com https://blog.insightdatascience.com/multi-armed-bandits-for-dynamic-movie-recommendations-5eb8f325ed1d

複数のcsvファイルを選択した場合は、シート別（ファイル名）にデータの貼り付けとグラフ（散布図）の描画を行う。 本コードは17行目(LIST_TARGET)に示す３つの系列のデータがcsvファイルに存在することを前提としている。 データ系列数が３つ以上の場合でもグラフに凡例が表示されなくなるが、動作はするはず。 #CSVのデータをExcelに貼り付けて散布図を作成 import tkinter as tk import tkinter.filedialog as fl import tkinter.messagebox as mb import numpy as np import openpyxl as px from openpyxl.chart import ScatterChart,Reference,Series import os root = tk.Tk() #ファイルダイアログを表示するディレクトリ INI_DIR = "D:/Python/Hello" #書き込み先のExcelの一行目に読み込み元ファイルがどんな物理量かを記載するためのリスト LIST_TARGET = ["時間[sec]", "角度[deg]", "角速度[deg/sec]"] #GUIのボタンを押した時の処理 def get(*args): #書き込み＆グラフ描画用のExcelを新規作成 wb = px.Workbook() #CSVファイルを選択する ※複数選択可 filetype = [("csv file",".csv")] path = fl.askopenfilenames(initialdir=INI_DIR, filetypes=filetype, title="select file") # 読み取った複数のCSVファイルを1つのExcelファイルにまとめる for i in range(len(path)): read_data = np.loadtxt(fname=path[i], skiprows=2, unpack=True, encoding="UTF-8", delimiter=",", dtype="float") #データ数、系列数を変数に格納 series_num = read_data.shape[0] data_num = read_data.shape[1] #ファイル名をシート名に設定 basename = os.path.basename(path[i]) #Excel内に新しいシートを作成してファイル名をつける new_sheet = wb.create_sheet(index=i, title=basename) ws = wb[basename] #Excelの1行目に上記リストを記載する for i in range(len(LIST_TARGET)): ws.cell(row=1, column=i+1, value=LIST_TARGET[i]) #Excelの2行目以降にCSVデータを代入する for i in range(data_num): for j in range(series_num): ws.cell(row=i+2, column=j+1, value=read_data[j][i]) #グラフとして散布図を作成する chart = ScatterChart() #グラフのタイトル chart.title = "時間 vs 角度、角速度" #x軸タイトル追加 chart.x_axis.title = "[sec]" #y軸タイトル追加 chart.y_axis.title = "[deg], [deg/sec]" #x軸とするデータ x_values = Reference(ws, min_col=1, min_row=2, max_row=data_num+1) for i in range(2, series_num+1): #y軸とするデータ y_values = Reference(ws, min_col=i, min_row=1, max_row=data_num+1) #x軸とy軸の組み合わせを決める con = Series(y_values, x_values, title_from_data=True) #マーカーの形 con.marker.symbol = "circle" #マーカーの塗りつぶし con.marker.graphicalProperties.solidFill = "ffffff" #マーカーの枠線の色 con.marker.graphicalProperties.line.solidFill = "000000" #散布図に系列を追加 chart.series.append(con) #シートのA1セルにグラフを追加 ws.add_chart(chart, "A1") #Excelデータをファイル名を指定して保存 typ = [("Excel", ".xlsx")] savefile = fl.asksaveasfilename(title="保存ファイル名を指定", filetypes=typ, initialdir=INI_DIR) wb.save(savefile + ".xlsx") #処理が完了したことをユーザーに知らせる mb.showinfo("確認","Excelファイルの作成が完了しました") message["text"] = "処理が完了しました" message = tk.Label(root,text="ファイルを選択してください", width=30) message.grid(row=0, column=0) button = tk.Button(text="開く", command=get) button.grid(row=0, column=1) root.mainloop()

この記事の概要 FFTの原理をなるべく数式で書いて、Juliaで実装します。 また、この記事を書くに当たって以下の本、記事をめちゃめちゃ参考にしました。 物理数学の直観的方法 (https://bookclub.kodansha.co.jp/product?item=0000194699 ) EMANの物理学（https://eman-physics.net/math/fourier01.html ） 安居院猛, 中嶋正之, 『FFT の使い方』, 産報出版株式会社 (1982) 樋口龍雄, 川又政征, 『ディジタル信号処理-MATLAB 対応-』, 昭晃堂 (2000) FFT(高速フーリエ変換):実装のための数式と実装例 (https://qiita.com/k-kotera/items/4fe20c48be56d2eb4ba4) DFTの動作 N個のデータが以下のようにベクトルの形で与えられるとする。 \boldsymbol{f} = \begin{pmatrix} f_0 \\ \vdots \\ f_{N-1} \\ \end{pmatrix} この時、$\boldsymbol{f}$は直交基底ベクトルの組{$\boldsymbol{e}_0,\boldsymbol{e}_1,\boldsymbol{e}_2,...$}と、その係数の組{$f_0,f_1,f_2,...$}で以下のように表せる。 \boldsymbol{f} = f_0\boldsymbol{e}_0+f_1\boldsymbol{e}_1+...+f_{N-1}\boldsymbol{e}_{N-1} ただし、{$\boldsymbol{e}_0,\boldsymbol{e}_1,\boldsymbol{e}_2,...$}の定義は以下の通りである。 \boldsymbol{e}_0= \begin{pmatrix} 1 \\ 0 \\ \vdots \\ 0 \\ \end{pmatrix} ,\boldsymbol{e}_1= \begin{pmatrix} 0 \\ 1 \\ \vdots \\ 0 \\ \end{pmatrix} ,..., \boldsymbol{e}_{N-1}= \begin{pmatrix} 0 \\ 0 \\ \vdots \\ 1 \\ \end{pmatrix} $\boldsymbol{f}$は別の直交基底ベクトルの組{$\boldsymbol{E}_0,\boldsymbol{E}_1,\boldsymbol{E}_2,...$}と、その係数の組{$F_0,F_1,F_2,...$}で以下のようにも表せる。 \boldsymbol{f} = F_0\boldsymbol{E}_0+F_1\boldsymbol{E}_1+...+F_{N-1}\boldsymbol{E}_{N-1} ただし、$\boldsymbol{E}_j$ $(j=0,1,...,N-1)$ の定義は以下の通りである。 \boldsymbol{E}_j = \begin{pmatrix} (e^{i2\pi\frac{0}{N}})^{j} \\ (e^{i2\pi\frac{1}{N}})^{j} \\ \vdots \\ (e^{i2\pi\frac{N-1}{N}})^{j} \\ \end{pmatrix} この時、{$\boldsymbol{E}_0,\boldsymbol{E}_1,..$}は以下のような性質を満たす。 $k\neq j$とし、 \boldsymbol{E}_j \cdot \boldsymbol{E}_k^*= N,　 \boldsymbol{E}_j \cdot \boldsymbol{E}_j^*= 0 よって以下のような式が成り立つ。 F_j= \frac{1}{N} \ \boldsymbol{f} \cdot \boldsymbol{E}_j^* = \frac{1}{N} {}^t\! \ \boldsymbol{E}_j^* \ \boldsymbol{f} ここで$\boldsymbol{E}_j^T$は$\boldsymbol{E}_j$を転置した横ベクトルである。 またこの係数$F_j$(フーリエ係数)を集めて$\boldsymbol{F}$というベクトルは以下のように定義できる。 \boldsymbol{F}= \begin{pmatrix} F_0 \\ F_{1} \\ \vdots \\ F_{N-1} \\ \end{pmatrix} = \frac{1}{N} \ \begin{pmatrix} {}^t\! \ \boldsymbol{E}_0^* \\ {}^t\! \ \boldsymbol{E}_1^* \\ \vdots \\ {}^t\! \ \boldsymbol{E}_{N-1}^* \\ \end{pmatrix} \ \boldsymbol{f} フーリエ係数のベクトルは元のデータベクトルに上記のような行列をかけることで求められる。 DFTの行列表現 前述した計算式を簡単に書くために$W_N$という$N \times N$行列を以下のように定義する。 W_N= \begin{pmatrix} {}^t\! \ \boldsymbol{E}_0^* \\ {}^t\! \ \boldsymbol{E}_1^* \\ \vdots \\ {}^t\! \ \boldsymbol{E}_{N-1}^* \\ \end{pmatrix} \ = \frac{1}{N} \ \begin{pmatrix} w_N^{0,0} & w_N^{0,1} & \ldots & w_N^{0,N-1} \\ w_N^{1,0} & w_N^{1,1} & \ldots & w_N^{0,N-1} \\ \vdots & \vdots & \ddots & \vdots \\ w_N^{N-1,0} & w_N^{N-1,1} & \ldots & w_N^{N-1,N-1} \\ \end{pmatrix} この時$W_N$の任意の成分$w_N^{j,k}$は以下のように表せる。 w_N^{j,k}=\exp(-i2\pi \frac{j*k}{N}) ,\quad(j,k=0,1,2,...N-1) この行列を用いてフーリエ係数を表すと以下のようになる。 \boldsymbol{F}= W_N \ \boldsymbol{f} FFTの仕組み 任意のフーリエ係数は元のデータベクトルと、基底ベクトルの複素共役との内積であった。 これを和の形に直して以下の式で表記する。 F_j=\boldsymbol{f} \cdot \boldsymbol{E}_j^* = \frac{1}{N} \big\{ \ f_0 w_N^{0,j}+f_1 w_N^{1,j}+..+f_{N-1} w_N^{N-1,j} \ \big\} =\frac{1}{N} \ \Sigma_{k=0}^{N-1} \ f_k w_N^{k,j} この和を$k$が奇数の場合と偶数の場合の二通の輪の取り方に分ける。 F_j =\frac{1}{N} \ \Sigma_{k=0}^{N-1} \ f_k w_N^{k,j} =\frac{1}{N} \big\{ \Sigma_{k'=0}^{ \frac{N}{2} -1} \ f_{2k'} w_N^{2k',j} +\Sigma_{k'=0}^{ \frac{N}{2} -1} \ f_{2k'+1} w_N^{2k'+1,j} \big\} =\frac{1}{N} \big\{ \Sigma_{k'=0}^{ \frac{N}{2} -1} \ f_{2k'} w_N^{2k',j} +w_N^{1,j} \ \Sigma_{k'=0}^{ \frac{N}{2} -1} \ f_{2k'+1} w_N^{2k',j} \big\} ここで w_{N}^{2k , j} = w_{\frac{N}{2}}^{k , j}= w_{\frac{N}{2}}^{k , j+\frac{N}{2}} より F_{j} =\frac{1}{N} \big\{ \Sigma_{k'=0}^{ \frac{N}{2} -1} \ f_{k'} w_{\frac{N}{2}}^{k',j} +w_{N}^{1,j} \ \Sigma_{k'=0}^{ \frac{N}{2} -1} \ f_{k'+1} w_{\frac{N}{2}}^{k',j} \big\} 上記の式は、元のデータ列から偶数番目の成分のフーリエ変換と奇数番目の成分のフーリエ変換が組み合わさったものである。 ここで、右辺第二項の $w_N^{1 , j} $ 注目すると以下のような性質がある。 w_N^{1 , j+ \frac{N}{2}} = \exp(-i2\pi \frac{j+\frac{N}{2}}{N}) = \exp(-i2\pi \frac{j}{N}+-i\pi) =- w_N^{1 , j} つまり以下のような性質がある。 F_{j'} =\frac{1}{N} \big\{ \Sigma_{k'=0}^{ \frac{N}{2} -1} \ f_{k'} w_{\frac{N}{2}}^{k',j'} +w_{N}^{1,j'} \ \Sigma_{k'=0}^{ \frac{N}{2} -1} \ f_{k'+1} w_{\frac{N}{2}}^{k',j'} \big\} \\ F_{j'+\frac{N}{2}} =\frac{1}{N} \big\{ \Sigma_{k'=0}^{ \frac{N}{2} -1} \ f_{k'} w_{\frac{N}{2}}^{k',j'} -w_{N}^{1,j'} \ \Sigma_{k'=0}^{ \frac{N}{2} -1} \ f_{k'+1} w_{\frac{N}{2}}^{k',j'} \big\} \\ （j'=0,1,..,\frac{N}{2} -1） これは大変だ。DFTではバカ正直に$F$の成分一つ一つについて(つまりN個の成分について)計算していたが、FFTでは前半分を計算してしまえば、後半はもうそれを使いまわせばいい。 もちろん掛け算や足し算の計算を余分にしなければならないが、それを考慮しても計算時間が削減できる。 さらに前述したとおり、二つの項はそれぞれ偶数番目、奇数番目の成分のフーリエ変換が組み合わさったものである。このフーリエ変換もFFTでやってしまえばもっともっと計算時間が削減できる。 実際にJuliaとPythonで実装してみてどのようにこのアルゴリズムが動くのか、またどのくらい時間が削減できるかみてみよう。 Juliaでの実装 まず元データ、実軸、周波数軸の用意から。 using Plots im=complex(0,1) #虚数単位 N=Int(2^6) #データ数は2の冪乗で用意 X=1 #xの最大値 dx=X/N #x軸上のデータ間隔 x=Vector(0:dx:X-dx) #x軸のデータ df=1/(X) #周波数軸上のデータ間隔 freq = Vector(0:df:1/dx-df) #周波数軸の生成 freq[Int(N/2)+1:N] .-= N #周波数の単位円上での区間を 0~2π から -π~π に y=sin.(2*pi*3*x) #波形データを用意 plot(x, y, xlabel ="x", ylabel ="y", fmt=:png,size=(500,300)) #波形をプロット こんな画像が出力されます。 次にDFTを実装。 function DFT(y) l = size(y)[1] F=zeros(Complex,Int(l)) for i in 1: l #どの周波数を計算するか指定 for j in 1: l. #シグマをとるためのfor F[i] += y[j]*exp(-2*pi*im*((i-1)*(j-1)/l)) end end return(F) end プロットしてみると以下のようなグラフになる。 F = DFT(y) #大きさの補正 F ./= (N/2) F[1] /= 2 plot(freq, abs.(DFT(y)), xlabel ="freq", ylabel ="amplitude", xlims =(0,N/2), ylims =(0,1.1), fmt=:png,size=(500,300)) つぎにFFTを実装してみる。 function W_n(n) #FFTの第二項で乗算する数値の列を作る。 X = Vector(0:1:n-1) return exp.(-2*pi*im*(X/n)) end function split_vec(X,n) #入力されたデータを奇数番目の成分と偶数番目の成分に分けるための関数 arr=[] for i in 1:n push!(arr,X[i:n:end]) end return arr end function FFT(y) N = size(y)[1] F=zeros(Complex,Int(N)) y_e , y_o=split_vec(y,2) #入力されたデータを奇数番目の成分と偶数番目の成分に分ける if N < 4 #ある程度まで分割したら、残りはDFTで計算。 A = DFT(y_e) B = W_n(N)[1:Int(N/2)] .* DFT(y_o) else A = FFT(y_e) B = W_n(N)[1:Int(N/2)] .* FFT(y_o) end F[1:Int(N/2)]= A.+B F[Int(N/2)+1:end]= A.-B return F end プロットしてみる。 F = FFT(y) F ./= (N/2) F[1] /= 2 plot(freq, abs.(F), xlabel ="freq", ylabel ="amplitude", xlims =(0,N/2), fmt=:png,size=(500,300)) 全く同じグラフが得られた。 一応差をとると、以下のようになる。 plot(freq, abs.(FFT(y) .- DFT(y)), xlabel ="freq", ylabel ="amplitude", xlims =(0,N/2), ylims =(0,10^(-12)), fmt=:png,size=(500,300)) これは結果が一致してると言いて良いだろう。 Juliaによる計算時間の比較。 さて実行時間を比較する。 以下のようなプログラムで計算時間を取得した。 time_DFT=[] #DFTの計算時間 time_FFT=[] #FFTの計算時間 n = []　　 #データ数 for i in 2:12 N=Int(2^i) push!(n,N) X=1 dx=X/N x=Vector(0:dx:X-dx) y=sin.(2*pi*3*x) push!(time_DFT,@elapsed DFT(y)) push!(time_FFT,@elapsed FFT(y)) end 以下はFFT、DETの計算時間をプロットしたもの データ数が多いほどFFTの方が速いことが確かに確認できる。また比をとると以下のようになる。 Pythonでの実装 まず元データ、実軸、周波数軸の用意から。 import numpy as np import matplotlib.pyplot as plt import time im=complex(0,1) N=int(2**6) #データ数は2の冪乗で用意 X=1 #xの最大値 dx=X/N #x軸上のデータ間隔 x=np.linspace(0,X-dx,N) #x軸のデータ df=1/(X) #周波数軸上のデータ間隔 freq =np.linspace(0,1/dx-df,N) #周波数軸の生成 freq[int(N*0.5):] -= N y=np.sin(2*np.pi*3*x) #波形データを用意 plt.plot(x,y) plt.title("Graph Title") plt.xlabel("X") plt.ylabel("Y") plt.show() こんな画像が出力されます。 次にDFTを実装。 def DFT(y): N = len(y) F=np.zeros(N,dtype="complex") for i in range(N): #どの周波数を計算するか指定 for j in range(N): #シグマをとるためのfor F[i] += y[j]*np.exp(-2*np.pi*im*((i-1)*(j-1)/N)) return(F) プロットしてみると以下のようなグラフになる。 F = DFT(y) #大きさの補正 F /= (N/2) F[0] /= 2 plt.plot(freq,np.abs(F)) plt.title("DFT_plot") plt.xlabel("freq") plt.ylabel("amplitude") plt.xlim(0,N/2) plt.show() つぎにFFTを実装してみる。 def W_n(n): #FFTの第二項で乗算する数値の列を作る。 X = np.linspace(0,N-1,N) return np.exp(-2*np.pi*im*(X/n)) def split_vec(X,n): #入力されたデータを奇数番目の成分と偶数番目の成分に分けるための関数 arr=[] N=len(X) for i in range(n): arr.append(X[i:N:n]) return arr def FFT(y): N = len(y) F = np.zeros(N,dtype="complex") y_e , y_o=split_vec(y,2) #入力されたデータを奇数番目の成分と偶数番目の成分に分ける if N < 4 : #ある程度まで分割したら、残りはDFTで計算。 A = DFT(y_e) B = W_n(N)[0:int(N/2)] * DFT(y_o) else: A = FFT(y_e) B = W_n(N)[0:int(N/2)] * FFT(y_o) F[0 : int(N/2)] = A + B F[int(N/2) : N] = A - B return F プロットしてみる。 F = FFT(y) #大きさの補正 F /= (N/2) F[0] /= 2 plt.plot(freq,np.abs(F)) plt.title("Graph Title") plt.xlabel("X") plt.ylabel("Y") plt.xlim(0,N/2) plt.show() 全く同じグラフが得られた。 これは結果が一致してると言いて良いだろう。 Pythonによる計算時間の比較。 さて実行時間を比較する。 以下のようなプログラムで計算時間を取得した。 time_DFT=[] #DFTの計算時間 time_FFT=[] #FFTの計算時間 n = [] #データ数 for i in range(10): N = int(2**(i+2)) #データ数は2の冪乗で用意 n.append(N) X = 1 #xの最大値 dx = X/N #x軸上のデータ間隔 x = np.linspace(0,X-dx,N) #x軸のデータ y = np.sin(2*np.pi*3*x) #波形データを用意 t0 = time.time() DFT(y) time_DFT.append(time.time()-t0) t0 = time.time() FFT(y) time_FFT.append(time.time()-t0) time_DFT = np.array(time_DFT) time_FFT = np.array(time_FFT) 以下はFFT、DETの計算時間をプロットしたもの データ数が多いほどFFTの方が速いことが確かに確認できる。また比をとると以下のようになる。 結論 FFTは、スペクトルを半分計算して残りにその計算結果を流用することで高速化している。 あと、確かにめっちゃ速い。

はじめに こんにちは。自分は現在SE1年目でJavaを業務で扱っています。 データ系の業務へのキャリアチェンジのために統計検定2級を取得しようと考えたのですが、 統計検定2級の勉強+pythonの勉強ができれば一石二鳥だと思い、統計WEBの練習問題をpythonを使って解くことにしました。 同じように統計学もpythonも初心者だという方の参考になれば幸いです。 また、書いてるコードがおかしかったり、もっと良い書き方や慣習等があれば経験者の方のご意見をお聞きしたいです。 目次 具体的な学習方法 練習問題の解答 参考文献 具体的な学習方法 統計WEB：統計学の時間（https://bellcurve.jp/statistics/course/ ） の練習問題をPythonで解いていきます。 学習のためにあえてライブラリを使わずに解いたり、NumPy等を使ってみたり色々な方法で問題を解いていきたいです。 他にもこういう書き方がある、等ありましたらご指摘いただければ幸いです。 ※学習が進み次第随時更新していきます。 練習問題の解答 Pythonで解けそうだと感じた問題のみ扱うので、順番は飛び飛びになりますがご了承ください。 3. さまざまな代表値 問題本文「引用元：練習問題（3. さまざまな代表値）」 https://bellcurve.jp/statistics/course/146.html 3-2 次の表はあるクラス50人の100点満点の国語のテストの結果をまとめたものである。この結果からクラスの平均点、中央値、モードを求めよ。 回答例 回答 import statistics as stat scores = [64, 56, 51, 77, 45, 73, 64, 44, 69, 66, 63, 53, 83, 72, 58, 68, 66, 74, 62, 60, 44, 48, 47, 78, 54, 59, 48, 57, 58, 71, 61, 51, 41, 64, 57, 42, 54, 45, 61, 52, 57, 54, 60, 52, 54, 70, 73, 72, 61, 73] mean = sum(scores) / len(scores) # mean = stat.mean(scores) scores.sort() off = len(scores) // 2 median = (scores[off] + scores[-(off + 1)]) / 2 # median = stat.median(scores) mode = stat.mode(scores) print(f'''\ 平均：{mean} 中央値：{median} 最頻値：{mode}''') 実行結果 平均：59.72 中央値：59.5 最頻値：54 3-3 次の表は、2011年から2015年までの日本の年間GDPの前年比である。 このデータを用いて、5年間のGDPの平均成長率を求めよ。 回答例 回答 import math GDP_list = [0.987, 1.000, 1.017, 1.015, 1.022] sum = 1 for GDP in GDP_list: sum *= GDP mean = math.pow(sum, 1/5) print(f'平均：{mean}') 実行結果 平均：1.008117163331727 3-4 ある100kmの道を、はじめの20kmは時速50kmで、次の50kmは時速40kmで、最後の30kmは時速60kmで走った。この時の平均時速を求めよ。 回答例 回答 mean = 100 / (20 / 50 + 50 / 40 + 30 / 60) print(f'平均時速：{mean}') 実行結果 平均時速：46.5116279069767 3-5 午前と午後でそれぞれ3000匹ずつのひよこを鑑定するノルマを課された 、一人のひよこ鑑定士について考える。ある日、午前中には1時間当たり1200匹のペースで鑑定し、午後は1時間当たり800匹のペースで鑑定した。 このひよこ鑑定士は、1日で平均すると一時間に何匹のペースでひよこを鑑定したことになるか。 回答例 回答 mean = 2 // (1 / 1200 + 1 / 800) print(f'調和平均：{mean}') 実行結果 調和平均：960.0 4. 箱ひげ図と幹葉表示 問題本文「引用元：練習問題（4. 箱ひげ図と幹葉表示）」 https://bellcurve.jp/statistics/course/5592.html 4-2 次のデータから、四分位範囲（IQR）を求めよ。 「56, 48, 78, 81, 86, 71, 72, 88, 46, 47, 89, 58, 43, 79, 48, 41」 回答例 回答 data_list = [56, 48, 78, 81, 86, 71, 72, 88, 46, 47, 89, 58, 43, 79, 48, 41] data_list.sort() def median_calc(target): off = len(target) // 2 median = target[off] + target[-(off + 1)] return median / 2 low_group = data_list[: len(data_list) // 2] high_group = data_list[-(-len(data_list) // 2 ) :] first_iq = median_calc(low_group) third_iq = median_calc(high_group) IQR = third_iq - first_iq print(f'IQR：{IQR}') 実行結果 IQR：32.5 参考文献 統計WEB：統計学の時間

概要 公式ドキュメントではMySQLでやってねっていうけど、sqliteの方が初心者開発者には優しいのでは？と思い、記事にしていきます。今回は短い記事になると思います。 開発環境 vscode mac djagno=3.2.1 python=3.9.1 きっかけ MySQLの中身を見たいけど、何かインストールしないといけない。でも面倒。。。 いつかは入れるけど、今はそんなに使わないしなぁと思い、色々調べたらvscodeでは、何もインストールせずにsqliteの中身を確認できるとのことだったので、試してみました。 内容 vscodeでコマンドパレットというのはご存知でしょうか？ pythonのインタプリタとかを指定するために使っているとは思いますが、DBにうまく値を入れられているかを確認するためにも使えます！コマンドパレットの開き方は、vscodeの「表示」をクリックして、すぐ下にコマンドパレットがあります。または、シフト+command+pでも開くことができます。開けたら、sqliteと打ってみましょう。そしたら、予測にSQLite:Open Databaseと出てくると思いますので、それを押します。そうしたら、db.sqlite3と出てくるので、それも押します。これで準備完了です。このSQLITE EXPLORERが今回の紹介したいものになります。 migrateしたら、ここに新しいtableができるので、右クリックでshow tableをして確認してみてください。私の場合はこんな感じです。 vscodeを使っている人はかなり簡単にできるので、試してみてください。

GoogleMapsApiはリクエストの回数で料金が上がるため、ユーザー毎にリクエストの回数を制限したかった。 GoogleMapsApiにはライブラリがあり、それを使って、python経由でデータを取得できるらしい。 以下サーバーで使えるAPIの種類。 ・Directions API ・Distance Matrix API ・Elevation API ・Geocoding API ・Geolocation API ・Time Zone API ・Roads API ・Places API ・Maps Static API メリットデメリット ・自動レート制限 ・失敗時の再試行 ・簡単な認証 ・非同期または同期 ではクライアントライブラリの使い方。 termianl git clone https://github.com/googlemaps/google-maps-services-python termianl pip install -U googlemaps djangoであればviews.pyに下記を書く import googlemaps from datetime import datetime gmaps = googlemaps.Client(key='ここにAPIkey') 住所検索 geocode_result = gmaps.geocode('東京都　渋谷区', lnaguage='ja') #レスポンス [{'address_components': [{'long_name': '渋谷区', 'short_name': '渋谷区', 'types': ['locality', 'political']}, {'long_name': '東京都', 'short_name': '東京都', 'types': ['administrative_area_level_1', 'political']}, {'long_name': '日本', 'short_name': 'JP', 'types': ['country', 'political']}], 'formatted_address': '日本、東京都渋谷区', 'geometry': {'bounds': {'northeast': {'lat': 35.6921778, 'lng': 139.7238593}, 'southwest': {'lat': 35.6415462, 'lng': 139.6613727}}, 'location': {'lat': 35.6619707, 'lng': 139.703795}, 'location_type': 'APPROXIMATE', 'viewport': {'northeast': {'lat': 35.6921778, 'lng': 139.7238593}, 'southwest': {'lat': 35.6415462, 'lng': 139.6613727}}}, 'place_id': 'ChIJ0Qgx67KMGGARd2ZbObLZHPE', 'types': ['locality', 'political']}] #緯度経度で検索 reverse_geocode_result = gmaps.reverse_geocode((35.6703827, 139.6939093), language='ja') 今後その他にも試してみよう。 インストールで参考にしたサイト https://google-maps-services-python.readthedocs.io/_/downloads/en/latest/pdf/ https://github.com/googlemaps/google-maps-services-python 使えるメソッド(リクエスト、レスポンス) https://developers.google.com/maps/web-services/client-library?hl=ja

この記事について kaggleに取り組むときに使ったライブラリについてまとめた備忘録。 取り組んだコンペは、CommonLit Readability Prize。 ライブラリ集 wandb 学習モデルの可視化や管理ができるサービス Weights & Biases のクライアントアプリ。学習結果のアップロードができる。 Missingno 欠損値を可視化するライブラリ。dataFrameを渡すと、欠損位置を2次元マトリクスで表示してくれる。 rich ターミナルの出力フォーマットを変更できる（richにできる）ライブラリ。日本語にも対応している。 Textstat テキストの統計情報を計算するライブラリ。読みやすさ、複雑度合いグレードレベル（アメリカでドキュメントにつけられる難易度の指標）を算出できる。

wiki の選択ソートを見たところ 最小最大を同時の確定させるロジックのことを double selection sort と呼ぶらしいことを知った 初めて知った・・・お勉強タイムなのでこれを今日は書いてみることにしました とらあえずは昇順ソートと降順ソートを書いてみまして これらを同時に行う感じですが・・・ 最小　と　最大　が入り組んだ場合の考慮が必要だったのですが もっとスマートに書けるのではないかな。 残課題となりました。 selection_sort_asc #昇順（最小値を先頭に移動） def selection_sort_asc(A): for i in range(len(A)): min = i for j in range(i+1, len(A)): if A[j] < A[min]: min = j A[i], A[min] = A[min], A[i] selection_sort_desc #降順（最大値を先頭に移動） def selection_sort_desc(A): for i in range(len(A)): max = i for j in range(i+1, len(A)): if A[j] > A[max]: max = j A[i], A[max] = A[max], A[i] selection_sort_asc2 #昇順（最大値を最後尾に移動） def selection_sort_asc2(A): for i in range(len(A)): max = len(A)-1-i for j in range(len(A)-1-i): if A[j] > A[max]: max = j A[len(A)-1-i], A[max] = A[max], A[len(A)-1-i] double_selection_sort #昇順（最小値を先頭に、最大値を最後尾移動、両方同時に移動） def double_selection_sort(A): for i in range(len(A)//2): min = i max = len(A)-i-1 for j in range(i,len(A)-i): if A[j] < A[min]: min = j if A[j] > A[max]: max = j A[len(A)-1-i], A[max] = A[max], A[len(A)-1-i] if min !=len(A)-1-i: A[i], A[min] = A[min], A[i] else: A[i], A[max] = A[max], A[i] 実行結果 LL = [7,4,3,17,14,13, 11, 12,15, 1, 2,5] L=LL print("src",L) # src [7, 4, 3, 17, 14, 13, 11, 12, 15, 1, 2, 5] double_selection_sort(L) print("double_selection_sort",L) # double_selection_sort [1, 2, 3, 4, 5, 7, 11, 12, 13, 14, 15, 17] L=LL selection_sort_asc(L) print("selection_sort_asc",L) # selection_sort_asc [1, 2, 3, 4, 5, 7, 11, 12, 13, 14, 15, 17] L=LL selection_sort_desc(L) print("selection_sort_desc",L) # selection_sort_desc [17, 15, 14, 13, 12, 11, 7, 5, 4, 3, 2, 1] L=LL selection_sort_asc2(L) print("selection_sort_asc2",L) # selection_sort_asc2 [1, 2, 3, 4, 5, 7, 11, 12, 13, 14, 15, 17] L=LL double_selection_sort(L) print("double_selection_sort",L) # double_selection_sort [1, 2, 3, 4, 5, 7, 11, 12, 13, 14, 15, 17] exit()

はじめに OpenAI GymとKeras-RLを用いた強化学習をしてみました。 調べているとOpenAI GymとKeras-RLを組み合わせた記事はたくさん見つかるのですが、CartPoleなど準備されている環境を試しているものが多いようです。私は自分で環境を作ってみたかったので、勉強のために簡単なゲーム環境を作ってみることにしました。本記事は私が勉強しながら試してみたことをメモしたものになります。簡単な例ですが、これを応用しくことでもっと難しい環境も作っていけるのではないかと期待しています。 自作ゲーム環境を作る ゲームのルール ゲームと呼ぶのもおこがましいですが、下記のような単純な環境を作ります。 Gはゴール、oは自分の位置です。 自分の位置oは左右に動かすことができます。動く量は左に1と10、右に1と10の4パターンです。 oがGに到達したら勝ちです(通過はダメ)。 100回動かしてGに到達できなかったら負けです。 上の動画は、自分の位置oをランダムに動かしたときのものです。このゲームは負けてますね。 ゲーム環境の実装 ゲーム環境はgym.Envから継承して作ります。 最低限、下記のメソッドをオーバーライドする必要があります。 メソッド 説明 reset ゲーム開始時に呼ばれ、環境を初期化する。 step ゲームのステップごとに呼ばれる。 render ゲームの状態を画面に表示する。画面を使用しない場合はpassとしておけば問題ないです。 MovingEnvクラスの定義を行っていきます。 まず定数部分の定義です。 class MovingEnv(gym.Env): # アクションの数。左に1と10、右に1と10移動できるので4となります。 ACTION_NUM = 4 # ゴールの位置と自分の位置の最小値と最大値の座標を0〜99としています。 MIN_POS = 0 MAX_POS = 99 # 最大ステップ数を100としています。このステップ数を超えるとゲームオーバーになります。 MAX_STEPS = 100 # いろいろなサンプルに倣ってレンダーモードを指定しています。 # 特に指定しなくても問題ありません。 metadata = {'render.modes': ['human']} __init__メソッド インスタンスの初期化を行います。 self.observation_spaceには観測値の最小値と最大値を指定します。MovingEnvでは、観測値をゴールの位置と自分の位置の2つです。座標の最小値と最大値はMIN_POSとMAX_POSで定義済みなので、これらの値をnumpy.ndarrayの配列として渡します。 def __init__(self): super().__init__() # Gym環境にアクション数を指定します。 self.action_space = gym.spaces.Discrete(self.ACTION_NUM) # 観測値の最小値と最大値をGym環境に指定します。 # この環境では、観測値としてゴールの位置と自分の位置を使用します。 # それぞれの最小値、最大値をnumpy.ndarrayで指定しています。 self.observation_space = gym.spaces.Box( low=np.array([self.MIN_POS, self.MIN_POS]), high=np.array([self.MAX_POS, self.MAX_POS]), dtype=np.int16) resetメソッド resetメソッドの実装です。ここではゲームの初期化を行います。 def reset(self): # 自分の位置を乱数から取得します。 self.pos = np.random.randint(self.MIN_POS, self.MAX_POS + 1) # ゴールの位置を乱数から取得します。 self.goal = np.random.randint(self.MIN_POS, self.MAX_POS + 1) # ステップ数を0にリセットします。 self.steps = 0 # 観測値を返します。 return np.array([self.goal, self.pos]) stepメソッド stepメソッドの実装です。このメソッドはステップごとに呼び出されます。 引数としてアクションを取ります。アクションの値により下記のように行動するようにします。 アクション 行動 0 自分の座標を-10移動します。 1 自分の座標を-1移動します。 2 自分の座標を+1移動します。 3 自分の座標を+10移動します。 自分の位置が最小値〜最大値から逸脱したときは範囲内に戻します。 rewardは報酬です。ゴールと自分の位置が同じになったとき100としています。 それ以外のときは-1としています。報酬を-1にすると、ゴールへ到達するまでの回数が多ければ多いほど報酬がマイナスになります。強化学習では報酬を多くするように調整されていくため、より少ない回数でゴールに到達できるようになります。 doneはゲームが終了したかどうかを表します。ゴールしたとき、または最大ステップ数を超えたときにTrueにしてゲームを終了します。 stepメソッドは観測値、報酬、終了判定、その他情報の4つの値を返します。 def step(self, action): # アクションによる動作の定義。 if action == 0: next_pos = self.pos - 10 elif action == 1: next_pos = self.pos - 1 elif action == 2: next_pos = self.pos + 1 elif action == 3: next_pos = self.pos + 10 else: next_pos = self.pos # 自分の位置が逸脱していたら調整する。 if next_pos < self.MIN_POS: next_pos = self.MIN_POS elif next_pos > self.MAX_POS: next_pos = self.MAX_POS self.pos = next_pos # ステップ数をカウントアップする。 self.steps += 1 # 報酬を計算する。ゴールに達したら100、それ以外では-1。 reward = 100 if self.pos == self.goal else -1 # ゲーム終了の判定。 done = True if self.pos == self.goal or self.steps > self.MAX_STEPS else False return np.array([self.goal, self.pos]), reward, done, {} renderメソッド renderメソッドの実装です。ゲームの状況を画面に表示します。 このゲームでは、100個の点(.)の上にゴールの位置(G)と自分の位置(o)を描画しているだけです。毎回同じ行に描画するため、print(f'\r〜', end='')のようにしています。 そのまま動かすと一瞬で終わってしまうので、ステップごとに0.1秒スリープさせて自分の位置の動きが見えるようにしています。 def render(self, mode='human'): a = ['.' for x in range(self.MIN_POS, self.MAX_POS + 1)] a[self.goal] = 'G' a[self.pos] = 'o' print(f'\r{"".join(a)}', end='') time.sleep(0.1) 全体のソースは下記のようになります。 import time import numpy as np import gym class MovingEnv(gym.Env): ACTION_NUM = 4 MIN_POS = 0 MAX_POS = 99 MAX_STEPS = 100 metadata = {'render.modes': ['human']} def __init__(self): super().__init__() self.action_space = gym.spaces.Discrete(self.ACTION_NUM) self.observation_space = gym.spaces.Box( low=np.array([self.MIN_POS, self.MIN_POS]), high=np.array([self.MAX_POS, self.MAX_POS]), dtype=np.int16) def reset(self): self.pos = np.random.randint(self.MIN_POS, self.MAX_POS + 1) self.goal = np.random.randint(self.MIN_POS, self.MAX_POS + 1) self.steps = 0 return np.array([self.goal, self.pos]) def step(self, action): ''' action 0: -10移動 1: -1移動 2: +1移動 3: +10移動 ''' if action == 0: next_pos = self.pos - 10 elif action == 1: next_pos = self.pos - 1 elif action == 2: next_pos = self.pos + 1 elif action == 3: next_pos = self.pos + 10 else: next_pos = self.pos if next_pos < self.MIN_POS: next_pos = self.MIN_POS elif next_pos > self.MAX_POS: next_pos = self.MAX_POS self.pos = next_pos self.steps += 1 reward = 100 if self.pos == self.goal else -1 done = True if self.pos == self.goal or self.steps > self.MAX_STEPS else False return np.array([self.goal, self.pos]), reward, done, {} def render(self, mode='human'): a = ['.' for x in range(self.MIN_POS, self.MAX_POS + 1)] a[self.goal] = 'G' a[self.pos] = 'o' print(f'\r{"".join(a)}', end='') time.sleep(0.1) __init__.pyへの登録 自作の環境をOpenAI Gymで動かすには、環境を登録しておく必要があります。登録は__init__.py内で行います。 私の環境では下記のようなフォルダ構成になっているので、 . ├── envs │ ├── __init__.py ・・・ ゲーム環境の登録を行う │ └── moving_env.py ･･･ ゲーム環境定義 └── moving_test1.py ・・・ ゲームを動かすプログラム __init__.pyは下記のように作りました。 from gym.envs.registration import register register( id='movingenv-v0', entry_point='envs.moving_env:MovingEnv', ) idは環境の名前です。OpenAI Gymを動かすときにはこの名前を指定します。 entry_pointは環境ファイルのパスとクラス名を指定します。私の環境だと./envs/moving_env.pyにファイルにMovingEnvという名前のクラスを定義しているので上記のようになります。 これでゲーム環境の作成は完了です。 ランダムに動かしてみる 上で作成したゲーム環境を動かしてみます。 まずは強化学習する前の状態です。env.action_space.sample()を呼び出すとランダムなアクションが取得できます。このアクションを使ってステップを実行してみます。 このプログラムを実行したものが「自作ゲーム環境を作る」にある動画となります。 import gym import envs.moving_env def main(): # ゲーム環境を作成します env = gym.make('movingenv-v0') # ゲーム環境を初期化します。 observation = env.reset() # 無限ループさせてステップを実行していきます。 i = 0 while True: i += 1 # env.action_space.sample()でランダムなアクションを取得できます。 action = env.action_space.sample() # stepに上記で取得したアクションを指定します。 observation, reward, done, _ = env.step(action) # 画面に表示します env.render() # 終了判定 # rewardが0より大きいときがゴールできたときです。 # それ以外はゲームオーバーです。 if done: if reward > 0: print(f' Goal!! steps={i}') else: print(' Game Over') break if __name__ == '__main__': main() 強化学習してみる Keras-RLを使って強化学習をさせます。 モデルはKerasを使って定義します。単純なゲームなので隠れ層が2つのシンプルなモデルとしてます。 import gym from tensorflow import keras from rl.agents.dqn import DQNAgent from rl.policy import BoltzmannQPolicy from rl.memory import SequentialMemory import matplotlib.pyplot as plt import envs.moving_env def main(): # 環境を作成します。 env = gym.make('movingenv-v0') # 環境からアクション数を取得します。このゲームでは4となります。 nb_actions = env.action_space.n # Kerasを使ってモデルを作成します。 model = keras.models.Sequential([ keras.layers.Flatten(input_shape=(1,) + env.observation_space.shape), keras.layers.Dense(16, activation="relu"), keras.layers.Dense(16, activation="relu"), keras.layers.Dense(nb_actions, activation="linear"), ]) # 経験値を蓄積するためのメモリです。学習を安定させるために使用します。 memory = SequentialMemory(limit=50000, window_length=1) # 行動ポリシーはBoltzmannQPolicyを使用しています。 # EpsGreedyQPolicyと比較して、こちらの方が収束が早かったので採用しています。 policy = BoltzmannQPolicy() # DQNAgentを作成します。 dqn = DQNAgent( model=model, nb_actions=nb_actions, memory=memory, target_model_update=1e-2, policy=policy) # DQNAgentのコンパイル。最適化はAdam,評価関数はMAEを使用します。 dqn.compile(keras.optimizers.Adam(lr=1e-3), metrics=['mae']) # 学習を開始します。100000ステップ実行します。 history = dqn.fit(env, nb_steps=100000, visualize=False, verbose=1) # 学習した重みをファイルに保存します。 dqn.save_weights('moving_test.hdf5') # ゲームごとのステップ数と報酬をグラフ化します。 plt.plot(history.history['nb_episode_steps'], label='nb_episode_steps') plt.plot(history.history['episode_reward'], label='episode_reward') plt.legend() plt.show() if __name__ == '__main__': main() これを実行すると下のグラフが表示されます。 グラフは横軸がエピソード数(ゲーム数)、縦軸は青色がゲーム終了までのステップ数、橙色が報酬です。 2000エピソードくらいからゲームオーバーとならずにゴールできているようです。 また、5000エピソードくらいから安定してますね。 では上記の結果を使用してゲームを実行してみます。 import gym from tensorflow import keras from rl.agents.dqn import DQNAgent from rl.policy import BoltzmannQPolicy from rl.memory import SequentialMemory import matplotlib.pyplot as plt import envs.moving_env def main(): env = gym.make('movingenv-v0') nb_actions = env.action_space.n model = keras.models.Sequential([ keras.layers.Flatten(input_shape=(1,) + env.observation_space.shape), keras.layers.Dense(16, activation="relu"), keras.layers.Dense(16, activation="relu"), keras.layers.Dense(nb_actions, activation="linear"), ]) memory = SequentialMemory(limit=50000, window_length=1) policy = BoltzmannQPolicy() dqn = DQNAgent( model=model, nb_actions=nb_actions, memory=memory, target_model_update=1e-2, policy=policy) dqn.compile(keras.optimizers.Adam(lr=1e-3), metrics=['mae']) # ↑ここまでは強化学習と同じソースです。 # 保存した重みを読み込みます。 dqn.load_weights('moving_test.hdf5') # 5エピソード(ゲーム）実行します。 dqn.test(env, nb_episodes=5) if __name__ == '__main__': main() 下が実行したときの動画です。いい感じで学習できています。 学習を早く終わらせるための改善 作成した環境では5000エピソード学習すると安定した動作となりました。 この単純なゲームで5000エピソードは多すぎると思い、環境を変更してみることにしました。変更したところは観測値の返し方です。このゲームではゴールと自分の位置の2つの観測値を返していました。しかしゴールは最初から最後まで変わることがないのでこれは冗長な感じがします。といって、自分の位置だけ返してもゴールの位置がわからないと学習速度が更に遅くなりそうです。そこで観測値として、ゴールの位置と自分の位置の差を返すようにしました。これだと1つの値で自分の位置からみたゴールの方向(左右)と距離を把握することができます。 新しい環境は下記のようになります。(観測値以外は同じです。) import time import numpy as np import gym class MovingEnv(gym.Env): ACTION_NUM = 4 MIN_POS = 0 MAX_POS = 99 MAX_STEPS = 100 metadata = {'render.modes': ['human']} def __init__(self): super().__init__() self.pos = None self.goal = None self.steps = 0 self.action_space = gym.spaces.Discrete(self.ACTION_NUM) self.observation_space = gym.spaces.Box( low=np.array([self.MIN_POS]), high=np.array([self.MAX_POS]), dtype=np.int16) def reset(self): self.pos = np.random.randint(self.MIN_POS, self.MAX_POS + 1) self.goal = np.random.randint(self.MIN_POS, self.MAX_POS + 1) self.steps = 0 return np.array([self.goal - self.pos]) def step(self, action): ''' action 0: -10移動 1: -1移動 2: +1移動 3: +10移動 ''' if action == 0: next_pos = self.pos - 10 elif action == 1: next_pos = self.pos - 1 elif action == 2: next_pos = self.pos + 1 elif action == 3: next_pos = self.pos + 10 else: next_pos = self.pos if next_pos < self.MIN_POS: next_pos = self.MIN_POS elif next_pos > self.MAX_POS: next_pos = self.MAX_POS self.pos = next_pos self.steps += 1 reward = 100 if self.pos == self.goal else -1 done = True if self.pos == self.goal or self.steps > self.MAX_STEPS else False return np.array([self.goal - self.pos]), reward, done, {} def seed(self, seed=None): return seed def render(self, mode='human'): a = ['.' for x in range(self.MIN_POS, self.MAX_POS + 1)] a[self.goal] = 'G' a[self.pos] = 'o' print(f'\r{"".join(a)}', end='') time.sleep(0.1) この環境を使って強化学習してみます。 import gym from tensorflow import keras from rl.agents.dqn import DQNAgent from rl.policy import BoltzmannQPolicy from rl.memory import SequentialMemory import matplotlib.pyplot as plt import envs.moving_env def main(): env = gym.make('movingenv-v0') nb_actions = env.action_space.n model = keras.models.Sequential([ keras.layers.Flatten(input_shape=(1,) + env.observation_space.shape), keras.layers.Dense(16, activation="relu"), keras.layers.Dense(16, activation="relu"), keras.layers.Dense(nb_actions, activation="linear"), ]) memory = SequentialMemory(limit=50000, window_length=1) policy = BoltzmannQPolicy() dqn = DQNAgent( model=model, nb_actions=nb_actions, memory=memory, target_model_update=1e-2, policy=policy) dqn.compile(keras.optimizers.Adam(lr=1e-3), metrics=['mae']) history = dqn.fit(env, nb_steps=10000, visualize=False, verbose=1) dqn.save_weights('moving_test_2.hdf5') plt.plot(history.history['nb_episode_steps'], label='nb_episode_steps') plt.plot(history.history['episode_reward'], label='episode_reward') plt.legend() plt.show() if __name__ == '__main__': main() これを実行すると下のグラフが表示されます。 700エピソードくらいで安定していることがわかります。 観測値をそのまま使用しても強化学習はできますが、環境を工夫することで短時間で学習できるようになることがわかります。 はまったところ 話が逸れますが、DQNAgentを強化学習を実行したところ下記のようなエラーが発生しました。 TypeError: Keras symbolic inputs/outputs do not implement `__len__`. You may be trying to pass Keras symbolic inputs/outputs to a TF API that does not register dispatching, preventing Keras from automatically converting the API call to a lambda layer in the Functional Model. This error will also get raised if you try asserting a symbolic input/output directly. 原因はインストールされているkeras-rlのバージョンが古かったことでした。少し前にTensorFlowのバージョンを2にアップデートしたことで前にインストールしていたkeras-rlとマッチしなくなったようです。keras-rlをアンインストールしてkeras-rl2をインストールすることで解決できました。 解決に時間がかかったのでここにメモしておきます。 $ pip uninstall keras-rl $ pip install keras-rl2 おわりに これでOpenAI GymとKeras-RLを使った基礎を身につけることができました。 強化学習はいろいろなところで使用されています。囲碁のチャンピオンに勝ったAophaGoもそうですし、最近よく耳にする自動運転もその一つです。個人レベルでそこまで壮大なものを作ろうとは思いませんが(作れない)、身近なところでいろいろ応用できそうだと思ってます。

例題5.：対数グラフ 対数グラフの設定方法を学ぶ． 例題演習matplotlibのまとめ https://qiita.com/s4s/items/026aa83b881b640b8b23 例(1)：サンプルコード Y軸が対数となる，片対数グラフを例としたサンプルコードを示す． plot05_example1.py import numpy as np # numpy import matplotlib.pyplot as plt # matplotlib plt import matplotlib.ticker as ticker # matplotlib ticker def main(): '''main program''' # setting f_out = 'graph05_example1.png' # output file (figure) # data x = np.arange(0, 10.1, 0.1) y1 = np.exp(x) #y2 = 1e5 * np.exp(-x) y2 = x data = {'p1': [x, y1, r'$y=\exp(x)$'], # packing 'p2': [x, y2, r'$y=x$']} # plot graph run_plot(f_out, data) def run_plot(f_out, data): '''plot graph [input] f_out: output file name data: data for plotting (x, y, label) ''' # (0) setting of matplotlib # graph gx = 14 # graph size of x [cm] gy = 6 # graph size of y [cm] dpi = 200 # graph DPI (100～600程度) # font f_family = 'IPApGothic' # font (sans-serif, serif, IPApGothic) font_ax = {'size': 10, 'color': 'k'} # 軸ラベルのフォント fs_le = 9 # 凡例のフォントサイズ fs_ma = 9 # 主目盛りのフォントサイズ # 軸ラベル x_label = 'X' # x軸ラベル y_label = 'Y' # y軸ラベル # 軸範囲 x_s = 0 # x軸の最小値 x_e = 10 # x軸の最大値 y_s = 0.1 # y軸の最小値 y_e = 1e5 # y軸の最大値 # 目盛り x_ma = 1 # x軸の主目盛り間隔 x_mi = 0.2 # x軸の副目盛り間隔 y_ma = 0.5 # y軸の主目盛り間隔 y_mi = 0.1 # y軸の副目盛り間隔 # (A) Figure # (A.1) Font or default parameter # https://matplotlib.org/stable/api/matplotlib_configuration_api.html plt.rcParams['font.family'] = f_family # (A.2) Figureの作成 # https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.figure.html # [Returns] # Figureオブジェクトが返される。 # [Parameter] # figsize=(6.4,4.8): グラフサイズ(x,y), 1[cm]=1/2.54[in] # dpi=100: 出力DPI tl = False # Figureのスペース自動調整 fc = 'w' # Figureの背景色 ec = None # 枠線の色 lw = None # 枠線の幅 fig = plt.figure(figsize=(gx/2.54, gy/2.54), dpi=dpi, tight_layout=tl, facecolor=fc, edgecolor=ec, linewidth=lw) # (A.3) グラフ間隔の調整(pltまたはfig) # bottom等はfigure()でtigit_layout=Trueにすれば不要 # https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.subplots_adjust.html bottom = 0.18 # グラフ下側の位置(0～1) top = 0.95 # グラフ上側の位置(0～1) left = 0.12 # グラフ左側の位置(0～1) right = 0.97 # グラフ右側の位置(0～1) hspace = 0.2 # グラフ(Axes)間の上下余白 wspace = 0.2 # グラフ(Axes)間の左右余白 fig.subplots_adjust(bottom=bottom, top=top, left=left, right=right, hspace=hspace, wspace=wspace) # (B) Axes # (B.1) Axesの追加 # 複数のグラフをタイル状に作成することができる。 # 複数グラフの配置を細かく設定する場合はadd_axes()やadd_gridspec()。 # https://matplotlib.org/stable/api/figure_api.html # [Returns] # Axesオブジェクトが返される。 # [Parameter] # *args=(1,1,1): 行番号,列番号,インデックス # 1つのFigureに複数のAxesを作る場合に指定する。 # add_subplot(2,1,1): 2行のAxesを作り、その1番目(上側) # add_subplot(211): 2,1,1と同じだが、簡略記法 # projection='rectilinear': グラフの投影法(polarなど) # sharex: x軸を共有する際に、共有元のAxisを指定する。 # sharey: sharexのy軸版 # label: Axesに対する凡例名（通常は使わない） n_row = 1 # 全グラフの行数 n_col = 1 # 全グラフの列数 n_ind = 1 # グラフ番号 fc = 'w' # Axesの背景色（通常は白） ax = fig.add_subplot(n_row, n_col, n_ind, fc=fc) # (B.2) 軸ラベル(axis label)の設定 # https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.set_xlabel.html # [Parameter] # xlabel: 軸ラベル # loc='center': 軸ラベル位置(x: left,center,right; y:bottom,center,top) # labelpad: 軸ラベルと軸間の余白 # fonddict: 軸ラベルのText設定 ax.set_xlabel(x_label, fontdict=font_ax) # x軸ラベルの設定 ax.set_ylabel(y_label, fontdict=font_ax) # y軸ラベルの設定 # (B.3) 軸の種類の設定 # https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.set_xscale.html xscale = 'linear' # x軸の種類(linear, log, symlog, logit) yscale = 'log' # y軸の種類(linear, log, symlog, logit) ax.set_xscale(xscale) # x軸の種類 ax.set_yscale(yscale) # y軸の種類 # (B.4) 軸の範囲の設定 # 自動で設定する場合は、auto=Trueのみにする # https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.set_xlim.html ax.set_xlim(x_s, x_e, auto=False) # x軸の範囲 ax.set_ylim(y_s, y_e, auto=False) # y軸の範囲 # (B.5) 軸のアスペクト比 # x軸とy軸の比率を強制的に設定する場合に行う。 # https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.set_aspect.html # https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.set_box_aspect.html asp_a = 'auto' # auto:自動, [数値]:x,yの比率 asp_b = None # None:自動, [数値]:x,yのグラフ長さの比率 ax.set_aspect(asp_a) # X軸とY軸の比率を設定する ax.set_box_aspect(asp_b) # X軸とY軸のグラフ長さの比率を設定する # (B.6) 目盛り関係の設定 # (B.6.1)主目盛りの位置(Locator) # https://matplotlib.org/stable/api/ticker_api.html mal_x = ticker.MultipleLocator(x_ma) # 等間隔目盛り # mal_x = ticker.IndexLocator(x_ma, x_off) # 等間隔目盛り(+offset) # mal_x = ticker.LogLocator(base=10) # 対数目盛り # mal_x = ticker.AutoLocator() # 自動目盛り # mal_x = ticker.NullLocator() # 目盛りなし # mal_x = ticker.LinearLocator(nx_ma) # 個数指定目盛り # mal_x = ticker.FixedLocator([0, 1, 3]) # 位置指定目盛り # mal_y = ticker.MultipleLocator(y_ma) # 等間隔目盛り # mal_y = ticker.IndexLocator(y_ma, y_off) # 等間隔目盛り(+offset) mal_y = ticker.LogLocator(base=10) # 対数目盛り # mal_y = ticker.AutoLocator() # 自動目盛り # mal_y = ticker.NullLocator() # 目盛りなし # mal_y = ticker.LinearLocator(ny_ma) # 個数指定目盛り # mal_y = ticker.FixedLocator([-1, 0, 1]) # 位置指定目盛り ax.xaxis.set_major_locator(mal_x) # x軸の主目盛り間隔の設定 ax.yaxis.set_major_locator(mal_y) # y軸の主目盛り間隔の設定 # (B.6.2)主目盛りの表記(Formatter) # https://matplotlib.org/stable/api/ticker_api.html maf_x = ticker.ScalarFormatter() # 数値 # maf_x = ticker.NullFormatter() # 目盛り表記なし # maf_x = ticker.FixedFormatter(['A','B','C','D','E','F']) # 指定表記 # maf_x = ticker.StrMethodFormatter('{x:.1f}m') # format記法 # maf_x = ticker.LogFormatterMathtext(base=10) # log記法(10^x) # maf_y = ticker.ScalarFormatter() # 数値 # maf_y = ticker.NullFormatter() # 目盛り表記なし # maf_y = ticker.FixedFormatter(['A','B','C','D','E','F']) # 指定表記 # maf_y = ticker.StrMethodFormatter('{x:.1f}m') # format記法 maf_y = ticker.LogFormatterMathtext(base=10) # log記法(10^x) ax.xaxis.set_major_formatter(maf_x) # x軸の主目盛り表記の設定 ax.yaxis.set_major_formatter(maf_y) # y軸の主目盛り表記の設定 # (B.6.3) 副目盛りの位置(Locator) # https://matplotlib.org/stable/api/ticker_api.html mil_x = ticker.MultipleLocator(x_mi) # 等間隔目盛り # mil_x = ticker.IndexLocator(x_mi, x_off) # 等間隔目盛り(+offset) # mil_x = ticker.LogLocator(base=10, subs=np.arange(2, 10)*0.1) # 対数目盛り # mil_x = ticker.AutoLocator() # 自動目盛り # mil_x = ticker.NullLocator() # 目盛りなし # mil_x = ticker.LinearLocator(nx_mi) # 個数指定目盛り # mil_y = ticker.MultipleLocator(y_mi) # 等間隔目盛り # mil_y = ticker.IndexLocator(y_mi, y_off) # 等間隔目盛り(+offset) mil_y = ticker.LogLocator(base=10, subs=np.arange(2, 10)*0.1) # 対数目盛り # mil_y = ticker.AutoLocator() # 自動目盛り # mil_y = ticker.NullLocator() # 目盛りなし # mil_y = ticker.LinearLocator(ny_mi) # 個数指定目盛り ax.xaxis.set_minor_locator(mil_x) # x軸の副目盛り間隔の設定 ax.yaxis.set_minor_locator(mil_y) # y軸の副目盛り間隔の設定 # (B.6.4) 副目盛りの表記(Formatter) # https://matplotlib.org/stable/api/ticker_api.html mif_x = ticker.NullFormatter() # 目盛り表記なし mif_y = ticker.NullFormatter() # 目盛り表記なし ax.xaxis.set_minor_formatter(mif_x) # x軸の副目盛り表記の設定 ax.yaxis.set_minor_formatter(mif_y) # y軸の副目盛り表記の設定 # (B.6.5) 目盛り関連の細部の設定 # https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.tick_params.html # [Parameter] # axis: 対象軸(x, y, both) # which: major, minor, both # labelsize: フォントサイズ # top, right, bottom, left: Tickを表記する場所(Trueで表記) # labeltop: Trueで上側にも目盛りを付ける(labelrightも同様) # labelrotation: 目盛りの表記角度 # labelcolor: 目盛り表記の色 # pad: 軸と目盛りの間隔 # color: 目盛り線の色 ax.tick_params(axis='both', which='major', labelsize=fs_ma, top=True, right=True) ax.tick_params(axis='both', which='minor', top=True, right=True) # (B.7) グリッドの設定 # https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.grid.html gb = True # 出力の有無(Falseかつ、lsなどの指定なしで、出力なし) which = 'major' # 目盛りの指定(major, minor, both) axis = 'y' # X,Y軸の指定(x, y, both) color = 'gray' # 線色 ls = '-' # 線種('-', '--', '-.', ':') lw = 0.5 # 線幅 # ax.grid(b=gb, which=which, axis=axis) # Falseの場合 ax.grid(b=gb, which=which, axis=axis, c=color, ls=ls, lw=lw) # (B.8) 直線 # axhline(), axvline(): 横線, 縦線（端から端まで） # axhspan(), asvspan(): 横線, 縦線（始点と終点指定） # axline(): 直線 # https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.axhline.html # [Parameter] # py = 0 # 横線の位置(y座標) # color = 'k' # 線色 # ls = '-' # 線種('-', '--', '-.', ':') # lw = 0.5 # 線幅 # ax.axhline(y=py, c=color, ls=ls, lw=lw) # (C) Plot # (C.1) plot(): 折れ線グラフ # 散布図はscatter() # https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.plot.html # [Returns] # Line2Dオブジェクト（配列）が返される。 # [Parameter] # *args: xデータ, yデータ, [fmt] # x,yデータはx軸とy軸のデータ（リストなど） # fmtは線や点の簡易的記法による指定 # *kwargs: Line2Dの設定値（多数あり） d_x = data['p1'][0] # x data（適宜変更） d_y = data['p1'][1] # y data（適宜変更） label = data['p1'][2] # 凡例（なしはNone） ls = '-' # 線種('-', '--', '-.', ':') lw = 0.7 # 線幅 color = 'b' # 線色(k, r, g, b, etc.) marker = 'None' # marker(o, s, v, ^, D, +, x, etc.; None=なし) mec = 'k' # markerの線色(k, r, g, b, etc.) mew = 0.5 # markerの線幅 mfc = 'w' # markerの色(none, k, r, g, b, w, etc.) ms = 2 # markerサイズ markevery = None # markerの出力間隔(None=全部, 2は2個毎) p2 = ax.plot(d_x, d_y, label=label, ls=ls, lw=lw, c=color, marker=marker, mec=mec, mew=mew, mfc=mfc, ms=ms, markevery=markevery) # (C.1) plot(): 折れ線グラフ d_x = data['p2'][0] # x data（適宜変更） d_y = data['p2'][1] # y data（適宜変更） label = data['p2'][2] # 凡例（なしはNone） ls = '-' # 線種('-', '--', '-.', ':') lw = 0.7 # 線幅 color = 'r' # 線色(k, r, g, b, etc.) marker = 'None' # marker(o, s, v, ^, D, +, x, etc.; None=なし) mec = 'k' # markerの線色(k, r, g, b, etc.) mew = 0.5 # markerの線幅 mfc = 'w' # markerの色(none, k, r, g, b, w, etc.) ms = 2 # markerサイズ markevery = None # markerの出力間隔(None=全部, 2は2個毎) p3 = ax.plot(d_x, d_y, label=label, ls=ls, lw=lw, c=color, marker=marker, mec=mec, mew=mew, mfc=mfc, ms=ms, markevery=markevery) # (D) Others: 凡例、テキストなど # (D.1) 凡例 # (D.1.1) 凡例データの一覧を作成 handles, labels = ax.get_legend_handles_labels() # (D.1.2) 凡例の作成 # https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.legend.html anc = (0, 1) # 凡例の配置場所(0～1の相対位置) loc = 'upper left' # 配置箇所:upper,center,lower; left,center,right; best ti = '' # 凡例のタイトル(なし=None or '') ti_fs = 9 # 凡例のタイトルのフォントサイズ ncol = 1 # 凡例の列数 fc = 'w' # 凡例の背景色（なし=none） ec = 'k' # 枠線の色（なし=none） fa = 1.0 # 透明度(0で透明, 1で透明なし) bp = 0.5 # 凡例の枠の余白(default=0.4) ms = 1.0 # マーカーの倍率(default=1.0) sha = False # True:影付き、False:影なし fb = True # 凡例の角：True:丸、False:四角 bap = 0.8 # 凡例の外側の余白(default=0.5) ls = 0.5 # 凡例間の上下余白 cs = 2.0 # 凡例間の左右余白 legend = ax.legend(handles, labels, title=ti, title_fontsize=ti_fs, loc=loc, bbox_to_anchor=anc, ncol=ncol, facecolor=fc, edgecolor=ec, framealpha=fa, fontsize=fs_le, borderpad=bp, shadow=sha, markerscale=ms, fancybox=fb, labelspacing=ls, columnspacing=cs, borderaxespad=bap) # (D.1.3) 凡例の細かい設定 legend.get_frame().set_linestyle('-') # 枠線の種類 legend.get_frame().set_linewidth(0.5) # 枠線の太さ # (D.2) テキスト # https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.text.html # txt = '(a) Test Plot' # テキスト文字列 # px = 0.97 # x位置(transformで基準座標の変更) # py = 0.05 # y位置(transformで基準座標の変更) # tf = ax.transAxes # 基準座標(ax.transAxes, fig.transFigure) # ha = 'right' # 左右位置(left, center, right) # va = 'bottom' # 上下位置(bottom, center, top, baseline) # rot = 0 # 回転角度(左回転) # color = 'k' # テキストの色 # font_tx = {'size': 10} # テキストのフォント # ax.text(x=px, y=py, s=txt, font=font_tx, ha=ha, va=va, color=color, # rotation=rot, transform=tf) # (E) save figure # https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.savefig.html # [Parameters] # fname: 出力ファイル名 # dpi: 出力DPI（Figure設定時のdpiが優先） print(f'write: {f_out}') plt.savefig(fname=f_out) # save figure # (F) close plt.close() if __name__ == '__main__': main() グラフは次のようになる． 解説 サンプルコードを基に，グラフ表記の変更方法などを説明する． 対数グラフへの変更方法 (1) 対数にする軸の範囲にマイナスが含まれないようにする． # 軸範囲 x_s = 0 # x軸の最小値 x_e = 10 # x軸の最大値 y_s = 0.1 # y軸の最小値 y_e = 1e5 # y軸の最大値 (2) 軸の種類をlinear（線形）からlog（対数）に変更する(B.3)． (B.3)で，軸の種類を変更する．対数はlogとなる． なお，特殊な対数グラフとしてsymlogとlogitがあるけれども，ほとんど使われないので，説明は割愛する． # (B.3) 軸の種類の設定 # https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.set_xscale.html xscale = 'linear' # x軸の種類(linear, log, symlog, logit) yscale = 'log' # y軸の種類(linear, log, symlog, logit) ax.set_xscale(xscale) # x軸の種類 ax.set_yscale(yscale) # y軸の種類 ひとまず，これで対数グラフにはなる．通常は，上記に加えて目盛りの変更が必要となる． (3) 目盛りの設定変更(B.6) サンプルコードのY軸の目盛り設定だけ抜き出して詳しくみてみる． まず，主目盛りの設定は以下の通り． # (B.6) 目盛り関係の設定 # (B.6.1)主目盛りの位置(Locator) mal_y = ticker.LogLocator(base=10) # 対数目盛り ax.yaxis.set_major_locator(mal_y) # y軸の主目盛り間隔の設定 # (B.6.2)主目盛りの表記(Formatter) maf_y = ticker.LogFormatterMathtext(base=10) # log記法(10^x) ax.yaxis.set_major_formatter(maf_y) # y軸の主目盛り表記の設定 目盛り場所を設定するLocatorには，対数用のLocatorであるLogLocator()を用いている．主目盛りの設定パラメータbaseは，対数の底を指定する． 目盛り表記を設定するFormatterには，LogFormatterMathext()を用いている．これは，上の図のように，$10^5$という表記になる．他の表記については，後ほど例を示す． 次に，副目盛りの設定は以下の通り． # (B.6.3) 副目盛りの位置(Locator) mil_y = ticker.LogLocator(base=10, subs=np.arange(2, 10)*0.1) # 対数目盛り ax.yaxis.set_minor_locator(mil_y) # y軸の副目盛り間隔の設定 # (B.6.4) 副目盛りの表記(Formatter) mif_y = ticker.NullFormatter() # 目盛り表記なし ax.yaxis.set_minor_formatter(mif_y) # y軸の副目盛り表記の設定 副目盛りの位置は主目盛りと同じく，LogLocator()だが，パラメータにsubsがある．例の場合は，(0.2, 0.3, ..., 0.9)の位置を指定しており，0.1～1.0の間に9本の副目盛りが作成されることとなる． なお，subs='auto'としておけば，自動で設定されるので，良く分からなければautoにしておけばよい． 副目盛りの表記は，例では表記なし（NullFormatter()）としている． 例(2)：サブグリッド線 対数グラフでは，サブグリッド（副グリッド）の線も示すことがある． その場合は，(B.7)で，サブグリッドの設定を加えれば良い． # (B.7) グリッドの設定 # https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.grid.html gb = True # 出力の有無(Falseかつ、lsなどの指定なしで、出力なし) which = 'major' # 目盛りの指定(major, minor, both) axis = 'y' # X,Y軸の指定(x, y, both) color = 'gray' # 線色 ls = '-' # 線種('-', '--', '-.', ':') lw = 0.5 # 線幅 # ax.grid(b=gb, which=which, axis=axis) # Falseの場合 ax.grid(b=gb, which=which, axis=axis, c=color, ls=ls, lw=lw) gb = True # 出力の有無(Falseかつ、lsなどの指定なしで、出力なし) which = 'minor' # 目盛りの指定(major, minor, both) axis = 'y' # X,Y軸の指定(x, y, both) color = 'gray' # 線色 ls = '-' # 線種('-', '--', '-.', ':') lw = 0.3 # 線幅 ax.grid(b=gb, which=which, axis=axis, c=color, ls=ls, lw=lw) グラフは次のようになり，Y軸にサブグリッドの線が加わっている． 例(3)：非常に広い範囲の対数目盛りの表記方法（目盛りの強制表示） Y軸の範囲の上限値を，$10^5$から$10^7$に変化させると，次のように，Y軸の主目盛りの数が間引きされ，副目盛りは表示されなくなる． これはmatplotlibが目盛りを自動調整するためだが，強制的に出力させたい場合もある． その場合は，LogFormatter()のパラメータnumticks（ticksの上限数）を設定する．サンプルコードでは，以下のように変更する． mal_y = ticker.LogLocator(base=10, numticks=9) # 対数目盛り 作成されるグラフでは，主目盛りが9個に増えている． さらに，副目盛りも強制出力するには，以下のように変更すれば良い． なお，numticksは出力する目盛り数の上限値なので，999などと設定しても結果は同じ． mil_y = ticker.LogLocator(base=10, subs=np.arange(2, 10)*0.1, numticks=8*9) # 対数目盛り グラフは以下のようになる． 例(4)：E表記法（コンピューター的指数表記） E表記法とは，「10^5」を「1e+05」と表記する方法で，プログラミング言語やExcelではおなじみの表記方法である．E表記法は，Formatterに，LogFormatter()を使えば良い（あまり一般的でないので，サンプルコードからは除外している）． サンプルコードでは，以下の部分を変更する． # (A.3) グラフ間隔の調整(pltまたはfig) left = 0.14 # グラフ左側の位置(0～1) # (B.6.2)主目盛りの表記(Formatter) maf_y = ticker.LogFormatter(base=10) # log記法(e表記) グラフは以下のようになる．1～$10^4$までは整数で表記され，1以下と$10^4$以上はE表記となる． 例(5)：副目盛りの位置変更 あまり一般的ではないが，常用対数の副目盛りを0.2と0.5の2つとする場合がある（倍半分が分かりやすい）． 下記のように，副目盛りのLogLocatorのパラメータsubsを変更する． # (B.6.3) 副目盛りの位置(Locator) mil_y = ticker.LogLocator(base=10, subs=[0.2, 0.5]) # 対数目盛り グラフは次のようになる． 演習 サンプルコードを基に，下記の演習問題に解答せよ． 解答はコメント欄を参照． 演習(1)：両対数グラフ 次の図のように，X軸も対数グラフにせよ． ※X軸の範囲を0.1～10，目盛りを対数用に変更． 演習(2)：グリッド線の変更 演習(1)の図から，グリッド線を追加し，以下の図のように変更せよ． ※グリッド線はX軸とY軸両方で作成，メイングリッド線は線幅0.5の実線，サブグリッド線は線幅0.3の破線．

はじめに 私はプログラミング歴１年の初心者です。 実務でWebサイトのコーディングを１年間行ってきました。 そろそろシステム開発もできるようになりたいということで LaravelやReactをこれから勉強していこうと思っております。 今回の目的 Pythonでできることを学びます！ Pythonにできること Pythonのアルゴリズム 線形探索を行える 2分探索を行える ハッシュ表探索を行える バブルソートができる 挿入ソートができる 選択ソートができる シェルソートができる ヒープソートができる 再帰関数が書ける クイックソートができる マージソートができる Pythonのオブジェクト指向 クラスを作成できる インスタンスを作成できる インスタンスメソッドを定義できる クラスを継承したサブクラスを定義できる ポリモーフィズムで抽象クラスと抽象メソッドを作成できる(抽象メソッドはサブクラスによって内容が変わる) カプセル化ができる パブリック変数とプライベート変数を定義できる ゲッターとセッターをセットすることでプライベート変数にアクセスできる 集約機能を使うことができる Pythonのクラス クラスには属性とメソッドを定義できる インスタンスを生成できる init()メソッドで初期値を実装する クラスを継承することができる メソッドはオーバーライドできる サブクラスにメソッドを追加することができる super()でスーパークラスを呼び出すことができる クラス変数はクラス全体で共有する値を扱う数 クラスメソッドはクラス全体で共有した値を扱ったメソッド 特殊メソッドでクラスに独自の処理を追加することができる 列挙型クラスで複数の定数をひとまとめにできる

Selenium(Python)で要素をクリックする直前に、毎回要素までスクロールさせる処理を行うためのコードです。 このGoogle検索でスクロール処理させる必要は本来ないはずですが、例として書いています。 test.py from selenium import webdriver from webdriver_manager.chrome import ChromeDriverManager from selenium.webdriver.support.events import EventFiringWebDriver, AbstractEventListener #要素をクリックする直前の処理を定義 class CustomListener(AbstractEventListener): def before_click(self, element, driver): # 要素までスクロールさせる driver.execute_script('arguments[0].scrollIntoView({behavior: "smooth", block: "center"});', element) # Webdriver ManagerでChromeDriverを取得 base_driver = webdriver.Chrome(ChromeDriverManager().install()) driver = EventFiringWebDriver(base_driver, CustomListener()) # Google検索で「Selenium」の検索結果を1〜10ページ目まで表示 driver.get('https://www.google.com/search?q=selenium') for i in range(2, 11): driver.find_element_by_xpath("//a[@aria-label='Page " + str(i) +"']").click() driver.quit() 参考 あまり知られてなさそうなメソッド element.scrollIntoView() Element.scrollIntoView() | MDN 【Python】before_click/after_click・・・要素がクリックされる直前/直後の処理を実施する | Selenium クイックリファレンス

製造業出身のデータサイエンティストがお送りする記事 今回は時系列解析手法のProphetをOptunaでハイパラチューニングをしてみました はじめに 過去に時系列解析に関しては何個か整理しておりますので、興味ある方は参照して頂けますと幸いです。 時系列解析 状態空間モデル 時系列予測をGBDTで実装してみた 時系列解析手法のSARIMAモデルを試してみた Prophetとは Prophet は Facebook が公開しているオープンソースの時系列解析ライブラリです。 Prophetのメリットは大きく下記があります。 1.モデル式設計の柔軟性が高い 2.サンプル間隔の柔軟性が高い 3.高速に計算可能 4.パラメータが理解しやすい 詳細はProphetの公式ドキュメントを参照してください。 Prophet実装 今回も使用したデータは、月ごとの飛行機の乗客数データです。 # ライブラリーのインポート import pandas as pd import numpy as np import seaborn as sns import optuna from matplotlib import pylab as plt %matplotlib inline from fbprophet import Prophet from fbprophet.plot import add_changepoints_to_plot import warnings warnings.filterwarnings("ignore") # https://stat.ethz.ch/R-manual/R-devel/library/datasets/html/AirPassengers.html df = pd.read_csv('../data/AirPassengers.csv') # float型に変換 df['#Passengers'] = df['#Passengers'].astype('float64') df = df.rename(columns={'#Passengers': 'Passengers'}) # datetime型に変換にする df.Month = pd.to_datetime(df.Month) # データの中身を確認 df.head() # データの可視化 fig, ax = plt.subplots() a = sns.lineplot(x="Month", y="Passengers", data=df) plt.show() Prophetでモデルを構築します。 def objective_variable(train, valid): def objective(trial): params = { 'changepoint_range' : trial.suggest_discrete_uniform('changepoint_range', 0.8, 0.95, 0.001), 'n_changepoints' : trial.suggest_int('n_changepoints', 20, 35), 'changepoint_prior_scale' : trial.suggest_discrete_uniform('changepoint_prior_scale',0.001, 0.5, 0.001), 'seasonality_prior_scale' : trial.suggest_discrete_uniform('seasonality_prior_scale',1, 25, 0.1), } # fit_model model = Prophet( changepoint_range = params['changepoint_prior_scale'], n_changepoints=params['n_changepoints'], changepoint_prior_scale=params['changepoint_prior_scale'], seasonality_prior_scale = params['seasonality_prior_scale'], ) model.fit(train) future = model.make_future_dataframe(periods=len(valid)) forecast = model.predict(future) valid_forecast = forecast.tail(len(valid)) val_mape = np.mean(np.abs((valid_forecast.yhat-valid.y)/valid.y))*100 return val_mape return objective def optuna_parameter(train, valid): study = optuna.create_study(sampler=optuna.samplers.RandomSampler(seed=10)) study.optimize(objective_variable(train, valid), timeout=500) optuna_best_params = study.best_params return study df = df.rename(columns={'Month':'ds','Passengers':'y'}) df = df[['ds','y']] df_train = df[df['ds'] < '1956-04-01'] df_valid = df[(df['ds'] >= '1956-04-01')&(df['ds'] < '1957-04-01')] df_test = df[df['ds'] >= '1957-04-01'] study = optuna_parameter(df_train, df_valid) あとは、チューニング結果のハイパーパラメータを使ってモデルの学習と予測を行います。 # fit_model best_model = Prophet( changepoint_range = study.best_params['changepoint_prior_scale'], n_changepoints=study.best_params['n_changepoints'], seasonality_prior_scale = study.best_params['seasonality_prior_scale'], changepoint_prior_scale=study.best_params['changepoint_prior_scale'], ) best_model.fit(df_train) feature_test = best_model.make_future_dataframe(periods=len(df_valid)+len(df_test), freq='M') forecast_test = best_model.predict(feature_test) forecast_test_plot = best_model.plot(forecast_test) fig, ax = plt.subplots() df.y.plot(ax=ax, label='Original', linestyle="dashed") forecast_test.yhat.plot(ax=ax, label='Predict') ax.legend() さいごに 最後まで読んで頂き、ありがとうございました。 今回は時系列解析手法のProphetをOptunaでハイパラチューニングをしてみました。 訂正要望がありましたら、ご連絡頂けますと幸いです。

指定したディレクトリ内の画像ファイルをすべてExcelに貼り付けるプログラム。 参考にしたページ import os import glob import imghdr import openpyxl import cv2 from tkinter import filedialog # 定数設定 INI_IMG_DIR = "D:\Python" # 貼り付ける画像を探す際の初期ディレクトリ SHEET_TITLE = '画像貼り付け' # 画像を貼り付けるエクエルシート名 # 変数 max_height = [] # 各行の画像の高さの最大値を保持 def get_file_names(set_dir_name): """ ディレクトリ内のファイル名取得（ファイル名のみの一覧を取得） """ file_names = os.listdir(set_dir_name) temp_full_file_names = [os.path.join(set_dir_name, file_name) for file_name in file_names if os.path.isfile(os.path.join(set_dir_name, file_name))] # ファイルかどうかを判定 return temp_full_file_names def attach_img(target_full_file_names, set_column_idx, set_dir_name): """ 画像を呼び出して、Excelに貼り付け """ set_row_idx = 1 column_letter = ws.cell(row=set_row_idx, column=set_column_idx).coordinate[0] ws.cell(row=1, column=set_column_idx).value = set_dir_name # 各列の1行目に、貼り付ける画像があるディレクトリ名を入力 max_width = 0 # 画像の幅の最大値を保持するための変数 target_full_file_names.sort() # ファイル名でソート for target_file in target_full_file_names: if imghdr.what(target_file) != None: # 画像ファイルかどうかの判定 img = openpyxl.drawing.image.Image(target_file) print('[' + column_letter + '][' + str(set_row_idx+1) + ']' + target_file + 'を貼り付け') # 画像のサイズを取得して、セルの大きさを合わせる（画像同士が重ならないようにするため） size_img = cv2.imread(target_file) height, width = size_img.shape[:2] if max_width < width: max_width = width if not max_height[set_row_idx-1:set_row_idx]: # 配列「max_height」において、「set_row_idx」番目の要素が存在しなければ、挿入 max_height.insert(set_row_idx-1, height) if max_height[set_row_idx-1] < height: max_height[set_row_idx-1] = height ws.row_dimensions[set_row_idx+1].height = max_height[set_row_idx-1] * 0.75 ws.column_dimensions[column_letter].width = max_width * 0.13 cell_address = ws.cell(row=set_row_idx + 1, column=set_column_idx).coordinate # セルの行列番号から、そのセルの番地を取得 img.anchor = cell_address ws.add_image(img) # シートに画像貼り付け set_row_idx += 1 # ワークブック設定 wb = openpyxl.Workbook() ws = wb.worksheets[0] # 1番目のシートを編集対象にする ws.title = SHEET_TITLE # 1番目のシートに名前を設定 # 貼り付ける画像を置いておくルートディレクトリを指定 ※日本語含むとうまくいかなかった input_img_file = filedialog.askdirectory( title = "貼り付ける画像のフォルダを指定", initialdir=INI_IMG_DIR) # 貼り付ける画像を置いておくルートディレクトリ内のディレクトリ名を再帰的に取得 dirs = glob.glob(os.path.join(input_img_file, '**' + os.sep), recursive=True) column_idx = 1 # 各ディレクトリについて操作 for dir_name in dirs: f_names = get_file_names(dir_name) # ファイル名取得 attach_img(f_names, column_idx, dir_name) # 画像貼り付け設定 column_idx += 1 # 次の列へ・・・ # ファイルへの書き込み save_file_path = filedialog.asksaveasfilename( title = "名前を付けて保存", filetypes = [("xlsx", ".xlsx"), ("xls", ".xls")], # ファイルフィルタ initialdir = "./", # 自分自身のディレクトリ defaultextension = ".xlsx" ) wb.save(save_file_path)

製造業出身のデータサイエンティストがお送りする記事 今回は勾配ブースティング決定木（XGBoost, LightGBM, CatBoost）でOptunaを使ってみました。 はじめに 勾配ブースティング木に関しては、過去に記事を書いておりますのでそちらを参照して頂けますと幸いです。 勾配ブースティング決定木（XGBoost, LightGBM, CatBoost）を実装してみた LightGBMのハイパーパラメータチューニング（Optuna）をしてみた XGBoostの実装 まずは、XGBoostをOptunaでハイパラチューニングを実装してみます。 今回もUCI Machine Learning Repositoryで公開されているボストン住宅の価格データを用いて予測モデルを構築します。 # ライブラリーのインポート import os import pandas as pd import numpy as np import seaborn as sns import matplotlib.pyplot as plt %matplotlib inline # ボストンの住宅価格データ from sklearn.datasets import load_boston # 前処理 from sklearn.preprocessing import StandardScaler from sklearn.model_selection import train_test_split # XGBoost import xgboost as xgb # Optuna import optuna from sklearn.model_selection import cross_val_score # 評価指標 from sklearn.metrics import r2_score from sklearn.metrics import mean_absolute_error from sklearn.metrics import mean_squared_error # データセットの読込み boston = load_boston() # 説明変数の格納 df = pd.DataFrame(boston.data, columns = boston.feature_names) # 目的変数の追加 df['MEDV'] = boston.target # データの中身を確認 df.head() 次にoptunaの事前準備をします。 # ランダムシード値 RANDOM_STATE = 10 # 学習データと評価データの割合 TEST_SIZE = 0.2 # 学習データと評価データを作成 x_train, x_test, y_train, y_test = train_test_split(df.iloc[:, 0:df.shape[1]-1], df.iloc[:, df.shape[1]-1], test_size=TEST_SIZE, random_state=RANDOM_STATE) def objective(trial): eta = trial.suggest_loguniform('eta', 1e-8, 1.0) gamma = trial.suggest_loguniform('gamma', 1e-8, 1.0) max_depth = trial.suggest_int('max_depth', 1, 10) min_child_weight = trial.suggest_loguniform('min_child_weight', 1, 40) max_delta_step = trial.suggest_loguniform('max_delta_step', 1e-8, 1.0) subsample = trial.suggest_uniform('subsample', 0.0, 1.0) reg_lambda = trial.suggest_uniform('reg_lambda', 0.0, 1000.0) reg_alpha = trial.suggest_uniform('reg_alpha', 0.0, 1000.0) model = xgb.XGBRegressor(eta = eta, gamma = gamma, max_depth = max_depth, min_child_weight = min_child_weight, max_delta_step = max_delta_step, subsample = subsample, reg_lambda = reg_lambda, reg_alpha = reg_alpha ) score = cross_val_score(model, x_train, y_train, cv=5, scoring="neg_mean_absolute_error" ) mae = score.mean() return mae 各パラメータについては、公式文章を参照してください。 次にoptunaで最適なハイパーパラメータを見つけます。 # optunaで最適値を見つける # 注：cross_val_scoreの出力は全て高いほど良い study = optuna.create_study(direction='maximize') study.optimize(objective, n_trials=500) あとは、チューニング結果のハイパーパラメータを使ってモデルの学習と予測を行います。 # チューニングしたハイパーパラメーターをフィット optimised_model = xgb.XGBRegressor(eta = study.best_params['eta'], gamma = study.best_params['gamma'], max_depth = study.best_params['max_depth'], min_child_weight = study.best_params['min_child_weight'], max_delta_step = study.best_params['max_delta_step'], subsample = study.best_params['subsample'], reg_lambda = study.best_params['reg_lambda'], reg_alpha = study.best_params['reg_alpha'] ) optimised_model.fit(x_train ,y_train) # XGBoost推論 y_pred = optimised_model.predict(x_test) # 評価 def calculate_scores(true, pred): """全ての評価指標を計算する Parameters ---------- true (np.array) : 実測値 pred (np.array) : 予測値 Returns ------- scores (pd.DataFrame) : 各評価指標を纏めた結果 """ scores = {} scores = pd.DataFrame({'R2': r2_score(true, pred), 'MAE': mean_absolute_error(true, pred), 'MSE': mean_squared_error(true, pred), 'RMSE': np.sqrt(mean_squared_error(true, pred))}, index = ['scores']) return scores scores = calculate_scores(y_test, y_pred) print(scores) 出力結果は下記のようになります。 R2 MAE MSE RMSE scores 0.85572 2.570919 15.08891 3.884445 LightGBMの実装 LightGBMは下記に実装サンプルがありますので、そちらの記事を参考にしてください。 LightGBMのハイパーパラメータチューニング（Optuna）をしてみた CatBoostの実装 最後に、CatBoostを実装します。 # ライブラリーのインポート import os import pandas as pd import numpy as np import seaborn as sns import matplotlib.pyplot as plt %matplotlib inline # ボストンの住宅価格データ from sklearn.datasets import load_boston # 前処理 from sklearn.preprocessing import StandardScaler from sklearn.model_selection import train_test_split # CatBoost import catboost as cb from catboost import CatBoost, Pool # Optuna import optuna from sklearn.model_selection import cross_val_score # 評価指標 from sklearn.metrics import r2_score from sklearn.metrics import mean_absolute_error from sklearn.metrics import mean_squared_error # データセットの読込み boston = load_boston() # 説明変数の格納 df = pd.DataFrame(boston.data, columns = boston.feature_names) # 目的変数の追加 df['MEDV'] = boston.target # データの中身を確認 df.head() 次にoptunaの事前準備をします。 # ランダムシード値 RANDOM_STATE = 10 # 学習データと評価データの割合 TEST_SIZE = 0.2 # 学習データと評価データを作成 x_train, x_test, y_train, y_test = train_test_split(df.iloc[:, 0:df.shape[1]-1], df.iloc[:, df.shape[1]-1], test_size=TEST_SIZE, random_state=RANDOM_STATE) def objective(trial): iterations = trial.suggest_int('iterations', 50, 300) depth = trial.suggest_int('depth', 4, 10) learning_rate = trial.suggest_loguniform('learning_rate', 0.01, 0.3) random_strength = trial.suggest_int('random_strength', 0, 100) bagging_temperature = trial.suggest_loguniform('bagging_temperature', 0.01, 100.00) od_type = trial.suggest_categorical('od_type', ['IncToDec', 'Iter']) od_wait = trial.suggest_int('od_wait', 10, 50) model = cb.CatBoostRegressor(iterations = iterations, depth = depth, learning_rate = learning_rate, random_strength = random_strength, bagging_temperature = bagging_temperature, od_type = od_type, od_wait = od_wait ) score = cross_val_score(model, x_train, y_train, cv=5, scoring="neg_mean_absolute_error" ) mae = score.mean() return mae 次にoptunaで最適なハイパーパラメータを見つけます。 # optunaで最適値を見つける # 注：cross_val_scoreの出力は全て高いほど良い study = optuna.create_study(direction='maximize') study.optimize(objective, n_trials=10) あとは、チューニング結果のハイパーパラメータを使ってモデルの学習と予測を行います。 # チューニングしたハイパーパラメーターをフィット optimised_model = cb.CatBoostRegressor(iterations = study.best_params['iterations'], depth = study.best_params['depth'], learning_rate = study.best_params['learning_rate'], random_strength = study.best_params['random_strength'], bagging_temperature = study.best_params['bagging_temperature'], od_type = study.best_params['od_type'], od_wait = study.best_params['od_wait'] ) optimised_model.fit(x_train ,y_train) # CatBoost推論 y_pred = optimised_model.predict(x_test) # 評価 def calculate_scores(true, pred): """全ての評価指標を計算する Parameters ---------- true (np.array) : 実測値 pred (np.array) : 予測値 Returns ------- scores (pd.DataFrame) : 各評価指標を纏めた結果 """ scores = {} scores = pd.DataFrame({'R2': r2_score(true, pred), 'MAE': mean_absolute_error(true, pred), 'MSE': mean_squared_error(true, pred), 'RMSE': np.sqrt(mean_squared_error(true, pred))}, index = ['scores']) return scores scores = calculate_scores(y_test, y_pred) print(scores) 出力結果は下記のようになります。 R2 MAE MSE RMSE scores 0.865739 2.541304 14.041161 3.747154 さいごに 最後まで読んで頂き、ありがとうございました。 今回は勾配ブースティング決定木の3つのアルゴリズム（XGBoost, LightGBM, CatBoost）でOptunaを使ってみました。 訂正要望がありましたら、ご連絡頂けますと幸いです。

はじめに 私はプログラミング歴１年の初心者です。 実務でWebサイトのコーディングを１年間行ってきました。 そろそろシステム開発もできるようになりたいということで LaravelやReactをこれから勉強していこうと思っております。 今回の目的 Djangoでアプリを実際に作っていく流れを学ぼうと思います。 公式チュートリアルを読んでいきます。 目次 アセット 管理画面のセット 実践 アセット スタイルシート # polls/staticディレクトリを作成する $ mkdir polls/static # polls/static/pollsディレクトリを作成する $ mkdir polls/static/polls polls/static/polls/style.css # スタイルをつける li a { color: green; } polls/templates/polls/index.html # 上部に追加する {% load static %} <link rel="stylesheet" type="text/css" href="{% static 'polls/style.css' %}"> 画像追加 # 画像を追加する background.gif polls/static/polls/style.css # 背景画像にセット body { background: white url("images/background.gif") no-repeat; } 管理画面のセット 管理画面のカスタマイズ polls/admin.py # 詳細画面の順番を入れ替えるカスタマイズ class QuestionAdmin(admin.ModelAdmin): fields = ['pub_date', 'question_text'] admin.site.register(Question, QuestionAdmin) polls/admin.py # 詳細画面のレイアウトを更に変更 class QuestionAdmin(admin.ModelAdmin): fieldsets = [ (None, {'fields': ['question_text']}), ('Date information', {'fields': ['pub_date']}), ] admin.site.register(Question, QuestionAdmin) polls/admin.py # リレーションを張ったオブジェクトの追加 from django.contrib import admin from .models import Choice, Question admin.site.register(Choice) ↓↓↓↓ さらに以下のようにすると ↓↓↓↓ class ChoiceInline(admin.StackedInline): model = Choice extra = 3 class QuestionAdmin(admin.ModelAdmin): fieldsets = [ (None, {'fields': ['question_text']}), ('Date information', {'fields': ['pub_date'], 'classes': ['collapse']}), ] inlines = [ChoiceInline] admin.site.register(Question, QuestionAdmin) polls/admin.py # 管理画面の一覧ページをカスタマイズする class QuestionAdmin(admin.ModelAdmin):内で list_display = ('question_text', 'pub_date', 'was_published_recently') list_filter = ['pub_date'] search_fields = ['question_text'] polls/models.py # 管理画面の一覧ページの表示を改善する from django.contrib import admin class Question(models.Model): # ... @admin.display( boolean=True, ordering='pub_date', description='Published recently?', ) def was_published_recently(self): now = timezone.now() return now - datetime.timedelta(days=1) <= self.pub_date <= now 設定ファイルをカスタムする djangosite/djangosite/settings.py # プロジェクトテンプレートをカスタマイズする TEMPLATES内の 'DIRS': [BASE_DIR / 'templates'], # djangositeディレクトリ直下にtemplatesディレクトリを作成 $ mkdir templates # templatesディレクトリ直下にadminディレクトリを作成 $ mkdir templates/admin # djangoのソースファイルの場所を探して $ python -c "import django; print(django.__path__)" # デフォルトDjango adminテンプレートディレクトリを探す $ cd djangoのソースファイルの場所 $ cd contrib/admin/templates # さらにその中のadmin/base_site.htmlを新しく作ったディレクトリにコピーする $ cp admin/base_site.html ~/my_study/djangosite/templates/admin templates/admin/base_site.html # base_site.htmlを編集する {% extends "admin/base.html" %} {% block title %}{% if subtitle %}{{ subtitle }} | {% endif %}{{ title }} | {{ site_title|default:_('Django site admin') }}{% endblock %} {% block branding %} <h1 id="site-name"><a href="{% url 'admin:index' %}">Polls Administration</a></h1> {% endblock %} {% block nav-global %}{% endblock %}

はじめに XenonPyのチュートリアルのDescriptor calculationについて解説します。XenonPyの概要とチュートリアルの一部（Sample dataとData access）の解説についてはこちらの投稿を参考にして下さい。 ※マテリアルズインフォマティクス関係の内容を他にも投稿していますので、よろしければこちらの一覧から他の投稿も見て頂けますと幸いです。 記述子計算 公式ドキュメントの"Descriptor calculation"に沿って解説します。 XenonPyには記述子計算のインターフェイスが含まれています 4つのタイプに分類される計15個のFeaturizerを計算することができる この投稿では"Composition"と"Structure"タイプの記述子計算について解説いたします。 環境 conda 4.9.2 python 3.7.1 xenonpy 0.4.4 Compositional descriptors（組成に関する記述子） XenonPyでは58個の元素レベルの特徴量（"elements_completed"）から290個の組成に関する特徴量（compositional features）を計算することが可能です。 from xenonpy.descriptor import Compositions cal = Compositions() cal 出力結果 Compositions: |- composition: | |- Counting | |- WeightedAverage | |- WeightedSum | |- WeightedVariance | |- GeometricMean | |- HarmonicMean | |- MaxPooling | |- MinPooling Compositionsには、Countingをはじめとする8つのfeaturizerが含まれていることが分かります。それぞれのfeaturizerの定義はこちらを参考にして下さい。この計算機を使うことで化合物の組成に関する情報を得ることができます。 サンプルデータを例に使い方を以下に示します。 from xenonpy.datatools import preset samples = preset.mp_samples comps = samples['composition'] descriptor = cal.transform(comps) descriptor 807個の化合物それぞれに対して、290個の記述子が得られたことが分かります。 ここで、変換前のサンプルデータを確認してみます。 samples 次に.transformに入力したデータを確認してみます。 samples['composition'] 出力結果 mp-1008807 {'Rb': 1.0, 'Cu': 1.0, 'O': 1.0} mp-1009640 {'Pr': 1.0, 'N': 1.0} mp-1016825 {'Hf': 1.0, 'Mg': 1.0, 'O': 3.0} mp-1017582 {'La': 1.0, 'Pt': 3.0, 'C': 1.0} mp-1021511 {'Cd': 1.0, 'S': 1.0} ... mp-9846 {'Rb': 2.0, 'Ca': 2.0, 'Sb': 2.0} mvc-4727 {'Ca': 4.0, 'Te': 2.0, 'O': 10.0} mvc-6767 {'Zn': 2.0, 'Ni': 1.0, 'W': 1.0, 'O': 6.0} mvc-6973 {'Cu': 6.0, 'W': 12.0, 'O': 42.0} mvc-89 {'Ca': 4.0, 'Nb': 4.0, 'Cu': 2.0, 'O': 16.0} Name: composition, Length: 807, dtype: object transformメソッドはpymatgen.Structureオブジェクトか上記のような組成の辞書形式のデータ（例：{‘H’:2,‘O’:1}）を受け取ることができます。ですので、以下のように自ら指定する任意の化合物に関して記述子を得ることが可能です。例として、H2o、LiF、CO2の記述子を取得してみます。 # 任意の化合物を定義 H2O = {'H': 2.0, 'O': 1.0} LiF = {'Li': 1.0, 'F': 1.0} CO2 = {'C': 1.0, 'O': 2.0} # 辞書形式データをtransformメソッドに受け渡す df = cal.transform([H2O, LiF, CO2]) df.index = ["H2O", "LiF", "CO2"] df Structural descriptors（構造に関する記述子） XenonPyでは構造に関する特徴量も計算可能です。 from xenonpy.descriptor import Structures cal = Structures() cal 出力結果 Structures: |- structure: | |- RadialDistributionFunction | |- OrbitalFieldMatrix Structuresには、RadialDistributionFunction(RDF:動径分布関数)とOrbitalFieldMatrix(OFM)の2つのfeaturizerが含まれていることが分かります。Compositional descriptorsと同様にサンプルデータに対して適用してみます。 descriptor = cal.transform(samples) descriptor 公式ドキュメントを参照するとmatminerと互換性があると記述があります。 まとめ XenonPyのチュートリアルのDescriptor calculationについて解説しました。 続編はまた投稿予定です。