リーマン多様体上の無制約最適化問題に対する共役勾配法によるアプローチのちょっとした歴史 (〜2020)． また，その Python での実装． 対象とする問題 リーマン多様体 $\mathcal{M}$ 上の無制約最適化問題 $$\mathop{\text{minimize}}_{x\in\mathcal M}\quad f(x)$$ を考える．ここで，$f:\mathcal M \to \mathbb R$ は滑らかな写像とする． リーマン多様体上の共役勾配法 [Absil et al. 2009] リーマン多様体 $\mathcal M$ 上の共役勾配法のアルゴリズムは， Absil らの多様体上の最適化の書籍 "Optimization Algorithms on Matrix Manifolds" にて紹介されている． $x^0$ を初期化する． $d^0=−\mathop{\mathrm{grad}} f(x^0)$ で探索方向を初期化する． for $k=0,1,\ldots$ do ステップ幅 $t_k$ を計算する． $x^{k+1}=R_{x^k}\left(t_kd^k\right)$ で更新する． $\beta_k=\frac{\left\langle\mathop{\mathrm{grad}}f(x^{k+1}),\mathop{\mathrm{grad}}f(x^{k+1})\right\rangle_{{x^{k+1}}}}{\left\langle\mathop{\mathrm{grad}}f(x^{k}),\mathop{\mathrm{grad}}f(x^{k})\right\rangle_{{x^{k}}}}$ を計算する． $d^{k+1}=−\mathop{\mathrm{grad}} f(x^{k+1})+\beta_k{\mathcal T_{t_kd^k}}(d^k)$ で探索方向を更新する． end for この手法については，以下の過去記事で詳細に述べている． 線型方程式と共役勾配法 - Qiita 【Python】リーマン多様体上の共役勾配法 - Qiita リーマン多様体上の共役勾配法 [Sato & Iwai 2015] 上記で紹介したリーマン多様体上の共役勾配法は，ベクトル輸送が非増加 \left\|\mathcal T_{t_kd^k}(d^k)\right\|_{x^{k+1}} \leq \left\|d^k\right\|_{x^k} である場合に収束することが証明されている [Ring & Wirth 2012]．この条件を常に満たすように次の係数 $s_k$ を導入した共役勾配法が提案された [Sato & Iwai 2015]． s_k = \min\left\{\frac{\left\|d^k\right\|_{x^k}}{\left\|\mathcal T_{t_kd^k}(d^k)\right\|_{x^{k+1}}},1\right\} $x^0$ を初期化する． $d^0=−\mathop{\mathrm{grad}} f(x^0)$ で探索方向を初期化する． for $k=0,1,\ldots$ do ステップ幅 $t_k$ を計算する． $x^{k+1}=R_{x^k}\left(t_kd^k\right)$ で更新する． $s_k$ を計算する． $\beta_k=\frac{\left\langle\mathop{\mathrm{grad}}f(x^{k+1}),\mathop{\mathrm{grad}}f(x^{k+1})\right\rangle_{x^{k+1}}}{\left\langle\mathop{\mathrm{grad}}f(x^{k}),\mathop{\mathrm{grad}}f(x^{k})\right\rangle_{x^{k}}}$ を計算する． $d^{k+1}=−\mathop{\mathrm{grad}} f(x^{k+1})+\beta_ks_k\mathcal T_{t_kd^k}(d^k)$ で探索方向を更新する． end for リーマン多様体上の共役勾配法 [Zhu & Sato 2020] $\mathbb R^n$ 上での共役勾配法の変数の更新式は $x^{k} = x^{k+1} - t_kd^k$ と変形できる． この式をあるレトラクション $R^\text{bw}$ を用いてリーマン多様体上へ拡張すると， $$x^k = R^\text{bw}_{x^{k+1}}(-t_kd^k)$$ とかける．これより，$d^k\in T_{x^k}\mathcal M$ の $T_{x^{k+1}}\mathcal M$ 上の自然な表現として， d^k \approx -t_k^{-1}{R^\text{bw}_{x^{k+1}}}^{-1}(x^k) が得られる．このとき，ベクトル輸送の非増加条件をみたすための係数 $s_k$ は， s_k = \min\left\{\frac{\left\|t_kd^k\right\|_{x^k}}{\left\|{R^\text{bw}_{x^{k+1}}}^{-1}(x^k)\right\|_{x^{k+1}}},1\right\} となる．以上より，2つのレトラクション $R^\text{fw}$, $R^\text{bw}$ を用いた次の共役勾配法が得られる． $x^0$ を初期化する． $d^0=−\mathop{\mathrm{grad}} f(x^0)$ で探索方向を初期化する． for $k=0,1,\ldots$ do ステップ幅 $t_k$ を計算する． $x^{k+1}=R^\text{fw}_{x^k}\left(t_kd^k\right)$ で更新する． $s_k$ を計算する． $\beta_k=\frac{\left\langle\mathop{\mathrm{grad}}f(x^{k+1}),\mathop{\mathrm{grad}}f(x^{k+1})\right\rangle_{x^{k+1}}}{\left\langle\mathop{\mathrm{grad}}f(x^{k}),\mathop{\mathrm{grad}}f(x^{k})\right\rangle_{x^{k}}}$ を計算する． $d^{k+1}=−\mathop{\mathrm{grad}} f(x^{k+1})-\beta_ks_kt_k^{-1}{R^\text{bw}_{x^{k+1}}}^{-1}(x^k)$ で探索方向を更新する． end for Python での実装 使用例 問題 次の Brockett cost function の最小化問題を考える． \begin{align} \mathop{\mathrm{minimize}}_{X\in \mathbb R^{n \times r}} &\quad f(X) :=\mathop{\mathrm{tr}}\left(X^\top A X N\right)\\ \text{subject to} &\quad X^\top X = I_r \end{align} ここで，$A\in\mathbb{R}^{n\times n}$ は対称行列，$N\in\mathbb{R}^{r\times r}$ は正値の対角行列 $N = \mathop{\mathrm{diag}}(\mu_1, \mu_2, \ldots, \mu_r)$ で，$0<\mu_1\leq\mu_2\leq\dots\leq\mu_r$ をみたす．最適解 $X_*$ は $A$ の小さい順から $r$ 個の固有値に対応する固有ベクトルを並べた行列となっている．また，ユークリッド空間における $f$ の勾配は， $$\nabla f(X) = 2AXN$$ で計算される． この最適化問題は，シュティーフェル多様体 \mathop{\mathrm{St}}(n, r) := \left\{X\in\mathbb{R}^{n\times r} \mid X^\top X = I_r\right\} 上の無制約最適化問題 \mathop{\mathrm{minimize}}_{X \in \mathop{\mathrm{St}}(n, r)} \quad f(X) へと変換することができる．このとき，目的関数 $f$ のシュティーフェル多様体上での勾配 $\mathop{\mathrm{grad}}f(X)$ は，点 $X$ における接空間 $T_X \mathop{\mathrm{St}}(n, r)$ への直交射影 $\mathrm{P}_{T_X \mathop{\mathrm{St}}(n, r)}$ を用いて， \mathop{\mathrm{grad}}f(X) = \mathrm{P}_{T_X \mathop{\mathrm{St}}(n, r)}(\nabla f(X)) で計算できることが知られている．ここで，接空間 $T_X \mathop{\mathrm{St}}(n, r)$ は T_X \mathop{\mathrm{St}}(n, r) = \left\{U\in\mathbb R^{n\times r}\mid X^\top U + U^\top X = O\right\} であるので，直交射影 $\mathrm{P}_{T_X \mathop{\mathrm{St}}(n, r)}$ は \mathrm{P}_{T_X \mathop{\mathrm{St}}(n, r)}(U) = U - \frac12 X \left(X^\top U + U^\top X\right) = U - X\mathop{\mathrm{sym}}\left(X^\top U\right) となる．これより，目的関数 $f$ のシュティーフェル多様体上での勾配 $\mathop{\mathrm{grad}}f(X)$ は， \mathop{\mathrm{grad}} f(X) = 2AXN - XX^\top AXN - XNX^\top AX で計算できる．詳細は "Optimization Algorithms on Matrix Manifolds" の 4.8 節にある． シュティーフェル多様体上の演算 内積 \left\langle U, V\right\rangle_X = \mathop{\mathrm{tr}}\left(U^\top V\right) レトラクション \begin{align} R^{\mathrm{qr}}_X(U) &= \mathop{\mathrm{qf}}(X+U)\\ R^{\mathrm{ct}}_X(U) &= \left(I_n - \frac12 \left(P_X UX^\top - XU^\top P_X\right)\right)^{-1} \left(I_n + \frac12 \left(P_X UX^\top - XU^\top P_X\right)\right) X \end{align} ここで，$\mathop{\mathrm{qf}}(X)$ は $X$ の QR 分解 $X = \mathcal{Q}\mathcal{R}$ の $\mathcal Q$ 行列を返す関数 ($\mathcal R$ の対角成分は正)，$P_X = I_n - XX^\top$ である． $R^\mathrm{ct}_X(U)$ は逆写像の計算が可能であり， {R^\mathrm{ct}_X}^{-1}(Y) = 2Y\left(I_r + X^\top Y\right)^{-1} + 2X\left(I_r + Y^\top X\right)^{-1} - 2X. ベクトル輸送 \mathcal{T}_U(V) = \mathrm{P}_{T_{R_X(U)} \mathop{\mathrm{St}}(n, r)}(V) 実行例 対称行列 $A$ を以下で生成する． import numpy as np n = 10 r = 3 np.random.seed(0) A = np.random.randn(n, n) A = 0.5 * (A + A.T) 上記で紹介した ２ 通りの共役勾配法による Brockett cost function を最小化するクラスを Brockett, BrockettInvR で実装した．詳細な実装はページ末に記載する． 以下のコードで，初期点を $X_0 = \left(I_r, O\right)^\top$ として実行した． x0 = np.eye(n, r) brockett = Brockett(A=A, r=r, max_iter=1000, tol=1e-6) xopt = brockett.optimize(x0, verbosity=1) # Terminated after 133 iterations. brockett_invr = BrockettInvR(A=A, r=r, max_iter=1000, tol=1e-6) xopt_invr = brockett_invr.optimize(x0, verbosity=1) # Terminated after 139 iterations. 得られた解の正しさは，numpy.linalg.eigh の結果と比較することで確かめられる． _, V = np.linalg.eigh(A) print(np.sum(xopt * V[:, :r][:, ::-1], axis=0)) # [1. 1. 1.] _, V = np.linalg.eigh(A) print(np.sum(xopt_invr * V[:, :r][:, ::-1], axis=0)) # [1. 1. 1.] 実装 リーマン多様体の抽象クラス from abc import ABCMeta, abstractmethod from typing import Union import numpy as np class Manifold(metaclass=ABCMeta): @abstractmethod def inner_product(self, x: np.ndarray, u: np.ndarray, v: np.ndarray) -> float: raise NotImplementedError('The function retraction is not implemented') @abstractmethod def norm(self, x: np.ndarray, u: np.ndarray) -> float: raise NotImplementedError('The function norm is not implemented') @abstractmethod def retraction(self, x: np.ndarray, v: np.ndarray) -> np.ndarray: raise NotImplementedError('The function retraction is not implemented') def distance(self, x: np.ndarray, y: np.ndarray) -> Union[float, np.ndarray]: raise NotImplementedError('The function distance is not implemented') def exp(self, x: np.ndarray, v: np.ndarray) -> np.ndarray: raise NotImplementedError('The function exp is not implemented') def log(self, x: np.ndarray, y: np.ndarray) -> np.ndarray: raise NotImplementedError('The function log is not implemented') def parallel_transport(self, x: np.ndarray, y: np.ndarray, u: np.ndarray) -> np.ndarray: raise NotImplementedError( 'The function parallel_transport is not implemented') def vector_transport(self, x: np.ndarray, y: np.ndarray, u: np.ndarray) -> np.ndarray: raise NotImplementedError( 'The function vector_transport is not implemented') def inverse_retraction(self, x: np.ndarray, y: np.ndarray) -> np.ndarray: raise NotImplementedError( 'The function inverse_retraction is not implemented') def gradient(self, x: np.ndarray, d: np.ndarray) -> np.ndarray: raise NotImplementedError('The function gradient is not implemented') シュティーフェル多様体のクラス import numpy as np class Stiefel(Manifold): def inner_product(self, x: np.ndarray, u: np.ndarray, v: np.ndarray) -> float: return np.vdot(u, v) def norm(self, x: np.ndarray, u: np.ndarray) -> float: return np.linalg.norm(u) def retraction(self, x: np.ndarray, v: np.ndarray) -> np.ndarray: q, r = np.linalg.qr(x + v) return q * np.sign(np.sign(np.diag(r)) + .5) def vector_transport(self, x: np.ndarray, y: np.ndarray, u: np.ndarray) -> np.ndarray: ytu = y.T @ u return u - y @ (ytu + ytu.T) * 0.5 def inverse_retraction(self, x: np.ndarray, y: np.ndarray) -> np.ndarray: xty = x.T @ y inv_matrix = np.linalg.inv(np.identity(xty.shape[0]) + xty) return 2 * y @ inv_matrix + 2 * x @ inv_matrix.T - 2 * x 共役勾配法のクラス from abc import ABCMeta, abstractmethod from typing import Tuple import numpy as np class ConjugateGradient(metaclass=ABCMeta): def __init__( self, manifold: Manifold, initial_step: float = 1.0, armijo_param: float = 1e-4, max_iter: int = 300, tol: float = 1e-6, min_step_size: float = 1e-16, extended_output: bool = False ): self.manifold = manifold self.initial_step = initial_step self.armijo_param = armijo_param self.max_iter = max_iter self.tol = tol self.min_step_size = min_step_size self.extended_output = extended_output self.objective_values_ = [] @abstractmethod def _f(self, x: np.ndarray) -> float: raise NotImplementedError('The function _f is not implemented') @abstractmethod def _df(self, x: np.ndarray) -> np.ndarray: raise NotImplementedError('The function _df is not implemented') def _beta(self, x1: np.ndarray, x2: np.ndarray, df1: np.ndarray, df2: np.ndarray, d: np.ndarray) -> float: inner1 = self.manifold.inner_product(x1, df1, df1) inner2 = self.manifold.inner_product(x2, df2, df2) return inner2 / inner1 def _stopping_criterion(self, x: np.ndarray) -> float: d = self.manifold.norm(x, self._df(x)) return d def _step_size(self, x: np.ndarray, d: np.ndarray, g: float, fx: float, fxold: float = None) -> Tuple[float, np.ndarray]: """ Armijo condition with back tracking. """ if fxold is not None: t = 2 * (fx - fxold) / g t *= 2 else: t = self.initial_step / self.manifold.norm(x, d) res = self.manifold.retraction(x, t * d) while self._f(res) > fx + self.armijo_param * t * g: t *= 0.5 res = self.manifold.retraction(x, t * d) # break when step size 't' becomes too small if t <= self.min_step_size: return 0, x return t, res def _vector_transport(self, x1: np.ndarray, x2: np.ndarray, d: np.ndarray, step_size: float) -> np.ndarray: v = self.manifold.vector_transport(x1, x2, d) norm_d = self.manifold.norm(x1, d) norm_v = self.manifold.norm(x2, v) s = norm_d / norm_v return v if s > 1 else s * v def optimize(self, x: np.ndarray, verbosity: int = 0) -> np.ndarray: fx = self._f(x) # store objective function value if necessary if self.extended_output: self.objective_values_.append(fx) # first iteration z = x.copy() df1 = self._df(z) d = -df1.copy() g = self.manifold.inner_product(x, df1, d) step_size, res = self._step_size(z, d, g, fx, None) fxold = fx fx = self._f(res) # store objective function value if necessary if self.extended_output: self.objective_values_.append(fx) # main loop i = 0 for i in range(self.max_iter - 1): # direction df2 = self._df(res) beta = self._beta(z, res, df1, df2, d) d = -df2 + beta * self._vector_transport(z, res, d, step_size) g = self.manifold.inner_product(res, df2, d) if g >= 0: if verbosity >= 2: print(f'Got an ascent direction (g = {g}), ' 'reset the search direction') d = -df2.copy() g = self.manifold.inner_product(res, df2, d) df1 = df2.copy() # update z = res.copy() step_size, res = self._step_size(res, d, g, fx, fxold) fxold = fx fx = self._f(res) # store objective function value if necessary if self.extended_output: self.objective_values_.append(fx) if step_size == 0: if verbosity >= 2: print(f'Got a zero step size.') break # break if convergence if self._stopping_criterion(res) < self.tol: break if verbosity >= 1: print(f'Terminated after {i + 2} iterations.') return res class ConjugateGradientInvR(ConjugateGradient, metaclass=ABCMeta): @abstractmethod def _f(self, x: np.ndarray): raise NotImplementedError('The function _f is not implemented') @abstractmethod def _df(self, x: np.ndarray): raise NotImplementedError('The function _df is not implemented') def _vector_transport(self, x1: np.ndarray, x2: np.ndarray, d: np.ndarray, step_size: float) -> np.ndarray: v = - self.manifold.inverse_retraction(x2, x1) / step_size norm_d = self.manifold.norm(x1, d) norm_v = self.manifold.norm(x2, v) s = norm_d / norm_v return v if s > 1 else s * v Brockett cost function のソルバークラス import warnings import numpy as np from sklearn.utils.validation import check_symmetric class Brockett(ConjugateGradient): def __init__( self, A: np.ndarray, r: int, initial_step: float = 1.0, armijo_param: float = 1e-4, max_iter: int = 300, tol: float = 1e-6, min_step_size: float = 1e-16, extended_output: bool = False ): if A.ndim != 2: raise ValueError( f'Expected 2D array, got {A.ndim}D array instead.') super().__init__( manifold=Stiefel(), initial_step=initial_step, armijo_param=armijo_param, max_iter=max_iter, tol=tol, min_step_size=min_step_size, extended_output=extended_output ) n = A.shape[0] if r > n: warnings.warn(f'Given r = {r} is larger than {n}. ' f'Use r = {n} instead.') r = n self.A = check_symmetric(A) self.N = np.arange(1, r + 1) / r def _f(self, x: np.ndarray) -> float: return np.sum((self.A @ x) * (x * self.N)) def _df(self, x: np.ndarray) -> np.ndarray: AX = self.A @ x AXN = AX * self.N return 2 * AXN - x @ (x.T @ AXN) - x @ (AXN.T @ x) class BrockettInvR(ConjugateGradientInvR): def __init__( self, A: np.ndarray, r: int, initial_step: float = 1.0, armijo_param: float = 1e-4, max_iter: int = 300, tol: float = 1e-6, min_step_size: float = 1e-16, extended_output: bool = False ): if A.ndim != 2: raise ValueError( f'Expected 2D array, got {A.ndim}D array instead.') super().__init__( manifold=Stiefel(), initial_step=initial_step, armijo_param=armijo_param, max_iter=max_iter, tol=tol, min_step_size=min_step_size, extended_output=extended_output ) n = A.shape[0] if r > n: warnings.warn(f'Given r = {r} is larger than {n}. ' f'Use r = {n} instead.') r = n self.A = check_symmetric(A) self.N = np.arange(1, r + 1) / r def _f(self, x: np.ndarray) -> float: return np.sum((self.A @ x) * (x * self.N)) def _df(self, x: np.ndarray) -> np.ndarray: AX = self.A @ x AXN = AX * self.N return 2 * AXN - x @ (x.T @ AXN) - x @ (AXN.T @ x) 参考までに，Pymanopt にも共役勾配法の実装がある．使用例については過去記事を参照されたい． 【Python】リーマン多様体上の共役勾配法 - Qiita 参考文献 Absil, P.-A., Mahony, R., & Sepulchre, R. (2009). Optimization algorithms on matrix manifolds. Princeton University Press. Ring, W., & Wirth, B. (2012). Optimization methods on Riemannian manifolds and their application to shape space. SIAM Journal on Optimization, 22(2), 596-627. Sato, H., & Iwai, T. (2015). A new, globally convergent Riemannian conjugate gradient method. Optimization, 64(4), 1011-1031. Zhu, X., & Sato, H. (2020). Riemannian conjugate gradient methods with inverse retraction. Computational Optimization and Applications, 77(3), 779-810.

yukicoder contest 293 参戦記 A 1491 銀将 さっぱり分からないので、まずはナイーブなコードを書いてみる. from itertools import product def f(n): a = [(-1, -1), (0, -1), (1, -1), (-1, 1), (1, 1)] result = set() for i in range(n + 1): for p in product(a, repeat=i): x, y = 0, 0 for xx, yy in p: x += xx y += yy result.add((x, y)) return len(result) 実際に実行してみる. >>> f(0) 1 >>> f(1) 6 >>> f(2) 18 >>> f(3) 38 >>> f(4) 66 >>> f(5) 102 >>> f(6) 146 >>> f(7) 198 >>> f(8) 258 >>> f(9) 326 どうやら f(n) = f(n-1) + (f(n-1) - f(n-2) + 8) になっているように見える. Wolfram|Alpha に f(n)=f(n-1)+(f(n-1)-f(n-2)+8), f(2)=18, f(3)=38 を入力したら f(n) = 4*n*2+2 と教えてくれたので解けた. K = int(input()) if K == 0: print(1) else: print(4 * K * K + 2)

PythonのWebアプリ開発用フレームワーク「Flask」の初歩の初歩、Hello, Worldをブラウザで表示するまでをまとめました。 Flaskクラスのインスタンスを作成 まずFlaskでアプリケーションを作る際にはFlaskクラスのインスタンスを作成します。 from flask import Flask app = Flask(__name__) Flaskクラスは変数にこのインスタンスがあるモジュールの名前をとります。 シンプルな（.pyファイルが１つしかないような）アプリケーションの場合は上記で良いのですが、 そうでない場合はモジュール名を渡す必要があるので、 app = Flask(__name__.split('.')[0]) などとすることがドキュメントでは奨励されています。 いくつかの拡張ライブラリではデバック時にappの正しい場所を変数として渡していないとライブラリ側で検知されず、デバックに苦労するみたいです。 （ドキュメントではFlask_SQLAlchemyなどが拡張パッケージの例として挙げられています） ルーティングとView Function 通常のWebアプリケーションではクライアントサイド（主にブラウザ）とサーバーがリクエスト/レスポンスのやりとりをしますが、Flaskではサーバー側を前項で作ったappが担います。 そのため、作成したappに何のリクエストが来た場合にどんな処理をするのか、を明示する必要があります。 このうち、「何のリクエスト」を示すのがルーティング、「どんな処理」を示すのがview functionです。 ルーティングはappのrouteメソッドを使って、デコレータとして記述し、View Functionはその名の通り関数として定義します。 from flask import Flask app = Flask(__name__) @app.route('/') def index(): return '<h1>Hello, World.</h1>' routeメソッドは引数にURLを取り、マッピングを行います。 上記ではルート（'/'）にアクセスした際にhtmlのh1タグをつけてHello, World.と表示するように処理してます。 ※リクエストのメソッドを指定できたり、URLにクエリパラメータを想定したりもできるのですが、一旦ここでは割愛します。 ローカル環境でFlaskアプリを起動し、Webブラウザで表示する Flaskでの処理は以上で、最後にローカル環境でアプリを起動すればOKです。 ※今更ですがFlaskはすでにインストール済みの想定です。 先ほど記述した.pyファイルをapp.pyとし、そのディレクトリで環境変数にてアプリファイル名を指定し、 flask runでアプリを起動します。 $ export FLASK_APP=app.py $ flask run * Serving Flask app "app.py" * Environment: production WARNING: This is a development server. Do not use it in a production deployment. Use a production WSGI server instead. * Debug mode: off * Running on http://127.0.0.1:5000/ (Press CTRL+C to quit) 上記の場合はhttp://127.0.0.1:5000/で起動しているので、ブラウザからアクセスすればきっとHello, World.と表示されていることでしょう。

ゴミ出しを忘れないよう、とりあえず作ってみた。 cron で1日1回一ヶ月間動かしてみて様子見。 ※url、token は伏せています。 #!/usr/bin/python3 # 市ゴミの日カレンダーの HTML を解析して # 今日が何のゴミの日かを LINE Notify に通知する # 該当なしの場合何もしない # import subprocess as sp import re import requests import datetime from bs4 import BeautifulSoup import urllib.request as req from string import digits SECRET = XXX def send_line_notify(notification_message): line_notify_token = 'yXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX' line_notify_api = 'https://notify-api.line.me/api/notify' headers = {'Authorization': f'Bearer {line_notify_token}'} data = {'message': f'message: {notification_message}'} requests.post(line_notify_api, headers = headers, data = data) yyyymm = str(datetime.datetime.today().strftime("%Y-%m")) url = "http://www.hoge.city.jp/gomical/calendar.html?date=" + yyyymm + "&id=" + str(SECRET) + "#gaba ge-calendar" response = req.urlopen(url) parse_html = BeautifulSoup(response,"html.parser") print (parse_html.title.string) table = parse_html.findAll("table",{"class":"calendar"})[0] # calendar class を探す rows = table.findAll("tr") dd = str(datetime.datetime.today().strftime("%d")) for row in rows: for cell in row.findAll(['td','th']): # カレンダーテーブルを読み込む if len(cell.get_text()) > 3: # ゴミの日以外は日付2文字しかない cal_dd = cell.get_text() m = re.search(r'\d+', cal_dd) # 1日はもやすごみ の1 を取得する if dd == m.group(): # dd (Today) と、カレンダーの日付が同じ（今日がゴミの日なら） tbl = str.maketrans('','', digits) # カレンダーから取得した日付を除去する str = cal_dd.translate(tbl) msg = "今日は" + str + "の日です" rc = sp.run("/var/www/html/alexa.sh -e 'speak:" + msg + "'", shell=True, stdout=True, stderr=sp.PIPE, text=True).returncode send_line_notify(msg) #end of Gabage.py

金融工学の本を読むと，しばしばボラティリティ・スマイルについて言及されています．ただ，自分の手で計算してみないとわかった気にならない性質なので，そもそもボラティリティ・スマイルとは何かから始め，実際にグラフで確認してみました．なお，すべての実装はgithub (https://github.com/yotapoon/volatility_smile )にアップしているので，実装だけ確認したい人はこちらをご覧ください． ボラティリティ・スマイルとは まず問題設定を明確にし，いくつか用語を定義してからボラティリティ・スマイルを説明します．問題設定は端折っている部分も多いため， 知識が足りない方はぜひhttps://qiita.com/yotapoon/items/30c1bc411c5a8f08d219 を参照してください．ご存知の方は飛ばしていただいて構いません． 問題設定 いま注目する原資産の価格が$S(t)$なる確率過程であらわされるものとします．期待収益率を$\alpha$，ボラティリティ(要するに変動率)を$\sigma$とすれば，$S(t)$は以下の幾何ブラウン運動に従うものとモデル化できます： $$ dS(t) = \alpha S(t) dt + \sigma S(t) dW(t) $$ ただし$W(t)$はブラウン運動で，$dW(t)$は平均$0$，標準偏差$\sqrt{dt}$の正規分布に従うといえます． また，単位時間当たり$r$の成長率を持つ安全資産の存在を仮定します．$r$は金利と考えることができます． 次に，ヨーロピアン・コールオプションを定義します．ヨーロピアン・コールオプションを購入するとは，「時刻$t = T$において，原資産の価格$S(T)$によらず，行使価額$K$で原資産を1単位取得できる"権利"を購入する」ということです．数式にまとめると，時刻$t = T$において $$ c(T, S(T)) = (S(T) - K)^+ $$ という価格となるような商品だと考えられます．ただし，$(x)^+$は$x \geq 0$のとき$x$を返し，$x < 0$のとき$0$を返す関数です． 同様にヨーロピアン・プットオプションは，時刻$t = T$において $$ p(T, S(T)) = (K - S(T))^+ $$ という価格となるような商品だといえます． $t = 0$時点においてこれらの商品を売買するためには，$c(t=0, S(0))$や$p(t=0, S(0))$の値(つまり価格)を決める必要があります．無裁定価格理論を用いると，一般に$t = 0$の時点以外においても価格を決定することができ， c(t, x) = x \Phi(d_+(T-t, x))) - K e^{-r(T-t)} \Phi(d_-(T-t, x)) \\ p(t, x) = -x \Phi(-d_+(T-t, x))) + K e^{-r(T-t)} \Phi(-d_-(T-t, x)) という形となります．ただし， d_{\pm} (\tau, x) = \dfrac{1}{\sigma \sqrt{\tau}} \left[ \log \dfrac{x}{K} + \left( r\pm \dfrac{1}{2} \sigma^2 \right) \tau \right] \\ \Phi(x) = \dfrac{1}{\sqrt{2\pi}} \int_{-\infty}^x e^{-x'^2/2} dx' を導入しました． インプライド・ボラティリティ コールオプションの価格$c(t, S(t))$を決定する式をよく見ると，計算するために必要な量は①満期までの時間$T-t$，②原資産価格$S(t)$，③行使価額$K$，④ボラティリティ$\sigma$，⑤金利$r$であることがわかります．この中で，④ボラティリティ$\sigma$の値は，契約時点$t = 0$において正確にはわからず，推定するしかありません．一方で，コールオプションの価格$c(0, S(0))$は市場を観察することで決められます．ゆえに，コールオプションの価格$c(t, S(t))$を決定する式の左辺に，市場で観察された値($c^*$としましょう)を代入することで，ボラティリティ$\sigma$を推定することができます．つまり，以下の式を$\sigma$に関する方程式として解くのです． c^* = x \Phi(d_+(T, x))) - K e^{-rT} \Phi(d_-(T, x)) \\ d_{\pm} (T, x) = \dfrac{1}{\sigma \sqrt{T}} \left[ \log \dfrac{x}{K} + \left( r\pm \dfrac{1}{2} \sigma^2 \right) T \right] \\ (細かい人のための追記)細かい人は解の一意性などが気になるかもしれませんが，$c(t, S(t))$が$\sigma$に関して単調増加であるという性質を考えると，解の一意性は直ちにわかります．さらに，これはニュートン法などのアルゴリズムが適用可能であることも示唆しています． 以上の手続きにより推定されたボラティリティ$\sigma$を，インプライド・ボラティリティと呼びます．これは，市場により「暗示された(implied)」ボラティリティであるということを意味しています． ボラティリティ・スマイル 長くなってしまいましたが，ようやく本題であるボラティリティ・スマイルを説明できます．前述のインプライド・ボラティリティ$\sigma$を，行使価額$K$の関数として描いてみたものが下図です．計算の詳細は後述します． グラフから，明らかに凸関数の形状をしていることが確認できます．これが人の笑っている顔に見えることから，ボラティリティ・スマイルという名前がついています(いうほど見えるか？)．もともとのモデルでは，原資産のボラティリティ$\sigma$は行使価額と全く無関係に決まっていたにも関わらず，このような非自明な関係が生じていることに注意します． さらに，①計算時点(2021年4月28日)での原資産価格(日経平均)$S(0) = 29053.97$付近から離れるほど，インプライド・ボラティリティが大きくなっていること，および②全体的に右下がりになっていることも重要です．これらの性質は，満期が短い株式オプションに典型的にみられるものです． ボラティリティ・スマイルの原因は，「市場参加者が想定している原資産価格$S(t)$の変動は，ブラックショールズモデルで仮定したものと一致していないため」だそうです．より詳しい解説をQuant Collegeの中の方がされているので，ぜひご覧ください( https://quantcollege.net/smile-pricing-volatility-smile-explained )． インプライド・ボラティリティの計算 まず，インプライド・ボラティリティを計算するための手続き，及びニュートン法によるアルゴリズムを示します．また，日本取引所の公開しているデータ( https://www.jpx.co.jp/markets/derivatives/option-price/index.html )を使って，実際にボラティリティ・スマイルを確認します． データの確認 日本取引所の公開しているデータから，日経225の2021年6限月(つまり原資産を日経225とし，2021年6月の第二金曜日に満期を迎えるオプション)の，2021年4月28日における終値を抽出しました．データの整形などは本題ではないため割愛します．行使価額$K$と，コール・プットそれぞれの価格の関係を以下に示します． なお，データによると原資産価格は$S(0) = 29053.97$，基準ボラティリティ(?)は$\sigma = 0.1544$だそうです．$\sigma = 0.1544$としてオプション価格をプロットしたものも同時に示しました．コールオプションに関しては良い近似になっているものの，プットオプションでは乖離が生じています．また，$\sigma$の値を変えてみても，プットオプションの乖離は改善されなかったことにも触れておきます． 目標は，行使価額$K$の値ごとに決まる市場価格$c^*$を用いて，インプライド・ボラティリティ$\sigma$を計算することです． ニュートン法によるインプライド・ボラティリティの計算 まず，コールオプションに関してインプライド・ボラティリティ$\sigma$をどのように求めるか考えてみます．そのために，コールオプションの価格$c(t, S(t))$を$c(K, \sigma)$のように書き直しておきます．いまの目標は，行使価額$K$とそれに対応する市場価格$c^*$から $$ c^* = c(K, \sigma) $$ を満たすようなボラティリティ$\sigma$を求めることです． ここで，$c(K, \sigma)$を$\sigma$で微分した値 $$ \nu(K, \sigma) = \dfrac{\partial c}{\partial \sigma} = S(0) \sqrt{T} \dfrac{1}{\sqrt{2 \pi}} \exp \left[ -\dfrac{1}{2} d_+^2 \right] > 0 $$ に注目します．ベガと呼ばれるこの値は明らかに正であり，ゆえに$c(K, \sigma)$はボラティリティ$\sigma$に対して単調増加です．よって，以下に示すニュートン法で数値解を得ることができます： $$ \sigma_t = \sigma_{t-1} - \dfrac{c(K, \sigma_{t-1}) - c^*}{\nu(K, \sigma_{t-1} )} $$ ただし，初期値$\sigma_0$が小さいとベガ$\nu(K, \sigma)$が$0$に近づいてしまい，うまく計算できないことがあるため，$\sigma_0$は大きくとることをお勧めします． なお，プットオプションについてもベガは同じ値を取ることが示せます(これはプット・コールパリティから直ちにわかります)．そのためほとんど同様の手続きで，プットオプションの市場価格からもインプライド・ボラティリティを計算することができます． 実装 インプライド・ボラティリティを計算する手順がわかったので，実装します．get_implied_volがインプライド・ボラティリティを計算するメソッドです． import numpy as np from scipy.stats import norm S0 = 29053.97 # initial stock price r = 0.00 # interest rate T = 44 / 365.0 # 1 year class European_call: def __init__(self, S0, r, T): self.S0 = S0 self.r = r self.T = T def get_d_list(self, K, vol): d_plus = (1.0 / (vol*np.sqrt(self.T))) * (np.log(self.S0/K) + (self.r+0.5*vol*vol)*self.T) d_minus = (1.0 / (vol*np.sqrt(self.T))) * (np.log(self.S0/K) + (self.r-0.5*vol*vol)*self.T) return [d_plus, d_minus] def get_value(self, K, vol): d_plus, d_minus = self.get_d_list(K, vol) return self.S0*norm.cdf(d_plus) - K*np.exp(-self.r*self.T)*norm.cdf(d_minus) def get_vega(self, K, vol): d_plus, d_minus = self.get_d_list(K, vol) return self.S0*np.sqrt(T)*(1.0/np.sqrt(2.0*np.pi))*np.exp(-0.5*d_plus*d_plus) def get_implied_vol(self, c_actual, K, num_iteration = 80): # c_actual:コールオプションの市場価格 implied_vol = 1.6 # 小さい値を取るとvegaが0になってしまい，ニュートン法に支障をきたす． for idx in range(num_iteration): # ニュートン法により，オプションの理論値と市場価格が等しくなるようなsigma_impliedを求める． implied_vol = implied_vol - (self.get_value(K, implied_vol)-c_actual) / self.get_vega(K, implied_vol) return implied_vol 以下の図は，このプログラムを利用して得られたインプライド・ボラティリティを示したものです．ただしcall option，put optionはそれぞれコールオプション，プットオプションの市場価格$c ^ *$，$p ^ *$から計算したインプライド・ボラティリティです．それぞれのオプションごとに，価格データが存在する行使価額の範囲が大きく異なっているため，カバーする領域も異なっています(「データの確認」の項を参照のこと)． おまけ もともとのデータに飛びがあるため，上のボラティリティ・スマイルの図は少しガタついています．せっかくなので，簡単な回帰によって図をきれいにしてみました． いい感じになっていることがわかります． なお，今回は行使価額$K$の関数としてのみ扱いましたが，満期までの時間$T$についても同様にインプライド・ボラティリティを計算することができます．このようにして得られる$\sigma(K, T)$のような平面を，ボラティリティ・サーフェスと呼びます．おそらく今回の実装を少しいじるだけでできると思うので，気が向いたらやってみます．

おはようございます。最近起きれてますか？ 私はコロナで完全に生活リズムが狂ってしまって、4時寝11時起きの生活をしています。 アラームはかけているのですがなかなか起きれない理由として、朝起きる理由がない、刺激がないのではと思い、それを解決する方法を考えてました。 その中で見つけました！それはキズナアイのおはようツイートです！！ だいたい6時から11時までの間にツイートがされていて、平日の朝は我々を励ます動画を毎日投稿してます。(じゃんけん勝負をしています) これがツイートされるタイミングで起きれば、起きれるのでは？？と思い、slackやdiscordに通知を入れてくれるBOTを作りました。 ちなみにwebhookのurlをGoogle フォームに入力したらみなさんのslackにも通知が入ります！！ 完成したもの googleフォームにwebhookURLを入力すれば使えます！ 仕様技術 Heroku Scheduler Python Tweety Google SheetsAPI Google form(入力側) 設計図 githubのリポジトリ また、このQiitaの内容は、2021/4/28に行われた#技育CAMPのLT大会で発表したアプリケーションです。 ステップ1 ツイートの取得 Heroku Schedulerは10分おきに実行をできるので、10分おきにキズナアイの新しいツイートがないかtweetyを使って検索をします。 肝となる部分は、リツイートを含まず、外部に保存しているtweetedより大きいキズナアイのツイートがあったら取得という形を取っています。 GetUserTweet.py def GetUserTweet_data( Account,tweet_id,option_time=10): ... tweets = api.user_timeline(Account, count=20, page=1) BaseURL = f"https://twitter.com/{ Account }/status/" tweet_data =[] for tweet in tweets: if tweet_id < tweet.id: TweetURL = BaseURL+str(tweet.id) tweet_text = f"{tweet.text}\n{TweetURL}" tweet:dict ={"id":tweet.id , "data":tweet_text} tweet_data.append(tweet) return tweet_data webhookで通知を送信する WebhookAction.py def WebhookApp(webhook_url,comment, username,image_URL): if "discord" in webhook_url: # Discord main_content = { 'username': username, 'avatar_url': image_URL, 'content':comment } headers = {'Content-Type': 'application/json'} response = requests.post(webhook_url, json.dumps(main_content), headers=headers) print(response.status_code) return response.status_code elif "slack" in webhook_url: # Slack response = requests.post(webhook_url, data = json.dumps({ 'text': comment, # 通知内容 'username': username, # ユーザー名 'icon_url': image_URL })) print(response.status_code) return response.status_code else: return Exception("Not found:Webuook") データベースはスプレットシート 今回、自分以外の人が支えるようにしたかったのと、フロントエンドの作成がめんどくさかったという理由からフロントエンドはGoogleフォームを使うことにしました。googleフォームに入力された結果はスプレットシートに保存することができるので、スプレットシートをデータベースみたいに利用しました。 Googleフォームを使うことで、URLの更新などもできるため、新規フォーム、更新フォームを個別に作る必要もないため簡単に作る時は便利です。 sheets.py def Auth_Sheet(credentials_name,SPREADSHEET_KEY,sheet_name,isTest=False): #2つのAPIを記述しないとリフレッシュトークンを3600秒毎に発行し続けなければならない scope = ['https://spreadsheets.google.com/feeds','https://www.googleapis.com/auth/drive'] #ダウンロードしたjsonファイル名をクレデンシャル変数に設定（秘密鍵、Pythonファイルから読み込みしやすい位置に置く） credentials = ServiceAccountCredentials.from_json_keyfile_name(credentials_name, scope) #OAuth2の資格情報を使用してGoogle APIにログインします。 gc = gspread.authorize(credentials) #共有設定したスプレッドシートキーを変数[SPREADSHEET_KEY]に格納する。 if isTest == True: SPREADSHEET_KEY = config.SPREADSHEET_TEST_KEY else: SPREADSHEET_KEY = config.SPREADSHEET_KEY #共有設定したスプレッドシートのシート1を開く Sheet = gc.open_by_key(SPREADSHEET_KEY).worksheet(sheet_name) return Sheet 上のコードは認証部分のコードです。 おはようツイートかどうかを判別する ここはそこそこ悩みました。 何も添付されていない系 おはよう！の進化系 げつようび系 総じて「おはよう」の変形なのですが、一貫性がなかったり文章なしの場合があるためおはようツイートをツイート内容から判別するのは確実ではないと感じました(機械学習で判定したかった…) そこで、他に気がついた法則性として、6時から11時までの「その日にはじめてツイートされた内容」をおはようツイートとして判別することにしました。 プログラムを起動した時間が5時代だったら、その日のツイートカウントを0にして、ツイートがあったら+1する原始的な方法を使いました。 この判定の数値もスプレットシート内に置いています。 isMorning.py def isMorningTweet(tweet,time_now,isTest=False): isMorningTweet = False #time_now = datetime.datetime.now() print(int(time_now.hour)) # 5時にリセットする！ if int(time_now.hour)== 5: Auth = SheetsAPI.sheet.Auth_Sheet(config.SheetAuthKeyName,config.SPREADSHEET_KEY,"Num",isTest) Auth.update_cell(2,1,0) Log = SheetsAPI.sheet.Auth_Sheet(config.SheetAuthKeyName,config.SPREADSHEET_KEY,"Log",isTest) Log.append_row([str(time_now.strftime("%Y/%m/%d %H:%M:%S")),f"5時代:0"]) elif 6 <= int(time_now.hour) and int(time_now.hour) < 12: # 時間での判別 print("朝時間") if tweet['data'].find("https://t.co") != -1 and tweet['data'] not in "RT": # 動画のツイート&RTじゃない Auth = SheetsAPI.sheet.Auth_Sheet(config.SheetAuthKeyName,config.SPREADSHEET_KEY,"Num",isTest) num :int = Auth.get_all_records()[0]['Num'] print(num) # その日の始まりのツイートを取得 if num == 0: isMorningTweet = True num += 1 Auth.update_cell(2, 1, num) Log = SheetsAPI.sheet.Auth_Sheet(config.SheetAuthKeyName,config.SPREADSHEET_KEY,"Log",isTest) Log.append_row([str(time_now.strftime("%Y/%m/%d %H:%M:%S")),f"最初のツイート{num}",True]) else: num += 1 Auth.update_cell(2, 1, num) Log = SheetsAPI.sheet.Auth_Sheet(config.SheetAuthKeyName,config.SPREADSHEET_KEY,"Log",isTest) Log.append_row([str(time_now.strftime("%Y/%m/%d %H:%M:%S")),f"その日の2個以上:{num}",False]) return isMorningTweet ログを取る 一応、確認のためにスプレットシート上にログを取るようにしました。 6月までに作ったファイルは無料で無制限なので、スプレットシート最強！！ という感じでキズナアイのおはようツイートを通知するBOTを作りました。 webhook URLを取得しGoogleフォームに入力することで、皆さんも使うことができます！ (ちなみに、これを作っても生活習慣は改善しなかったので、次はBOTではなくアラームアプリでも作ればなんとかなる？？？と思いました) ということで是非使ってみてください

これは競プロ形式のPython学習サービス「TOYPRO(トイプロ)」の練習問題のうちの一つ，「割り込み厳禁」 (500点問題)の解説です． 問題文（一部省略） 並んでいる人の人数$N$と、彼らがどのような順番で並んでいるのかを表したリストが与えられます。割り込んだ人を特定するプログラムを作成してください。具体的には、割り込んだ人が持っている整理券の番号を、小さい順にすべて出力してください。割り込んだ人がいない場合は"割り込んだ人はいません"と出力してください。 制約 $1 \leq N \leq 100$ 解説 順番のリストを前から見ていくと大変です．(大変とは言っていない(後述のO(N^2))) しかし，割り込んだ人の後ろには必ずそれより小さい番号がある，ということを利用すると，まず最小値を大きな値でおき，順番のリストを後ろから最小値を更新，または割り込んだ人リストに追加しながら見ていく，という方法で解くことが出来ます． そして，割り込んだ人リストを小さい順にソートし，出力します． 実装例 # 入力 N = 3 a = [1, 3, 2] m = 101 # 十分に大きな値で最小値を初期化 f = [0] * (N+1) # f[i]が1ならi番目の人は割り込んだ人，0ならそうでない人 for i in a[::-1]: # 後ろから見ていく if i > m: # 後ろに自身より小さい番号がある場合 f[i] = 1 # 割り込んだ人としてマーク else: m = i # 最小値を更新 # 出力 if 1 in f: # 割り込んだ人が1人でもいる場合 for i in range(1, N+1): if f[i]: # f[i]がTrueの場合 print(i) else: # 割り込んだ人が0人の場合 print("割り込んだ人はいません") 計算量 計算量は$O(N)$です． 他の解法として，割り込んだ人の番号そのものをメモしてソートして出力する，というようなものがありますが，ソートがボトルネックとなって計算量は$O(N \log N)$となってしまいます．が，正直大差ないのでこれでも大丈夫です． というかそもそも制約が小さいので前から見ていって$O(N^2)$でもいけます．（は？） まあ，制約が大きかったときの練習になるのでヨシ！ 終わりに 最後までご覧いただきありがとうございました． 記事のミス，直すべき箇所等があった場合，コメントまたはTwitter(@MUSUKA0605)に報告をして頂けると嬉しいです．それでは〜

結果から原因を予測する問題を逆問題と言います。 逆問題では順問題とは異なり問題が非適切になっていることが大半で、安定した解が得られません。 非適切性が原因で生じる誤差鋭敏性という現象について説明します。 環境 Windows10 Python 3.8.6 逆問題とは？ モデルの出力から入力を推定する問題のことです。 逆問題ではモデルに組み込まれていない観測値に含まれる未知量が原因で非適切な問題となることが多いです。 出力から入力を求める過程を逆過程と言い、逆問題を解く方法を逆解析と言います。 逆問題の適切性 逆問題は以下の3つの適切性を満たさないとき安定した解が得られません。 解の存在性 解の一意性 解の安定性 1 と 2 は条件を満たしている、あるいは満たしていなくても問題にならないことが多いです。 3 が満たされていないとき、誤差鋭敏性という現象が生じて解が安定しなくなります。 観測値にはノイズが含まれていますのでほとんどの逆問題で発生します。 誤差鋭敏性とは？ 観測誤差が含まれる観測値で逆解析をすると解が大きく変動する現象のことをいいます。 これは特異値が小さな特異ベクトルで生じます。特異値の詳細はwikiを参照してください。 誤差鋭敏性のメカニズムは以下の通りです。 偏回帰係数を $A$ として $r$ 番目の特異ベクトルの特異値が小さな値だったとします。観測誤差が含まれた観測値を $A$ の逆行列に入力すると観測誤差の $r$ 番目の特異ベクトル成分と特異値の逆数が掛け合わされて解が大きく変動します。 では特異値が小さいときに誤差鋭敏性が起きるのか確認してみましょう。 ... # 観測値b1から0.1ずつずらした観測値をb2として解の変動を確認します。 >>> A = np.array([[7.1, 5.9], [7.7, 6.5]]) >>> b1 = np.array([3.6, 3.8]) >>> b2 = np.array([3.5, 3.9]) >>> np.linalg.inv(A) @ b1 array([1.3611111, -1.0277778]) >>> np.linalg.inv(A) @ b2 array([-0.3611111, 1.0277778]) ... # 観測値が0.1しかずれていないにもかかわらず解が大きく変動しています。 $A$ を特異値分解すると特異値が小さいのが分かります。 >>> u, s, vh = np.linalg.svd(A) >>> s array([13.66591, 0.052686]) 誤差鋭敏性の大きさは特異値の最大値と最小値の比の条件数で決まります。 条件数が大きいときは誤差鋭敏性が大きい悪条件になります。 >>> s[0]/s[1] 259.385034 $r$ 番目の特異ベクトルを削除する低ランク近似あるいは特異値を正則化項で大きくするチホノフ正則化によって誤差鋭敏性は適切化されます。 まとめ 逆過程で解を陽に求めると誤差鋭敏性が生じます。 誤差鋭敏性は特異値の小さな特異ベクトルで生じその大きさは条件数で決まります。 悪条件になっている場合低ランク近似あるいはチホノフ正則化で適切化する必要があります。 適切化手法は別の記事に書きます。

はじめに 「株価分析（SMA) - Backtestingを使ってリターンを計算する」の続編となります。今回はRSIを使って株価のリターンをシミュレートしたいと思います。 RSIは、Relative Strength Indexの略で、日本語だと相対力指数と呼ばれています。RSIの値が高いと買われすぎ、低いと売られすぎの指標となるそうです。 さっそくRSIを計算してみます。 まずは株価の取得です。こちらは今まで通りpandas_datareaderを使用して取得しています。一度取得したものは覚えるようにしています。 def get_stock(ticker, start_date, end_date): ''' get stock data from Yahoo Finance ''' dirname = '../data' os.makedirs(dirname, exist_ok=True) period = f"{start_date.strftime('%Y%m%d')}_{end_date.strftime('%Y%m%d')}" fname = f'{dirname}/{ticker}_{period}.pkl' df_stock = pd.DataFrame() if os.path.exists(fname): df_stock = pd.read_pickle(fname) start_date = df_stock.index.max() + datetime.timedelta(days=1) if end_date > start_date: df = pandas_datareader.data.DataReader( ticker, 'yahoo', start_date, end_date) df_stock = pd.concat([df_stock, df[~df.index.isin(df_stock.index)]]) df_stock.to_pickle(fname) return df_stock 2011年〜2020年の10年間の日経平均株価を取得してRSIを計算します。RSIの計算期間は14日としています。(Talibのデフォルト値と同じ。） START_DATE = datetime.date(2011, 1, 1) END_DATE = datetime.date(2020, 12, 31) TICKER = '^N225' df = get_stock(TICKER, START_DATE, END_DATE) df['RSI'] = talib.RSI(df['Adj Close'].values, timeperiod=14) 日経平均株価と計算したRSIをグラフにしてみます。 fig, axes = plt.subplots(nrows=2, sharex=True, figsize=(10, 6)) axes[0].set_title('RSI') axes[0].plot(df.index, df['Adj Close'], label=f'{TICKER} Close') axes[0].grid() axes[0].legend(loc='upper left') axes[1].plot(df['RSI'], label='RSI 14days') axes[1].axhline(y=20, color='blue', linestyle=':') axes[1].axhline(y=30, color='red', linestyle=':') axes[1].axhline(y=70, color='red', linestyle=':') axes[1].axhline(y=80, color='blue', linestyle=':') axes[1].grid() axes[1].legend(loc='upper left') plt.show() RSIが高いと買われすぎ。つまりRSIが高くなった後は株価が下がる可能性が高くなる。 RSIが低いと売られすぎ。つまりRSIが低くなった後は株価が上がる可能性が高くなる。 80%と70%、30%と20%のところに横線を入れていますが、前者はRSIが高いと判断するときによく使われる値、後者はRSIが低いと判断するときによく使われる値です。 でも、このグラフだけではよくわかりませんね。 RSI閾値に達したときの株価変化率の推移 では、RSIが上限基準値を超えたとき、あるいは下限基準値を下回ったとき、その後の株価の変化はどうなっているのでしょうか？ RSIの計算期間を5日〜40日、RSI上限値を60%〜89%、RSI下限値を10%〜39%の総当たりで確認してみます。RSIの上限値あるいは下限値を超えたときから30日間の株価の平均変化率を求めます。30日間の計算は歴日で行っているので、実際の日数は休日分だけ少なくなってます。計算結果はdf_upper、df_lowerに保存します。 MAX_TIMEPERIOD = 40 df_stock = get_stock(TICKER, START_DATE, END_DATE) df_upper = pd.DataFrame() df_lower = pd.DataFrame() for timeperiod in np.arange(5, MAX_TIMEPERIOD + 1): df = df_stock.copy() df['RSI'] = talib.RSI(df['Adj Close'].values, timeperiod=timeperiod) df = df.dropna() # RSI上限 results = [] for upper_limit in np.arange(60, 90): # RSIが上限を超えた日付を探す dates = df[(df['RSI'].shift() < upper_limit) & (df['RSI'] >= upper_limit)].index # RSIが上限を超えた日を基準に株価の変化率を求める changes = calc_change(df, dates) # 結果に追加 results.append( {'limit': upper_limit, 'average': np.average(sum(changes, []))}) # 結果をDataFrameに変換 d = pd.DataFrame(results) d['timeperiod'] = timeperiod # 変化率を%にし基準を0とする d['average'] = d['average'] * 100 - 100 df_upper = pd.concat([df_upper, d]) # RSI下限 results = [] for lower_limit in np.arange(10, 40): # RSIが下限を下回った日付を探す dates = df[(df['RSI'].shift() > lower_limit) & (df['RSI'] <= lower_limit)].index # RSIが下限を下回った日を基準に株価の変化率を求める changes = calc_change(df, dates) # 結果に追加 results.append( {'limit': lower_limit, 'average': np.average(sum(changes, []))}) # 結果をDataFrameに変換 d = pd.DataFrame(results) d['timeperiod'] = timeperiod # 変化率を%にし基準を0とする d['average'] = d['average'] * 100 - 100 df_lower = pd.concat([df_lower, d]) 結果をグラフにしてみます。 fig, axes = plt.subplots(nrows=2, figsize=(8, 8)) for timeperiod in df_upper['timeperiod'].unique(): d = df_upper[df_upper['timeperiod'] == timeperiod] axes[0].plot(d['limit'], d['average'], label=timeperiod) for timeperiod in df_lower['timeperiod'].unique(): d = df_lower[df_lower['timeperiod'] == timeperiod] axes[1].plot(d['limit'], d['average'], label=timeperiod) for i in range(2): axes[i].axhline(y=0, color='red', linestyle=':') axes[i].grid() axes[i].set_xlabel('RSI閾値') axes[i].set_ylabel('変化率(%)') fig.suptitle('RSI閾値と変化率') plt.subplots_adjust(top=0.9, bottom=0.1, hspace=0.3) axes[0].set_title('RSI上限値を超えた後') axes[1].set_title('RSI下限値を下回った後') plt.show() 計算期間ごとに線を引いています。RSI上限値のグラフでは70%手前から株価が下落する傾向が現れています。RSI下限値のグラフでは25%より少し先まで株価が上昇している傾向があります。ただこのグラフはわかりにくいですね。ヒートマップで作成し直してみます。 fig, axes = plt.subplots(nrows=2, figsize=(8, 16)) sns.heatmap(df_upper.pivot('timeperiod', 'limit', 'average'), square=True, cmap='Blues_r', ax=axes[0]) sns.heatmap(df_lower.pivot('timeperiod', 'limit', 'average'), square=True, cmap='Reds', ax=axes[1]) for i in range(2): axes[i].invert_yaxis() axes[i].set_xlabel('RSI閾値') axes[i].set_ylabel('RSI期間(日)') axes[0].set_title('RSI上限閾値と株価変化率') axes[1].set_title('RSI下限閾値と株価変化率') plt.subplots_adjust(top=0.9, bottom=0.1, hspace=0.3) plt.show() RSI上限値のグラフでは濃い青色の部分が株価が下落している部分です。このグラフだとRSI上限値が73,74%あたりとなります。このときのRSI計算期間は34〜40日となっています。RSI下限値のグラフでは濃い赤色の部分が株価が上昇している部分です。RSI下限値としては10〜24%、RSI集計期間は11〜38日の範囲で線状となっています。かなりはっきりした傾向が見られます。 Backtestingを使ってリターンを計算する 続いてBacktestingを使用してどれくらいのリターンが得られるか計算してみます。 Strategyクラスは下記のようにしました。RSIが上限値を超えたときに売り、RSIが下限値を下回ったときに買い、の単純な内容です。 class RsiStrategy(Strategy): ''' RSI Strategy ''' timeperiod = 14 rsi_upper = 80 rsi_lower = 20 def init(self): close = self.data['Adj Close'] self.rsi = self.I(talib.RSI, close, self.timeperiod) def next(self): ''' RSIが上限値を超えたら買われすぎと判断し売る RSIが下限値を下回ったら売られすぎと判断し買う ''' if crossover(self.rsi_lower, self.rsi): self.buy() elif crossover(self.rsi, self.rsi_upper): self.sell() シミュレートは下記で実行できます。 df = get_stock(TICKER, START_DATE, END_DATE) bt = Backtest( df, RsiStrategy, cash=INIT_CASH, trade_on_close=False, exclusive_orders=True ) 結果は25%のマイナスでした。 Start 2011-01-04 00:00:00 End 2020-12-30 00:00:00 Duration 3648 days 00:00:00 Exposure Time [%] 97.9992 Equity Final [$] 741527 Equity Peak [$] 1.24246e+06 Return [%] -25.8473 Buy & Hold Return [%] 163.934 Return (Ann.) [%] -3.03027 Volatility (Ann.) [%] 21.3582 Sharpe Ratio 0 Sortino Ratio 0 Calmar Ratio 0 Max. Drawdown [%] -70.8682 Avg. Drawdown [%] -9.02723 Max. Drawdown Duration 3131 days 00:00:00 Avg. Drawdown Duration 297 days 00:00:00 # Trades 10 Win Rate [%] 40 Best Trade [%] 42.2201 Worst Trade [%] -42.0814 Avg. Trade [%] -2.95518 Max. Trade Duration 876 days 00:00:00 Avg. Trade Duration 358 days 00:00:00 Profit Factor 0.93255 Expectancy [%] -0.567648 SQN -0.489027 _strategy RsiStrategy _equity_curve ... _trades Size EntryBa... dtype: object optimizeを使用して最適解を見つけます。（私のPCで2時間かかりました。） df = get_stock(TICKER, START_DATE, END_DATE) bt = Backtest( df, RsiStrategy, cash=INIT_CASH, trade_on_close=False, exclusive_orders=True ) stats, heatmap = bt.optimize( timeperiod=range(5, 41), rsi_upper=range(60, 91), rsi_lower=range(10, 41), return_heatmap=True) 最適時のリターンは214%です。(資産が3倍になったということですね。） このときのRSI計算期間は19日、RSI上限値が90%、RSI下限値が21%となっています。 Start 2011-01-04 00:00:00 End 2020-12-30 00:00:00 Duration 3648 days 00:00:00 Exposure Time [%] 97.9992 Equity Final [$] 3.14228e+06 Equity Peak [$] 3.14331e+06 Return [%] 214.228 Buy & Hold Return [%] 163.934 Return (Ann.) [%] 12.5035 Volatility (Ann.) [%] 23.5335 Sharpe Ratio 0.531307 Sortino Ratio 0.865442 Calmar Ratio 0.393282 Max. Drawdown [%] -31.7927 Avg. Drawdown [%] -4.09891 Max. Drawdown Duration 840 days 00:00:00 Avg. Drawdown Duration 56 days 00:00:00 # Trades 2 Win Rate [%] 100 Best Trade [%] 107.403 Worst Trade [%] 51.5613 Avg. Trade [%] 77.2973 Max. Trade Duration 3285 days 00:00:00 Avg. Trade Duration 1789 days 00:00:00 Profit Factor NaN Expectancy [%] 79.4823 SQN 462.381 _strategy RsiStrategy(time... _equity_curve ... _trades Size EntryBa... dtype: object 下記でトレード結果のグラフを表示します。 ほとんどトレードをしてないようです。ただ、株価が下がったときに一気に買っています。これって新型コロナで株価が暴落したときですね。 bt.plot() 下記でヒートマップのグラフも表示できます。 plot_heatmaps(heatmap, agg='mean', plot_width=2048, filename='heatmap') 変数が3つあるので、それぞれの組み合わせのヒートマップとなってます。 これでも傾向はわかりますが、3D散布図を作成してみることにします。 from mpl_toolkits.mplot3d import Axes3D d = heatmap.reset_index() fig = plt.figure() ax = Axes3D(fig) ax.scatter(d['timeperiod'], d['rsi_upper'], d['rsi_lower'], c=d['SQN']) ax.set_xlabel('RSI期間(日)') ax.set_ylabel('RSI上限') ax.set_zlabel('RSI下限') plt.show() 図ではイメージと逆ですが、薄い丸の方がSQNが高いことを表しています。（SQNが高い=良いということです。） 最後に 今回はRSIを使用した取引をシミュレートしてみました。前回のSMAではリターンが140%でしたが、RSIを使用した場合は214%となりました。なかなかの成果です。それにしても新型コロナが買いのチャンスだったとは… ソースはGitHubに置いています。

今回実施したこと 仮想環境上にDjangoをインストールした。 学んだ知識 仮想環境の構築方法とDjangoを用いたプロジェクトの作成。 PowerShellによる仮想環境の作成 ① WindowでPowerShelを開き、管理者による実行を押下する。 ②作業用のディレクトリをエクスプローラーから作成しておく。 ※今回はProgramingというディレクトリを作成しておく。 ③作成したディレクトリにPowerShell上で移動する。 cd C:\Users\Programing ④以下のコマンドを実行し、仮想環境を構築・立ち上げを実施する。 python -m venv test test\Scripts\activate 仮想環境に入れると先頭に(環境名)で表示される。 ※dectivateを入力することで仮想環境から抜けることができる。 djangoのインストール ①以下のコマンドを実行しdjangoをインストールする。 ※コマンドに失敗した場合はpipをUpdateする必要があることが考えられる。 pip install django ②仮想環境にプロジェクトを構築する。 ※今回は「mysite」というプロジェクトを構築する。 python test\Scripts\django-admin.py startproject mysite 悩んだポイントでアプリケーションとプロジェクトの概念が存在するが、 1つのプロジェクト内に複数のアプリケーションを構築することができる認識で問題なし。 ※アプリケーションの構築は次の記事に記載することにする。 ・アプリケーションとは、実際の処理を行うWeb アプリケーションを指します。 　　ex.)ブログシステムや公開レコードのデータベース、単純な投票アプリ ・プロジェクトとは、あるウェブサイト向けに設定とアプリケーションを集めたものです。 一つのプロジェクトには複数のアプリケーションを入れられる。 イメージ図 ここまでのディレクトリ構成は以下の通り。 mysite/ ├── manage.py | └──mysite(プロジェクトと同じディレクトリ名) | ├── __init__.py ├── settings.py ├── urls.py └── wsgi.py ③サーバを立ち上げる。 cd .\mysite\ $ python manage.py runserver ④ブラウザに「http://127.0.0.1:8000 」を入力しアクセスする。

TouchDesugnerで外部モジュールがimportできない!!!! ※この記事はメモ感覚で書いているので少々文が荒いところがあるかもしれません。 問題の概要 TouchDesignerにおいてposenetをimportしようとして、Python 64-bit Module Pathにそのモジュールが含まれるフォルダを指定しているがmoduleが見つからないと怒られる。 今回の解決方法 TouchDesignerの再起動!! Python 64-bit Module PathはTouchDesignerを再起動しないと更新されないそうです。 これはひどい これであなたの不具合も解消されることを祈ります... 参考 TouchDesignerで外部pythonのライブラリを参照する時には3.5系を指定する この記事のおかげで治りました!ありがとうございます!

Line botを作成するまでの手順をまとめる。 botはpythonで実装した。 LINEの開発概要は以下の解説を参照。 https://developers.line.biz/ja/docs/line-developers-console/overview/#page-title Line Developersのアカウントとチャネル作成 以下の手順に従い、Line Developersのアカウントとプロバイダー、チャネルを作成する。 プロバイダーの作成手順は省略。 チャネルは「Messageing API」を選択する。 チャネル作成時の項目は以下のような感じ。 チャネル名は適当に。 業種は個人のお試しでやっているので、大業種：個人、小業種：個人(IT・コンピュータ) にした。 作成したMessaging APIチャネルの「Messaging API設定」タブのQRコードをLINEアプリで「友だち追加」から読み込むと、LINEに友達追加できる。 応答メッセージの設定(Line Developers) 簡単なものであれば、Line Developers上で応答メッセージの設定ができる。「Messaging API設定」タブの応答メッセージの「編集」から編集する。 特定のメッセージに対して応答したり、時間帯別での応答も設定できる。 Webhookを利用した応答メッセージ LINEに投稿されたメッセージを外部サーバに転送し、外部サーバで作成した応答をLINEに応答することができる。 外部サーバはHerokuを使用する。Herokuの対応言語はここに記載されてる。 Herokuでアプリ作成 Herokuでアカウントを作成し、ログインする。 Herokuで 「New」→「Create new app」を選択。 「Create app」を押してApp作成。 App nameはHeroku全体(ログインアカウント内ではない)で一意な名前である必要があるため、既に使われている名前は指定できない。名前にこだわりがなければApp name入力なしで「Create app」を押すと、適当な名前のAppをHerokuが作成してくれる。 Botアプリ実装 PythonでBotアプリを実装する。 なお、LINEでは、Botアプリ向けにLine bot sdkを提供している。公式のSDKとサードパーティ製SDKがある。Pythonは公式SDKが提供されている。 他に、Webサーバとしてflaskを使用する。 $ mkdir line-bot-python $ cd line-bot-python $ python --version Python 3.6.4 $ python -m venv line-bot $ source line-bot/bin/activate $ echo line-bot/ > .gitignore $ pip install flask line-bot-sdk $ pip freeze > requirements.txt $ echo python-3.6.13 > runtime.txt $ echo web: python main.py > Procfile runtime.txtの仕様はここを参照。 Procfileの仕様はここを参照。 続いて、Botプログラムを作成。 main.py from flask import Flask, request, abort from linebot import ( LineBotApi, WebhookHandler ) from linebot.exceptions import ( InvalidSignatureError ) from linebot.models import ( MessageEvent, TextMessage, TextSendMessage, ) import os from argparse import ArgumentParser app = Flask(__name__) # 環境変数取得 # LINE Developersで設定されているチャネルアクセストークンとチャネルシークレットを設定 YOUR_CHANNEL_ACCESS_TOKEN = os.environ["CHANNEL_ACCESS_TOKEN"] YOUR_CHANNEL_SECRET = os.environ["CHANNEL_SECRET"] line_bot_api = LineBotApi(YOUR_CHANNEL_ACCESS_TOKEN) handler = WebhookHandler(YOUR_CHANNEL_SECRET) @app.route("/callback", methods=['POST']) def callback(): """ Webhookからのリクエストの正当性をチェックし、ハンドラに応答処理を移譲する """ # リクエストヘッダーから署名検証のための値を取得します。 signature = request.headers['X-Line-Signature'] # リクエストボディを取得します。 body = request.get_data(as_text=True) app.logger.info("Request body: " + body) # handle webhook body try: handler.handle(body, signature) # 署名検証で失敗した場合、例外を出す。 except InvalidSignatureError: app.logger.warn("Invalid Signature.") abort(400) # handleの処理を終えればOK return 'OK' @handler.add(MessageEvent, message=TextMessage) def handle_message(event): """ LINEへのテキストメッセージに対して応答を返す Parameters ---------- event: MessageEvent LINEに送信されたメッセージイベント """ line_bot_api.reply_message( event.reply_token, TextSendMessage(text="応答です。 " + event.message.text)) if __name__ == "__main__": # app.run() arg_parser = ArgumentParser( usage='Usage: python ' + __file__ + ' [--port <port>] [--help]' ) # Herokuは環境変数PORTのポートで起動したWeb Appの起動を待ち受けるため、そのポート番号でApp起動する arg_parser.add_argument('-p', '--port', type=int, default=int(os.environ.get('PORT', 8000)), help='port') arg_parser.add_argument('-d', '--debug', default=False, help='debug') arg_parser.add_argument('--host', default='0.0.0.0', help='host') options = arg_parser.parse_args() app.run(debug=options.debug, host=options.host, port=options.port) Heroku CLIインストール ここの内容に従い、Heroku CLIをインストールする。 Macでのインストールは以下コマンドで。 $ brew tap heroku/brew && brew install heroku チャネルアクセストークンの設定 チャネルアクセストークンをHerokuの環境変数に設定する。 Line DevelopersのMessaging API設定の、チャネルアクセストークン（長期）の「発行」を押すと、トークンが発行されるので、これをコピーする。 コマンドラインより、以下を実行する。 $ heroku config:set CHANNEL_ACCESS_TOKEN=[Line Developersで発行したチャネルアクセストークン] -app=[HerokuでCreate appした名前] チャネルシークレットの設定 チャネルシークレットをHerokuの環境変数に設定する。 Line Developersのチャネル基本設定の、チャネルシークレットをコピーする。 コマンドラインより、以下を実行する。 $ heroku config:set CHANNEL_SECRET=[Line Developersで発行したチャネルシークレット] -app=[HerokuでCreate appした名前] Herokuにアップロード コマンドラインからプログラムをアップロードする。 $ heroku login $ git init $ heroku git:remote -a [HerokuでCreate appした名前] $ git add . $ git commit -am "line bot app" $ git push heroku master Webhook設定 Line DevelopersよりWebhookを設定する。 Messaging API設定より、Webhook URLの「編集」を押す。 Webhook URLを設定する画面になるので、herokuのURLを指定する。 URLを設定したら、「検証」を押して、Webhookが動作するかを確認する。 動作する場合は以下のように表示される。 以下のようなエラーが表示された場合は、コマンドラインからheloku logsコマンドを実行するなどして、エラー内容を確認して対処する。 確認後、Webhookの利用のスイッチを有効にする。 設定後、LINEからメッセージ送信すると、Webhookで設定した応答もLINEに返されるようになる。 上記は、LINEの応答メッセージとWebhookが両方返されている。 Webhookのみを有効にする場合は、Line Developersの設定で応答メッセージをオフにする。 参考

はじめに PythonでNumPyによる配列の操作方法で、比較と論理演算の使い方をまとめました。配列の前要素について比較・論理演算子が適用されるので、コーディングが簡易になり非常に便利だと思います。 動作環境 OS:MacOS 11.2.3 (Big Sur) Python:3.9.0 NumPy:1.20.2 比較演算子を利用した例 sample01 >>> import numpy as np >>> a = np.arange(15) >>> h = np.reisize(a, (5,3)) >>> h array([[ 0, 1, 2], [ 3, 4, 5], [ 6, 7, 8], [ 9, 10, 11], [12, 13, 14]]) >>> h > 8 array([[False, False, False], [False, False, False], [False, False, False], [ True, True, True], [ True, True, True]]) 簡単な例ですが、”h > 8"とすることで、要素のうち8より大きいものについてはTrueで、それ以外ではFalseであるブール値の配列が返されます。 次の例のようにインデックス参照によるデータ抽出方法もあります。戻り値は１次配列となります。 sample02 >>> h[h > 8] array([ 9, 10, 11, 12, 13, 14]) >>> h[h > 8].ndim 1 >>> h[h > 8].shape (6,) >>> h[(h>4)|(h>12)] >>> h[(h<4)|(h>12)] >>> array([ 0, 1, 2, 3, 13, 14]) >>> h[(4<=h)&(h<=12)] >>> array([ 4, 5, 6, 7, 8, 9, 10, 11, 12]) where()関数による対象要素の抽出 numpyにはwhere()関数が用意されており、SQL文のように要素の中から、比較・論理演算子によって必要なデータを抽出することができます。また、TrueかFalse以外にも出力する値を設定することができます。 sample03 >>> np.where(h > 8, 1, 0) array([[0, 0, 0], [0, 0, 0], [0, 0, 0], [1, 1, 1], [1, 1, 1]]) >>> np.where(h % 2 == 0, h, 0) array([[ 0, 0, 2], [ 0, 4, 0], [ 6, 0, 8], [ 0, 10, 0], [12, 0, 14]])

はじめに PythonのNumPyで多次元配列を操作する方法をまとめました。 別の記事で１次元配列の操作方法についてまとめましたので、こちらも併せてご覧いただければと思います。 動作環境 OS: MacOS 11.2.3 (Big Sur) Python: 3.9.0 NumPy: 1.20.2 多次元配列の基本的な操作方法 多次元配列の生成 以下のように、多次元のlist配列を引数に指定することで生成できます。 sample01 >>> import numpy as np >>> a = np.arange(8, dtype=float) >>> a array([0., 1., 2., 3., 4., 5., 6., 7.]) >>> b = np.array([a, a * 2]) #引数は[]カッコで囲むことに注意。 >>> b array([[ 0., 1., 2., 3., 4., 5., 6., 7.], [ 0., 2., 4., 6., 8., 10., 12., 14.]]) 多次元配列になってもスライス機能は使うことができます。 sample02 >>> b[:,:4] array([[0., 1., 2., 3.], [0., 2., 4., 6.]]) >>> b[1,2:5] array([4., 6., 8.]) zeros, ones, emptyによる初期化 zeros（全ての要素が0）、ones（全ての要素が1）、empty（全ての要素が空）によって多次元配列の初期化を行うことができます。 sample03 >>> import numpy as np >>> a = np.zeros((2,3),dtype=int) >>> a array([[0, 0, 0], [0, 0, 0]]) >>> b = np.ones((2,3),dtype=float) array([[1., 1., 1.], [1., 1., 1.]]) >>> c = np.eyes(3) #対角線が全て1となる正方行列の生成 array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]]) 多次元配列のメタ情報 sample04 >>> import numpy as np >>> a = np.zeros((2,3,4),dtype='i') >>> a.size #全要素数 24 >>> a.ndim #次元 3 >>> a.shape #行列の型 (2, 3, 4) >>> a.dtype #要素のデータ型 dtype('int32') 多次元配列で利用できる数理関数 １次元配列と同様に、総和(sum)、平均(mean)、分散(var)、標準偏差(std)等が多次元配列でも利用できます。 1次元配列と違い、多次元配列の場合は計算する次元を指定できます。わかりやすく説明するために２次元配列で例を示します。 sample05 >>> import numpy as np >>> a = np.arange(0, 8, 2) >>> b = np.array([a, 2*a]) >>> b array([[ 0, 2, 4, 6], [ 0, 4, 8, 12]]) >>> b.sum() #全要素の総和 36 >>> b.sum(axis=0) #列ごとに要素の総和を算出 array([ 0, 6, 12, 18]) >>> b.sum(axis=1) #行ごとに要素の総和を算出 array([12, 24]) 上記のように要素の計算では、次元（軸）を指定することができます。 形とサイズの変更 reshape()関数、resize()関数 reshape()関数やresize()関数によって、もとの多次元配列の次元や型を変更することができます。 sample06 >>> import numpy as np >>> a = np.arange(15) >>> a array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]) >>> a.reshape((3,5)) array([[ 0, 1, 2, 3, 4], [ 5, 6, 7, 8, 9], [10, 11, 12, 13, 14]]) >>> a.reshape((5,3)) array([[ 0, 1, 2], [ 3, 4, 5], [ 6, 7, 8], [ 9, 10, 11], [12, 13, 14]]) >>> np.resize(a, (3,5)) array([[ 0, 1, 2, 3, 4], [ 5, 6, 7, 8, 9], [10, 11, 12, 13, 14]]) >>> np.resize(a, (5,3)) array([[ 0, 1, 2], [ 3, 4, 5], [ 6, 7, 8], [ 9, 10, 11], [12, 13, 14]]) >>> スタッキング(hstack()、vstack()) 同じ型の行列を水平にスタックしたり、垂直にスタックすることができます。具体的に例を見ていきます。 sample07 >>> import numpy as np >>> a = np.arange(15) >>> b = np.resize(a, (3, 5)) >>> b array([[ 0, 1, 2, 3, 4], [ 5, 6, 7, 8, 9], [10, 11, 12, 13, 14]]) >>> c = np.hstack((b, 2 * b)) >>> c array([[ 0, 1, 2, 3, 4, 0, 2, 4, 6, 8], [ 5, 6, 7, 8, 9, 10, 12, 14, 16, 18], [10, 11, 12, 13, 14, 20, 22, 24, 26, 28]]) >>> d = np.vstack((b, 2 * b)) array([[ 0, 1, 2, 3, 4], [ 5, 6, 7, 8, 9], [10, 11, 12, 13, 14], [ 0, 2, 4, 6, 8], [10, 12, 14, 16, 18], [20, 22, 24, 26, 28]]) フラット化 今までは１次元の配列を多次元配列に型変換する例を見てきましたが、こちらでは多次元の配列を１次元配列にする（フラット化）する例を見ていきます。 sample08 >>> import numpy as np >>> a = np.arange(15) >>> b = a.reshape((5, 3)) >>> b array([[ 0, 1, 2], [ 3, 4, 5], [ 6, 7, 8], [ 9, 10, 11], [12, 13, 14]]) >>> b.flatten() array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]) >>> b.flatten(order='C') array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]) >>> b.flatten(order='F') array([ 0, 3, 6, 9, 12, 1, 4, 7, 10, 13, 2, 5, 8, 11, 14])

概要 今回はPythonで出るエラーについてまとめていきたいと思います。 初学者の時は、エラーが出ると何から手をつければ良いかわからなくなる時はよくあると思います。 ググろうにもググり方がわからないという状況になって、エラーを極端に怖がってしまいます。 しかし、本来エラーはありがたいものなのです。 なぜかというと、どこかダメなのかを丁寧に教えてくれるので、慣れてくると良い友達のような感覚になります。 なのでこの記事を読んだ後はエラーを友達のように考えるようにしましょう！ それでは早速見ていきましょう！ エラーが出た時の対策 まずはエラーが出た時の最初の対策を見ていきましょう。 まずは一番最後の行に出力されている「〜Error」の部分をみましょう。 次にそのエラーがどんなものかを理解して、「〜Error」の後ろ部分をみてみましょう。 そこをみてもすぐに解決しない場合は、「〜Error」の後ろ部分をググってみましょう。（ここで気をつけるのが、パスや個人情報部分を取り除いてググってみてください。） それでも解決しない時は、エラーを上に戻って、下にスクロールしながらエラーをみてみましょう。これでどこでエラーが出ているのか、どのようなエラーが出ているのかの細かい部分を見ることができます。 これでダメだったら質問投稿サイトなどで質問してみましょう！ エラー一覧 エラーの対処法を確認したところで、早速エラーの種類の一覧を見ていきましょう！ SyntaxError 「構文エラー」です。 「Pythonの構文に沿ってないよ！」怒られている状態です。 print"かるでね") ^ SyntaxError: invalid syntax この場合は「(」が足りていませんね。 矢印の部分を確認するとわかりやすいですが、今回はちょっとずれていますね…。 よくある「SyntaxError」の種類を紹介します。 クオーテーションのつけ忘れ。 「[]」や「()」の両方または、片方のつけ忘れ。 「=」と「==」の使い間違い。 NameError 「名前」に対してのエラーです。 「こんな名前どこにも定義されていないよ！」と怒られています。 ----> 1 print(かるでね) NameError: name 'かるでね' is not defined 「かるでね」という変数を定義していないのに、「かるでね」という変数を出力しようとしてしているために起こるエラーです。 エラーの部分の名前を定義するか、文字列にするなど、自分の状況に合わせて対処してください。 IndentationError 「インデント」に対するエラーです。 「変なところにインデント入ってるよ？」というエラーです。 インデントがずれているために起こるエラーです。 def prints(): print("かるでね1") print("かるでね2") print("かるでね3") print("かるでね2") ^ IndentationError: unexpected indent 「かるでね2」のインデントが１つずれているのが確認できますね。 単純にインデントを整えればエラーがなくなります。 TypeError 「型」に対してのエラーです。 「いやいやその型でそれはできんよ。」と怒られています。 cardene = 1234 for i in cardene: print(i) 1 cardene = 1234 ----> 2 for i in cardene: 3 print(i) TypeError: 'int' object is not iterable 「int型（数字型）」だから「for文」では回せないよと言われています。 「cardene」変数を文字列にするなどすれば問題なく出力できます。 AttributeError 「属性」に関するエラーです。 「そんな属性ないよ？」と注意されています。 import math print(math.PI) ----> 3 print(math.PI) AttributeError: module 'math' has no attribute 'PI' 「math.pi」というものはありますが、「math.PI」というものが存在しないために起こるエラーです。 単純にスペルミスを無くせば治せます。 cardene = "かるでね" cardene.append(100) 1 cardene = "かるでね" ----> 2 cardene.append(100) AttributeError: 'str' object has no attribute 'append' ModuleNotFoundError 「モジュール」に対するエラーです。 「そんなモジュールなくね？」と怒られています。 自分の環境に「pip」や「conda」でインストールできていないのにインポートしていたり、モジュールの名前を間違えたときに起こるエラーです。 import cardene ----> 1 import cardene ModuleNotFoundError: No module named 'cardene' 「cardene」というモジュールは存在しないので、適切なモジュールをインストールすれば解決します。 cardene = "かるでね" cardene.append(100) ----> 2 cardene.append(100) AttributeError: 'str' object has no attribute 'append' 「cardene」という文字列に対して、「100」を付け足そうとしているのでエラーが起きています。 「cardene」をリストにすれば解決できます。 ValueError 「型」に対しての「値」に関するエラーです。 「その値じゃ受け取れないよ〜」と怒られています。 cardene = "かるでね" print(int(cardene)) ----> 2 print(int(cardene)) ValueError: invalid literal for int() with base 10: 'かるでね' 「cardene」は「かるでね」という文字が入っているので、数字に変換することはできないために起きているエラーです。 「かるでね」を「100」などに変更すれば解決します。 ZeroDivisionError 「ゼロ」で割り算などをしたときに出るエラーです。 「ゼロで割っちゃいかんでしょ」と怒られています。 print(100 / 0) ----> 1 print(100 / 0) ZeroDivisionError: division by zero 「ゼロ」で割ってしまっているので、別の数字で割れば解決できます。 ImportError こちらも「モジュール」に対するエラーです。 「ファイル見てみたけど、そんなモジュールないよ？」と怒られています。 「from ~ import ~」で２つ目の「~」の部分が存在しないときに起こるエラーです。 from tensorflow import cardene ImportError: cannot import name 'cardene' from 'tensorflow' 「tensorflow」というライブラリは存在しますが、「tensorflow」の中には「cardene」というモジュールが存在しないのでエラーが出ています。 「import」部分を存在する「モジュール」に変更すれば解決します。 KeyError 「辞書」に指定の「キー」がないために起こるエラーです。 「そんなキーは存在しないよ？」と怒られています。 cardene = {"cardene1": "cardene1", "cardene2": "cardene2", "cardene3": "cardene3"} cardene["かるでね"] ----> 2 cardene["かるでね"] KeyError: 'かるでね' 「cardene」という辞書には「かるでね」という「キー」が存在しないのが確認できます。 「“cardene1”」のように指定すれば解決します。 IndexError 「リスト」などの「配列」で間違った「インデックス番号」を指定したときに起こるエラーです。 「そんなインデックス番号はないよ？」と怒られています。 cardene = [1, 2, 3, 4, 5] print(cardene[10]) ----> 2 print(cardene[10]) IndexError: list index out of range 「cardene」という「リスト」は要素数が「５」なのに、「インデックス番号10番」にアクセスしようとしているのが確認できます。 「10」を「3」などに指定し直せば解決します。 FileNotFoundError 「ファイル」が存在しないときに出るエラーです。 「君、存在しないファイルにアクセスしようとしてるよ」と言われています。 with open('file.txt') as f: print(f.read()) ----> 1 with open('file.txt') as f: 2 print(f.read()) FileNotFoundError: [Errno 2] No such file or directory: 'file.txt' 「file.txt」というファイルを開こうとしていますが、そんなファイルはないよと言われているのが確認できます。 「パス」や「ファイル名」を修正することで解決できます。 FileExistsError すでに存在しているファイルを作成しようとするときに起こるエラーです。 「もうそのファイル名は存在しているから作れんよ」と怒られています。 import os os.mkdir('cardene') ----> 3 os.mkdir('cardene') FileExistsError: [Errno 17] File exists: 'cardene' 今回は「ディレクトリ」ですが、すでに存在しているよと言われているのが確認できます。 存在しない名前に変えたり、存在するファイルやディレクトリを削除することで解決できます。 最後に 今回はPythonで出るエラーについてまとめてきました。 結構色々な種類があって最初は混乱してしまいますが、何回も遭遇すると勝手に覚えてきます。 まずはどんなエラーがあるのかをこの記事で確認してもらって、エラーとどんどん仲良くなっていってください。 もしお友達になれるようなことがあれば報告してもらえると嬉しいです！ 何か質問などがあれば下のTwitterのDMを送るか、ここのコメントで聞いてください！ それでは！

本書は筆者たちが勉強した際のメモを、後に学習する方の一助となるようにまとめたものです。誤りや不足、加筆修正すべきところがありましたらぜひご指摘ください。継続してブラッシュアップしていきます。 © 2021 NPO法人AI開発推進協会 本書はTransformerを使った物体検出モデルであるDETRを説明します。（Transformerを理解している前提で記載しています。まだ理解していない方は当法人で作成している別冊「Seq2Seq&Transformer」を先に読んでください。） 【参考文献、サイト】 論文　End-to-End Object Detection withTransformers End-to-End Object Detection with Transformers(DETR)の解説 Transformerを物体検出に採用！話題のDETRを詳細解説！ １．はじめに DETRはFacebook AI Research（FAIR）が2020年5月に公開した、Transformerを使った初めての物体検出モデルです。これまでの物体検出は、NMS（Non Maxinum Suppression）、Anchor boxの数、アスペクト比、Bounding boxのデフォルト座標等を人手でチューニングする必要があり実質的にEnd to Endが実現されていませんでした。 ※物体検出モデルの詳細は別冊のSSDまたはFaster R-CNNを参照してください。 DETRは、物体検出問題を直接集合予測問題（a direct set prediction problem）として捉え、NMSやRPN等の複雑な仕組みや上記人手のチューニングをなくした「CNN ＋ Transformer」のシンプルな構成でEnd to Endな物体検出を実現します。 ※（筆者のイメージでは）DETRでのTransformerのシーケンス（時系列）は「画像内のオブジェクトの関連性」を把握するために使われていると思われます。例えば、馬がいればその上に乗っているオブジェクトは人間である可能性が高い等です。 DETRはTansformerを活用することでシンプルな構成ながらFaster R-CNN等に匹敵する性能を達成しています。 ２．DETRアーキテクチャ 従来の物体検出モデルは、複数の物体を検出してNMSで絞り込むというプロセスを介していましたが、DETRでは入力した画像内に写っているすべてのオブジェクトを一括して推論します。また、画像内の推論したすべてのオブジェクトとすべての正解ラベルとをハンガリアン法（後述、推論と正解を比べ最も当てはまりの良いペアを見つけてそれらに対して損失をとるものです）で対応付けを実施しています。 (1) DETRの推論フロー 推論フローは図１のように非常にシンプルです。 図１　DETR推論時のフロー 画像をCNN（ResNet等）に入力する。 CNNで出力されたチャネル数（ex.2048）を１×１畳み込みでチャネル数をより小さい次元d(ex.256)に縮小する。 特徴マップが生成される（d個） Transformerにd個の特徴マップを入力する N個のオブジェクト情報（クラス、位置、サイズ）が出力される※N個が一括して推論されます。なお、Nは１００個等入力画像に含まれるよりも十分に多い数を選びます (2) ハンガリアン法（推論結果と正解ラベルの自動対応付け） 直接集合予測問題を解く上で、推論と正解とのマッチングが正しく行われているのか正しく判断できる必要があります（人手でなく機械で）。DETRでは、二部マッチング問題として考え、ハンガリアン法で適切な対応付けを行います。 二部マッチング問題 異なる２つの集合の要素をどのように組み合わせると目的を最大値で達成することができるの かを解くものです。 ハンガリアン法 二部マッチング問題を効率的に解くための手法。全検索すると n! の計算量が必要になるとこを　n3 で計算できます。　DETRではこの手法を用いてオブジェクトの推論結果（N、 ex.100）と正解ラベル（Nとの差分はすべてNo Object）の対応付けを行います。 図２　２部マッチングの例 　 実際の対応付けイメージは下記にようになります。 まず正解ラベル y と推論結果 y^ を突き合わせて、コスト行列（ex. N=100であれば、100×100行列に各ペアのロスが入る）を作ります。 ハンガリアンアルゴリズムでこの行列を解いて、総和が最小なマッチングパタンσ^（重複なしのN（100）ペア）を求めます。 このマッチングパターン σ^ が正解ラベル y (n個はオブジェクト、「N-n」個はNo Object)と推論結果 y^ の対応付けになります。 (3) ハンガリアンアルゴリズムを用いたロスの算出 ロス計算は以下の流れになります。 DETRに画像を入力し、推論結果y^(N個分（ex.100個）のオブジェクト推論)を出力※yiには分類クラス情報ciと矩形情報biを持つ 正解ラベルyと推論結果y^をもとに、ハンガリアンアルゴリズムで最適なマッチングパターンσ^を見つけ出す 正解ラベルyと推論結果y^の最適なマッチングパターンσ^におけるロスを計算する 　 マッチングパターンσ^の計算は以下になります。 Lbox : Generalized IoUロス（IoUに距離の概念を追加したもの）と回帰ロス（位置・サイズ）を足したものです。回帰ロス（位置・サイズ）はスケール変化に弱いため、IoUロスを追加しています。 Generalized IoUロスと回帰ロス（位置・サイズ）の和を矩形ロス（Lbox）としています。 Lmatch : 予測と正解ラベルをマッチングしたときに生じるロスです。 クラスロス(そのクラスである確率 X - 1 )と、矩形ロス（ Lbox ）の和をLmatchとしています。 σ^ : Lmatchをすべて計算し、その結果から物体の対応関係一覧を得ます。 各マッチングパターン σ∈ €Nにおいて、各ペア（正解ラベル yi と推論結果 yi^ ）のロス Lmatch の総和が最小となるマッチングパターン σ^を求め ています LHungarian : σ^によって得た物体の対応関係一覧から最適な予測と正解ラベルのマッチングを算出します。 クラスロスと矩形ロス（Lbox）の和の全オブジェクト総和を LHungarian としています。 ３．DETRアーキテクチャ DETR全体のアーキテクチャは図３のようにシンプルで、「CNNバックボーン」、「エンコーダー・デコーダートランスフォーマ」、最終的な検出予測をシンプルな「FFN（フィードフォワードネットワーク）の３つのコンポーネントで構成されています。 図３　DETRアーキテクチャ (1) 処理の流れ CNN出力の特徴マップ（ResNetだと2048）を、 1 × 1 Convolutionでdチャネルに圧縮する Transformer encoderへ、サイズH * Wの特徴マップd個を、H * W個のd次元特徴として入力する。このときpositional encodingで位置情報を補足する。 Transformer encoderは、H * W個の中間オブジェクト特徴（各d次元）を出力する Transformer decoderへ、 H * W個の中間オブジェクト特徴とN個のobject query（※１）を入力する Transformer decoderは、N個のオブジェクト特徴（各d次元）を出力する N個のオブジェクト特徴をそれぞれＦＦＮを通して、Ｎ個のオブジェクト推論結果（クラス、矩形）が出てきます ※１　画像内から「物体検出」と「分類」をするために学習される値。初期値は完全にランダムなベクトル数値。分類するクラスに対して十分に大きな値にします。（イメージ的には画像内に映っている物体の関連性をみながら対応付けていくと思われます） 図４は原論文に記載されているＮ＝１００のときの内２０個のボックスの予測を可視化したものです。画像サイズごとにポイントは色分けされており、緑の色が小さい画像、赤は大きな横長のボックス、青は大きな縦長のボックスを示します。各Ｎそれぞれで見方の傾向に違いがあることが見て取れます。この違いが適切に物体検出と分類を行うことを可能にしています。 図４　COCO 2017 val画像セットでボックス予測の可視化 (2) Auxiliary decoding loss オブジェクトの識別性能をあげるために、補助（Auxiliary）ロスを用います。デコーダも各層の出力にFFN（最終層のコピー）を設けて、各層で最終層と同様の推論をします。そして、各層でロス計算を同様にして、各層で誤差伝搬を走らせます。 ４．詳細なアーキテクチャ DETRにおいて、各アテンション層に展開される位置エンコーディングを用いたトランスフォーマーの詳細は図５のとおりです。　CNNバックボーンからの画像特徴量は、トランスフォーマー・エンコーダーを通じて、空間位置エンコーダとともに渡されます。空間位置エンコーダ（図５のSpatial positional encoding）は、マルチヘッド・セルフアテンション層ごとにQueryとKeyに追加されます。 その後、デコーダはQuery(初期値はゼロに設定されています)、出力位置エンコーディング（オブジェクトクエリ）、エンコーダメモリを受信し、予測したクラスラベルとバウンディングボックスの最終的なセットを生成します。　この処理は、複数のマルチヘッドのセルフ・アテンションとデコーダエンコーダ・アテンションを介して行います。最終デコーダの最初のセルフアテンション層はスキップ可能です。 図５　DETRのトランスフォーマーの構成図 ５．識別結果（論文より） (1) Eecoder最終層のattentionマップ 図６は、Encoderの256スロットの中から、顕著なattentionマップを持つ4スロットを可視化した図です。これからEncoderの段階で既にオブジェクト分離できているのがわかります。 図６　　エンコーダのセルフアッテンションの可視化 (2) Decoder最終層のattentionマップ 図７は、Decoderの100スロットの中から、顕著なattentionマップを持つ2スロットを可視化したものです。上述のようにEncoderで大まかなオブジェクト分離は完了しているため、Decoder側ではオブジェクト境界（足や鼻先）を集中して見ることで、精密なオブジェクトの境界線を判断しています。 図７　デコーダアテンションの可視化 ６．シンプルさ 論文上に記載されているPytorchで記載された推論コードは下記のように約３０行程度のコードになっています。 図８　DETRコード ７．おわりに 本稿では、直接集合予測を行うためトランスフォーマーと二部マッチングを使い物体検出を行うDETRについて説明しました。DETRは実装が簡単でセグメンテーション領域への拡張も可能です。さらにFaster R-CNNよりもラージオブジェクトに対しては高い性能が出ます。一方でスモールオブジェクトの精度等では課題がありますが、今後改善がされていくものと思いますので継続して「注視」していきます。 物体検出（SSDやFaster-RNN）や、Transformerについても別冊で執筆しておりますので興味のある方は、ぜひ参照ください。

日本語化する方法 変更方法 設定を開く Application -> Advanced settings -> Langage おわりに 調べても全然わからなかったのでメモとして残しておく。

初めに この記事は、前回投稿した記事の続きである。 確認されることを推奨する。 wsl2上で、python3が実装されている前提で話を進める。 前回記事から初めてプログラムをする人にとっては、正直順序を飛ばしている感が非常にある。 人によっては、pythonのfor文やif文などの使い方から学んだほうが肌に合う人もいるだろう。 ただ、自分は、とりあえず動いて楽しいものを、コードのコピペでもなんでもいいから実装して、そのあとで、中身を理解していく、という順序のほうが肌に合っているので、print("hello world")の次に、"OpenJTalkを用いた音声合成"を持ってきた。 以下に、サンプルコードや環境構築のコマンドなどを載せているが、一回で全部理解する必要はない。できたおもちゃで遊んでいくうちになんとなく理解していけばいいだろう。 お詫び 前回記事から、時間が経ってしまった。 筆が進まなかった。 ゴメンネ。 目次 OpenJTalkについて OpenJTalkのインストール wslから、.wavファイルを再生できるようにするための作業 pythonを使った、OpenJTalkのhello world OpenJTalkについて OpenJTalkは、「入力された日本語テキストに基づいて自由な音声を生成するHMMテキスト音声合成システム」である。 公式サイトより抜粋。 平たく言うと、日本語の文字列をもとに、.wavファイルを作成してくれる、ありがたい道具である。 感謝しながら使いましょう。 .wavファイルとは、.mp3や.m4aと同じように、音声データを表現するファイルの拡張子の一種である。 OpenJTalkを使用するために、必要なものは、大きく以下の３つである。 ・OpenJTalk本体 ・辞書モデル ・音響モデル OpenJTalk本体: 公式サイトで配布されている。 windowsで使用するときは、ここからダウンロードしてくればよいが、今回は、wsk上にコマンドプロンプトでダウンロードするので、特にwebページに足を運ぶ必要はない。 辞書モデル： 今回は、一般に使用される日本語辞書の"MeCab"を使用する。 音響モデル： 音響モデルは、.htsvoiceファイルなどで表現される。 今回は、代表的な音響モデルの一つである"mei"を使用する。 なお、meiには、"angry"、"bashful"、"happy"、"normal"、"sad"、の感情モデルがそれぞれ存在する。 これらの扱いについても後程述べる。 また、音響モデルは、様々な有志が作成し、ネットで公開してくださっている。(このサイトなど) 利用規約を守って使用してみるとよい。 また、日本語の.wavファイルを生成するだけなら、"SofTalk"や"SHABERU_Light"などのアプリが配布されている。 非常に便利なアプリである。 ただ、これらは完成されたアプリで、自分のコードに基本的に組み込めないため、今回は特に触れない。 OpenJTalkのインストール 以下の作業はすべてwsl上での作業である。 注意せよ。 ただ、 以下のコマンドを、wsl上の好きなところで実行せよ。 ただ、前回も述べた通り、"workspace"などの名前のディレクトリを作成して、その中で行うほうがマナーが良い。 おまじない sudo apt update sudo apt upgrade 環境構築に使用する道具のインストール sudo apt install -y wget sudo apt install -y unzip wgwt : web上から、URL指定でファイルをダウンロードする道具 unzip : .zipファイルを解凍するための道具 OpenJTalk本体のインストール sudo apt install -y open-jtalk 実体としては、wsl上の、/usr/bin に、"open_jtalk"として存在する。 (環境によっては、実体の場所は変わるのかも。詳しくないので、有識者の方、いらっしゃいましたら教えてください。) なお、この"実体として存在する場所"は、実際にpythonで使用するときに必要になってくるので、よく覚えておこう。 辞書モデルのインストール sudo apt install -y open-jtalk-mecab-naist-jdic 実体としては、wsl上の、/var/lib/mecab/dic/open-jtalk/naist-jdic に存在する。 この場所は覚えておこう。 音響モデル実装装置(みたいなもの)のインストール sudo apt install -y hts-voice-nitech-jp-atr503-m001 これの場所は特に覚えておかなくてよい。 音響モデル(mei)本体の.zipファイルのダウンロード wget http://sourceforge.net/projects/mmdagent/files/MMDAgent_Example/MMDAgent_Example-1.7/MMDAgent_Example-1.7.zip 上記.zipファイルの解凍 sudo unzip MMDAgent_Example-1.7.zip この中の、MMDAgent_Example-1.7/Voice/meiが、"mei"音響モデルの全体であり、その下の階層に、 前述した、"angry"、"bashful"、"happy"、"normal"、"sad"の感情の音響モデル(.htsvoiceファイル)が存在する。 この場所(~/workspace/MMDAgent_Example-1.7/Voice/meなど)は、覚えておく必要がある。 ただ、この音響モデル("mei")は、よく使うので、どこからでもアクセスできる場所にコピーして、その場所を覚えたほうが、風情がある。 以下のコマンドを実行し、ゲットした"mei"を、別のところに移動させよう。 meiを扱いやすい場所に移動。 sudo cp -r MMDAgent_Example-1.7/Voice/mei/ /usr/share/hts-voice/ /usr/share という場所に、"hts-voice"という場所を作成し、そこに"mei"を保管。 この場所はちょっと特殊な場所なので、cpコマンドの前に、sudoをつける必要がある。 この、/usr/share/hts-voice/mei/の場所は、覚えておく必要がある。 以上で、OpenJTalkのインストールは終了。 wslから、.wavファイルを再生できるようにするための作業 ここでは、wsl側から、音を再生できるようにするための作業である。 おそらく、macを使用している場合は、ここでの作業はスキップしてよいと思う。 この章の作業は、windows側の作業と、wsl側の作業の二つがある。注意せよ。 音を再生させるには、ハードウェアが必要で、ハードウェアを管理しているのは、windowsなので、windows側でも作業する必要がある。 なお、ここでのコマンドなどはすべて理解しなくてもよい。難しいから。 作業は、この神サイトに沿って行う。 非常に分かりやすくて丁寧なサイトである。 これより丁寧に説明できる気がしないので、上記の記事を読んでm(_ _)m この章の作業は、正直めんどい。 そして、音声合成(.wavファイル生成)自体には、本質的に関係ない。 そのため、この章の作業を飛ばして、次の章へ進んでも、特に問題はないと思う。（この章を飛ばしても、ロボットなどの別のハードウェアにしゃべらせる場合は、問題ない。） この章を飛ばす場合は、生成した.wavファイルを、windows側から指定して、再生させる必要がある。 pythonを使った、OpenJTalkのhello world 詳しい話をする前に、以下のコマンドを、wslのマナーの良い場所(workspace上など)で、実行してみてほしい。 前章を飛ばした人は、window側からアクセスできるように、/mnt 以下の場所で実行してみてほしい。 以下のコマンドは、このサイトを参照して、作成した。 OpenJTalkを使用して、.wavファイルを生成するための、最低限のコマンド。 echo 'おはようございます。僕は元気です。' | open_jtalk -x /var/lib/mecab/dic/open-jtalk/naist-jdic -m /usr/share/hts-voice/mei/mei_happy.htsvoice -ow test.wav うまくいけば、カレントディレクトリ（今いる場所）に、"test.wav"という名前の、.wavファイルができていると思う。 これを再生すれば、meiがhappyな感情で"おはようございます。僕は元気です。"としゃべってくれる。 ただし、このような使い方は基本的にしない。だるいから。 上記のコマンド（シェルコマンドという）を、pythonを使用して、可読性と汎用性を高めつつまとめたのが、以下のコード。 ここなどにも類似コードが書かれている。まとめ方の参考にすると良い。 このpythonファイルを実行するためには、"subprocess"と"os.path"というライブラリが必要である。 一般に、このようなライブラリは、外部からインストールしなければ使えないが、これらは、"python標準ライブラリ"と飛ばれるもので、外部のモジュールのような下準備なしで、いきなり使用(import)できる。 make_wav.py import subprocess #シェルコマンドをpythonから命令するのに必要なもの import os.path # OpenJTalkをインストールしたパス OPENJTALK_BINPATH = '/usr/bin' #OpenJTalk本体の場所 OPENJTALK_DICPATH = '/var/lib/mecab/dic/open-jtalk/naist-jdic' #辞書モデルの場所 OPENJTALK_VOICEPATH = '/usr/share/hts-voice/mei/mei_{emotion}.htsvoice' #音響モデルの場所 def make_wav(text, speed=1.0, emotion='normal', output_file='__temp.wav', output_dir=os.getcwd()): """ OpenJTalkを使ってmeiの声で音声合成してwavファイルを作る関数。 引数emotionには'normal', 'happy', 'bashful', 'angry', 'sad'のいずれかを指定可能。 """ #以下は、上記のだるいコマンドをまとめたもの。 open_jtalk = [OPENJTALK_BINPATH + '/open_jtalk'] mech = ['-x', OPENJTALK_DICPATH] htsvoice = ['-m', OPENJTALK_VOICEPATH.format(emotion=emotion)] speed = ['-r', str(speed)] outwav = ['-ow', os.path.join(output_dir, output_file)] cmd = open_jtalk + mech + htsvoice + speed + outwav c = subprocess.Popen(cmd,stdin=subprocess.PIPE) c.stdin.write(text.encode('utf-8')) c.stdin.close() c.wait() if __name__ == '__main__': make_wav("おはようございます。僕は元気です。", speed=0.8, emotion='happy', output_file='test.wav', output_dir=os.getcwd()) これを、実行すれば、.wavファイルが生成される。 このpythonファイルに書かれていることは簡単で、上記のだるいコマンドを、ただまとめただけである。 ついでに、発話のスピードと、感情を変更させることができるようにしてある。 OpenJTalkインストール時に、場所を覚えておけ、といったのは、ここで必要になってくるからである。 pythonを日ごろから扱っている人にとっては、"pythonっぽくない泥臭いコードだなぁ"と思われるかもしれないが、多分これがpython + OpenJTalkの使い方だと思う。(有識者の方、いらっしゃいましたら、コメントください。) 以上で、OpenJTalkのことを完全に理解したと思う。

はじめに NumPyについて、基本的な使い方をまとめました。 動作環境 OS：MacOS 11.2.3(Big Sur) Python:3.9.0 NumPy:1.20.2 基本操作 NumPy配列の生成 以下の例のように、list配列からNumPy配列を生成する方法があります。 sample01 >>> import numpy as np >>> a = np.array([0, 0.5, 1.0, 1.5, 2.0]) #list配列から生成 >>> type(a) numpy.ndarray >>> a = np.array(list('abc')) #文字列のリスト配列から生成 >>> a array(['a', 'b', 'c'], dtype='<U1') 単純な数値のシーケンスを要素として生成するためにarange()関数が用意されています。 numpy.arange(start, end, step) startからend - 1までの数字の順列をstep飛ばしで作成します。startは何も指定しないと0となります。stepは何も指定しないと1です。 sample02 >>> import numpy as np >>> a = np.arange(8) >>> a array([0, 1, 2, 3, 4, 5, 6, 7]) >>> a = np.arange(1, 8) >>> a array([1, 2, 3, 4, 5, 6, 7]) >>> a = np.arange(1, 8, 2) array([1, 3, 5, 7]) >>> a = np.arange(8, dtype=float) #データ型を指定できる。 array([0., 1., 2., 3., 4., 5., 6., 7.]) list配列と同様に、スライス機能もあります。 sample03 >>> import numpy as np >>> a = np.arange(8) array([0, 1, 2, 3, 4, 5, 6, 7]) >>> a[:5] array([0, 1, 2, 3, 4]) >>> a[5:] array([5, 6, 7]) ndarrayクラスの組み込みメソッド 統計学でよく使われる平均や分散、標準偏差など、多数の関数をメソッドとして持っています。いくつか例をあげます。 sample04 >>> import numpy as np >>> a = np.arange(8, dtype=float) >>> a array([0., 1., 2., 3., 4., 5., 6., 7.]) >>> a.sum() #総和 28.0 >>> a.mean() #平均 3.5 >>> a.var() #分散　 5.25 >>> a.std() #標準偏差 2.29128784747792 >>> a.cumsum() #先頭要素からの累積和 array([ 0., 1., 3., 6., 10., 15., 21., 28.]) ベクトル演算 NumPyでは、ベクトル演算が定義されています。これによって複雑な演算処理のコーディングもシンプルに作成できます。 sample05 >>> import numpy as np >>> a = np.arange(8, dtype=float) >>> a array([0., 1., 2., 3., 4., 5., 6., 7.]) # 各要素に定数を加算 >>> a + 3 array([ 3., 4., 5., 6., 7., 8., 9., 10.]) # 定数を各要素に加算 >>> 3 + a array([ 3., 4., 5., 6., 7., 8., 9., 10.]) # 同じオフセットの要素の和 >>> a + a array([ 0., 2., 4., 6., 8., 10., 12., 14.]) # 各要素を定数倍 >>> a * 3 array([ 0., 3., 6., 9., 12., 15., 18., 21.]) # 各要素をベキ乗 >>> a ** 2 array([ 0., 1., 4., 9., 16., 25., 36., 49.]) # 同じオフセットの要素の乗算 >>> a ** a array([ 0., 1., 4., 9., 16., 25., 36., 49.]) NumPyに定義されているスタティックの数理関数 数学で使われる指数関数や平方根、logなどもNumPyパッケージに含まれています。 sample06 >>> import numpy as np >>> a = np.arange(8, dtype=float) >>> a >>> np.exp(a) #要素ごとに指数演算 array([1.00000000e+00, 2.71828183e+00, 7.38905610e+00, 2.00855369e+01, 5.45981500e+01, 1.48413159e+02, 4.03428793e+02, 1.09663316e+03]) >>> np.sqrt(a) #要素ごとに平方根 array([0. , 1. , 1.41421356, 1.73205081, 2. , 2.23606798, 2.44948974, 2.64575131]) >>> np.log(a) #要素ごとにlog関数を適用。ただし、要素が'0'についてはWarningが表示され、計算結果も’-inf'となります array([ -inf, 0. , 0.69314718, 1.09861229, 1.38629436, 1.60943791, 1.79175947, 1.94591015]) Pythonでの数理関数にはmathパッケージもあります。単一の要素を計算する場合は、mathの方が計算速度が早いです。pythonのコマンドインタプリタでは、%timeitを実行する関数の前に記載することで、処理速度を計測できます。 sample07 >>> import numpy as np >>> import math >>> %timeit np.sqrt(4) 1.05 µs ± 13.4 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each) #µs = マイクロ秒 (1秒 = 1,000,000µs) >>> %timeit math.sqrt(4) 129 ns ± 6.05 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each). #ns = ナノ秒　　（1µs = 1,000ns) >>> 上記の例を見ても、mathパッケージに比べるとNumPyの方が処理速度が落ちますので、注意しなければなりません。 等間隔の１次元配列の生成 aからbをx等分した場合の配列を作成することができます。 sample08 >>> import numpy as np >>> a = np.linspace(0,5,12) #0から5までを12等分する >>> a array([0. , 0.45454545, 0.90909091, 1.36363636, 1.81818182, 2.27272727, 2.72727273, 3.18181818, 3.63636364, 4.09090909, 4.54545455, 5. ])

Python Slack SDK を選んだ理由 この記事を書いている2021年4月現在、Slack を操作するために Python を使うなら Python Slack SDK が最善ではないかと思います。また変わるかもしれませんが、わざわざ古いものを使うメリットはありません。 Python Slack SDK はドキュメントが非常に充実しています。 また、https://api.slack.com/ にある Python のサンプルコードも Python Slack SDK を使うものになっています。 やりたいこと 古い Slack の投稿を一括削除したい また、前提条件として次の三つがあります。 対象チャンネルの全ての投稿について既に削除権限がある 対象チャンネルはあらかじめ特定されている 使うのは自分だけ コードを書く前の準備 Slack API のトークンを手に入れる必要があります。 古い投稿を削除するので、どれが古いのかを探す権限が必要です。 また、当然に削除権限も必要です。 そうなると、Web API だけ使えれば十分で、User として動けばいいことになります。 必要なメソッドは、 conversations.history https://api.slack.com/methods/conversations.history と chat.delete https://api.slack.com/methods/chat.delete だけです。 そして、この二つを使うために必要な権限は、channels:history と chat:write の二つです。 User Token Scopes にこの二つの権限を設定し、ワークスペースにインストールすると、User OAuth Token を入手できます。 コードを書く エラー処理やログの記録などは省略しています。必要に応じて追加しましょう。 import os import time from slack_sdk import WebClient # トークンは環境変数に入れておくものとする client = WebClient(token=os.environ.get("SLACK_USER_OAUTH_TOKEN")) # チャンネル ID も実行時に与える方がいいかもしれない channel_id = "C1234567890" # 30日よりも古いものを削除する場合 late_ts = str(int(time.time()) - 60 * 60 * 24 * 30) # 指定したタイムスタンプより古いものだけを取得する result = client.conversations_history(channel=channel_id, latest=late_ts) conversation_history = result['messages'] for x in conversation_history: client.chat_delete(channel=channel_id, ts=x['ts']) print("{} deleted".format(x['ts'])) エラー処理を省いたとはいえ、非常に短いコードで済みます。

概要 今回は「Django」で「HTML」ファイルと「Views.py」ファイルでデータを受け渡しする方法を見ていきます。 これを学ぶことで、自分の思い通りにデータを扱うことができ、アプリ開発の幅を広げることができます。 ここではコードのみ載せているので、解説が必要な方はぜひ下記の記事を参考にしてください。 [Django]DjangoでHTMLとViews.pyでデータを受け渡しする 中級者向けの知識なので、「Djangoまだ初心者だよ…」という方は以下の記事などを参考にして、基礎部分の知識をつけてから読んでください！ [Django初心者]DjangoでTodoリストを作成する 早速見ていきましょう！ こーど ではコードを見ていきましょう。 入力用.html . . . <form action="{% url '~' %}" method="post"> {% csrf_token %} <div class="field"> <label>入力１</label><input type="text" name="nyuryoku1" id="nyuryoku1"> <label>入力２</label><input type="text" name="nyuryoku2" id="nyuryoku2"> <label>入力３</label><input type="text" name="nyuryoku3" id="nyuryoku3"> </div> <input type="submit" value="送信"> </form> . . . ```python: views.py def get_nyuuryoku(request): message = 'データ受け取ったよ！' if request.method == 'POST': nyuryoku1 = request.POST['nyuryoku1'] nyuryoku2 = request.POST['nyuryoku2'] nyuryoku3 = request.POST['nyuryoku3'] params = { "nyuryoku1": nyuryoku1, "nyuryoku2": nyuryoku2, "nyuryoku3": nyuryoku3, "message": message, } return render(request, '~', params) ```html: 表示用.html . . . <div class="field"> <h4>{{ nyuryoku1 }}</h4> <h4>{{ nyuryoku2 }}</h4> <h4>{{ nyuryoku3 }}</h4> <h4>{{ message }}</h4> </div> . . . 終わりに 今回はDjango」で「HTML」ファイルと「Views.py」ファイルでデータを受け渡しする方法を見てきました。 具体的に動きを確認できなかったので、わかりにくかったと思います。 今後実際の動きを確認できるようにアップデートしていこうと思うので、待っていてください！ それでは！

data = get_query_result(293) result = {} for row in data['rows']: add_result_row(result, {'hoge' : row['email']}) add_result_column(result, 'hoge', '', 'string') とてもシンプルに書いた↑ ドキュメントは見つからないので、ここを見るしかなさそう https://github.com/getredash/redash/blob/master/redash/query_runner/python.py

元記事：https://www.octoparse.jp/blog/build-a-web-crawler-with-selenium-and-python/ Webサイトから大量のデータをできるだけ早く取得する必要があるとします。それぞれのWebサイトに手動でアクセスして、コピペでデータを取得することなく、どうやって自動的にデータを取得するのでしょうか？その答えが「Webスクレイピング」です。Webスクレイピングに通じて、この作業を自動化にしてくれます。 今回はPythonでWebサイトからデータをクローニングして、CSVファイルに書き込むというようなWebクローラーを実際に作成してみましょう。 一、必要なPython開発環境を導入 Pythonには、プログラムを組むために便利な標準ライブラリが数多くあります。今回は、以下のライブラリを使用しています。 ・Selenium ーー ブラウザを自動的に操作するライブラリです。主にWebアプリケーションのテストやWebスクレイピングに利用されます。 ・BeautifulSoup ーー HTMLおよびXMLドキュメントを解析するためのPythonパッケージです。 ・csv ーー CSVフォーマットで保存するために使用されます。 したがって、プログラミングを実戦する前に、以下の準備が必要となります。 ・Python 2.xまたはPython 3.xの環境 ・Selenium、BeautifulSoup、csvライブラリの導入 ・Google Chrome ブラウザ それでは、早速始めましょう！ 二、ページ要素（HTMLドキュメント）の検証 Webサイトから、要素（HTMLドキュメント）を検証し、データがどのように構成されているかを分析する必要があります。HTML基礎知識はこちらのリンクで参照ください。今回はテーブルのデータを抽出するように試してみましょう。 Ctrl + Shift + I を押して、テーブルの要素を指定すると、HTMLのソースコードに表示されます。 したがって、テーブルの要素名は「table」と分かっています。 なお、Webクローラーを構築するたびに、HTMLドキュメント内の要素を定位するため、xPathの知識が必要となります。xPathのチュートリアルはこちらのリンクからアクセスできます。 三、コードを書く 1.まず、必要なライブラリをすべてインポートしましょう。 web-scraping.py import csv # csvモジュールをインポート from selenium import webdriver # selenium webdriver をインポート from bs4 import BeautifulSoup # BeautifulSoup をインポート 2.Webdriverを使用する前に、chromedriverへのパスを設定する必要があります。 ※/path/to/chromedriverをWebdriverのパスに変更してください。 web-scraping.py driver = webdriver.Chrome("/path/to/chromedriver") 3.以下のコードを参照してURLを開いてください。 web-scraping.py driver.get("http://test-sites.octoparse.com/?page_id=192") 4.URLを開くためのコードが書けたところで、いよいよWebサイトからデータを抽出します。 先に述べたように、抽出したいデータは要素に入っています。データを持つ要素を探し出し、データを抽出します。以下のコードを参照してください。 web-scraping.py content = driver.page_source BS = BeautifulSoup(content, "html.parser") table = BS.findAll("table", {"class":"wp-block-table is-style-stripes"})[0] # テーブル"wp-block-table is-style-stripes"を指定 rows = table.findAll("tr") # テーブル中<tr>要素の内容を抽出 print(rows) # 抽出したHTMLドキュメントを検証 最後に、web-scraping.pyで保存します。 四、コードを実行してデータを抽出する コードを実行して、必要なHTMLドキュメントを正しく抽出するかどうかを確認します。 五、データを必要なフォーマットで保存 データを抽出した後、抽出したデータをCSV（Comma Separated Value）形式で保存します。そのため、コードに以下の内容を追加します。 web-scraping.py with open("web-scraping.csv", "w", encoding='utf-8') as file: # ファイル名は「web-scraping.csv」を指定する writer = csv.writer(file) for row in rows: csvRow = [] for cell in row.findAll(['td', 'th']): # tdとth要素をループでファイルに書き込む csvRow.append(cell.get_text()) writer.writerow(csvRow) 六、Pythonでスクレイピングしましょう それは最終的なコードです。 追加した後、もう一度コード全体を実行してみてください。 抽出結果は「web-scraping.csv」というファイル名が作成され、このファイルに抽出されたデータが格納されます。 七、Octoparseでスクレイピングする方法 プログラミングが苦手、或いは英語のコードばかりなので苦手意識を持っている方は、スクレイピングツールのOctoparseはおすすめします。 Octoparseは「自動識別」機能があるので、ページのURLを入力するだけで、Webページ上各項目のデータ（テキストとリンクを含む）、「次のページ」ボタン、「もっと見る」ボタン、およびページのスクロールダウンを自動的に検出し、タスク（Webクローラー）を自動的に生成することができます。 早速ですが、Octoparseで自動化の魅力を体験してみましょう。 1.Octoparseを起動して、スクレイピングしたいWebページのURLを入力します。 「抽出開始」 ボタンをクリックして進みます。 2.Octoparseでページが読み込みされたら、自動的にページ上の内容を識別します。 自動識別とは、自動的にページ上の必要なデータを検出して識別するという役立つ機能です。ポイント&クリックをする必要はなく、Octoparseは自動的に処理します。 14.Octoparse_自動識別.png 3.識別が完了すると、データプレビューで識別したデータを表示され、確認してから「ワークフローの生成」ボタンを押します。 4.これで簡単にWebクローラーが作成しました！ 上の「実行」ボタンをクリックして、すぐデータを抽出できます。簡単ではないでしょうか。 八、まとめ Pythonでスクレイピングはそんなに簡単ではないので、学ぶ時間がなく、効率的にスクレイピングがしたい、プログラミングが苦手、或いは英語のコードばかりなので苦手意識を持っている方はスクレイピングツールはおすすめです。 関連記事 PHPで簡単なWebクローラーを作ってみた PythonによるWebスクレイピングを解説 Python vs Octoparse！初心者向きのスクレイピング方法はどっち？ 【初心者向け】ExcelとVBAでWebスクレイピング実戦！

本書は筆者たちが勉強した際のメモを、後に学習する方の一助となるようにまとめたものです。誤りや不足、加筆修正すべきところがありましたらぜひご指摘ください。継続してブラッシュアップしていきます。 © 2021 NPO法人AI開発推進協会 本書は物体検出モデルであるFaseter R-CNNについて説明します。(CNNの基礎を理解している前提で記載しています。まだ理解していない方は別冊のCNNの基礎を先に読んでください） 【参考文献、サイト】 + Faseter R-CNN原論文 + Faster R-CNNにおけるRPNの世界一分かりやすい解説 + 物体検出についての歴史のまとめ（１） + Faster R-CNN: Down the rabbit hole of modern object detection 物体検出の代表的なものにR-CNN、YOLO、SSDなどがあげられます。R-CNNはMicrosoft社によって開発されたアーキテクチャで、R-CNN、Fast R-CNN、Faster R-CNNと発展してきました。説明の際は、R-CNN、Fast R-CNN、Faster R-CNNと順を追って説明していくものが多いですが、R-CNNやFast R-CNNは精度・性能面からもはや使うことはないため、本書ではFaster R-CNNに特化して説明します。過去のモデルの技術を踏襲しているものについては補記して説明します。また物体検出の一般的な説明については、本シリーズ別冊のSSDの説明書を参考にしてください。 １．Faster R-CNN概要 Faster R-CNNは、2015年にMicrosoft社が開発した、Deep LearningによるEnd-to-Endの学習（※１）に初めて成功した物体検出モデルです。（かなりおおまかですが）流れは以下の通りです。 ① 画像上のある矩形の中身が物体なのか背景なのかを学習 ② ①で物体を検出した場合は、具体的に何（ex.犬、猫）が写っているのかを学習 上記①でRPN(Region Proposal Network)と呼ばれるCNN構造を用いている点がFaster R-CNNの特徴です。これにより従来モデルでネックとなっていたSelective Search（※２）処理が不要となり、大幅なスピードアップと精度向上を実現しました。 ※１：End-to-Endの学習とは、入力（この場合、画像）から目的の結果（出力）を直接単一モデルで学習するる学習方法です。Faster R-CNNでは、物体の境界予測と分類を（事前にRPMを学習したうえで）１つのネットワークとして学習を行います。 ※２：Selective Searchとは、候補領域を切り出すためのアルゴリズムの１つです。ピクセルレベルで似た特徴をもつ領域をグルーピングしていくもので、高負荷な処理となります。R-CNNやFast R-CNNは、Selective Searchを使って物体検出を行っていました。 ２．Faster R-CNNの構造 図 1　 Faster R-CNN　全体の構造（原論文から引用加筆） 　 Faster R-CNNは、図1の通り以下の４つの処理から構成されています。 　①　入力画像から特徴マップを出力する処理（学習済みVGG16などを流用） 　②　RPNと呼ばれる物体が写っている場所と、その矩形の形を得る処理 　③　ROIと呼ばれる入力を固定長に変換する処理 　④　何が写っているかを判断する分類処理 (1) 特徴マップを出力 最初のステップでは、分類のタスク用に事前トレーニングされたCNN（本書ではVGG16）を使用し、中間層の出力を特徴マップとして使用します。 図 ２．VGG16のネットワーク概略図（参考サイトから引用） 　 上図２のVGG16で、特徴マップは最後のプーリング層の１つ手前にある14×14×512の層を指します（特徴マップを出力した以降の層、つまり7×7×512以降の層はVGG16の画像分類タスクを行う層のためRPNでは使用しません）。 また、特徴マップのサイズは、VGG16の場合、４回プーリングを実施しているため、元の画像の１/16の縦横をもった512次元のデータになります。例えば、元の画像が300×400×3であれば、特徴マップは18×25×512になります。 図 3　VGG16の特徴マップのサイズ (2) RPNで物体の位置を検出 RPNはそれ自体が「ある画像のどこに物体が写っているのか（物体が写っている場所と、その矩形の形）」を検出する機械学習モデルです。何が写っているかまではRPNでは判断しません。 ⅰ) AnchorとAnchor boxes 特徴マップを生成したら次にAnchorを設定します。Anchorは特徴マップ上の各点になり、それぞれのAnchorに対して９つ（ハイパーパラメータで指定）のAnchor boxesを作っていきます。各Ancho boxesと正解（Ground Truth）を見比べた情報がRPNを学習させる際の教師データとなります。 Anchorは矩形を描く際の中心となります。特徴マップが18×25×512の場合、18×25=450個のAnchorが設定され、下図４のように等間隔に設定されます。 図 ４　Anchor boxesのイメージ 　 Anchorは上図４のように等間隔に設定されるため、「物体」がどこにあってもどれかしらのAnchorがその「物体」の中心となることが可能となります（＝矩形の中心候補）。物体の中心位置が決まるとあとは矩形の縦横が決まれば物体の検出枠が描画できますが、この矩形の縦横を決めるのがAnchor boxです。 各Anchorから「基準の長さ」と「縦横比」をそれぞれ決めることで、複数のAnchor boxesを作り出します。例えば、基準の長さ=(64, 128, 256)、縦横比=(1:1, 1:2, 2:1)とすると、下図5のようなAnchor boxesが生成されます。 このとき１つのAnchorに対して、基準の長さ(3要素) ×　縦横比(3要素)=9個の矩形が生成されます。また、各Anchor boxesの面積は揃える必要があります（ex.基準の長さ=64であれば、縦横比＝1:1のときは64×64=4096なので、縦横比=1:2のときは縦45×横91(=4096)の矩形、縦横比=2:1のときは縦91×横45(=4096)の矩形）。 実際には画像からはみ出したAnchor boxesは無視するので、実際は以下のようになります。 図 ５ Anchor boxesの決定 　 Anchor boxesは最大で18×25×9=4050個作られます。ここで様々な形のbox(矩形)が提案されることにより、正解（ground truth）がどんな形であっても正解と似ているboxの候補を生成することができます。あとは一つ一つのAnchor boxと正解（ground truth）を見比べて、その結果をRPNの出力とします。 ⅱ） RPN層での出力について RPN層では、Anchor boxの中身が背景なのか物体か（下図６のcls layerの2k scores）、物体だった場合は正解（ground truth）とどのくらいズレているのか（下図６のreg layerの4k coordinates）」の２項目を出力します。 図 ６　VGG16の場合のRPN層の出力（原論文から引用加筆） 　 上図６において、「k」はAnchor boxesの数です（今回の例では９）。 ① 「背景が物体か」については、ground truthとAnchor boxesのIoUを計算し、「IoU < 0.3」なら「背景」、「IoU > 0.7」なら「物体」とラベルを付けます。そのため、9(Anchor boxesの数) × 2（ラベル） = 18クラスが作られます。なお、「0.3 ≦ IoU ≦ 0.7」のboxについては、背景とも物体とも言えないので学習には使いません。（IoUの閾値もハイパーパラメータで指定） 　 図 ７　 cls layer出力層のイメージ 　 ② 「ground truthであるAnchor boxのズレ」については、「中心ｘ座標のズレ」、「中心y座標のズレ」、「横の長さのズレ」、「縦の長さのズレ」という４つの指標で評価します。そのため、9(Anchor boxesの数) × 4（ズレ) = 36クラスが作られます。 図 ８　 reg layer出力層のイメージ ３．Faster R-CNNの学習 Faster R-CNNの学習は、RPN部分の学習（図1．Aの部分）とFaster R-CNN全体の学習（図1．Bの部分）を交互に行い、それぞれに精度を高めながら学習を進めます。 (1) RPNの学習 ① 入力画像をVGG16等に通して特徴マップを得る ② ①に対し、さらに３×３の畳み込み（512チャネル）を行い、さらに1×1の畳み込みをかけることで、「背景か物体か」用のアウトプットと、「ground truthとのズレ」用のアウトプットの２つを作る ①を入力、②を出力のモデルとして学習します。「物体か背景か」の２値分類問題と、「ground truthとのズレ」の回帰問題を同時に解いていきます。前者はバイナリクロスエントロピー、後者はＬ１ノルム（絶対値誤差）をベースにした誤差関数を適用して学習を行います。 (2) Faster R-CNN全体の学習 ① 入力画像をVGG16等に通して特徴マップを得る（RPNでの学習と同じもの） ② ROI Pooling（※３）を使用し、特徴マップを固定長のサイズにする ③ ②で得たベクトルを4096の全結合層に2回通し、最後にクラス分類用と、矩形ズレ回帰用の２種類の出力を得る ①を入力、②出力のモデルとして、学習します。なお、RPNを学習させるために用いた「3×3の畳み込み ＋ 1×1の畳み込み」は、RPN特有の構造のため、Faster R-CNNで何かを予測する際には使用しません。また、③の全結合層部分の入出力は常に固定サイズとなるため、Faster R-CNNでは入力画像のサイズに依存しない学習、予測が可能になります。 ※３：ROIプーリングとは、Fast R-CNNで採用された不均一なサイズの入力に対して最大プーリングを実行し「固定サイズの特徴マップ」を得る手法です。 ① 領域候補を同じサイズのセクション（断片）に分割する（このとき、セクションの数は出力の次元と同じ） ② 各セクションで最大値を見つける ③ これらの最大値を出力バッファにコピーする この結果、サイズの異なる長方形のリストから、固定サイズの対応する特徴マップのリストを素早く取得できます（出力サイズの次元は、入力の特徴マップのサイズや領域提案のサイズには依存せず、領域候補を分割するセクション（断片）の数だけによって決定される）。 ROIプーリングのメリットの１つは処理速度です。フレームに複数の物体候補がある場合、それらのすべてに対して同じ入力特徴マップを使用でき、一般的に非常にコストがかかるモデルネットワーク処理の初期段階での畳み込み計算を大幅に削減可能なためです。実際の動作は以下のようなになる。 例）8×8の単一の特徴マップの１つの関心領域（左上、右下の座標）=(0, 3)、(7, 8)に対して、 2×2のRoIプーリングを実行 【参考サイト】Region of interest pooling explained ４．おわりに Faster R-CNNの説明は以上になります。Faster R-CNNはYOROやSSDとならび代表的な物体検出モデルです。　是非マスタしましょう。SSDについては、別冊にて執筆しておりますので、興味のある方はぜひ参照ください。

自分用メモ pip ではpip uninstall tensorflow*のようにwildcardが使えなかったから、ワンライナーで書いた。 pip list | grep tensorflow | awk '{print $1}' | xargs -I{} pip -y uninstall {} tensorflowの部分を任意の文字列に変更してください。

概要 少額の軍資金から自動売買をプログラムして増やしていく企画の論文読み編です． バックテストのための記録コードが完了し，プログラムが回っているので その間に論文を読んでいきます． DEEPLへの丸投げなので変な日本語がありますが，ご容赦ください． 元論文 Bitcoin Price Prediction Based on Deep Learning Methods https://www.researchgate.net/publication/339143532_Bitcoin_Price_Prediction_Based_on_Deep_Learning_Methods Bitcoin Price Prediction Based on Deep Learning Methods 深層学習法に基づくビットコイン価格の予測 学習法に基づくビットコイン価格予測 概要 ビットコインは現在人気のある暗号通貨で, 将来性が期待されています. それは時系列の株式市場のようなもので, 指数化されたデータポイントのシリーズです. 私たちは, さまざまな深層学習ネットワークと, 精度を向上させる方法を調べました. min-max正規化, Adam optimizer, windows min-max正規化などです. 1分ごとのビットコイン価格のデータを収集し, ビットコイン価格を時間単位で反映させるように並べ替え, 合計56,832ポイントを取得しました. 24時間分のデータを入力とし, 次の時間のビットコイン価格を出力しました. 各モデルを比較した結果, メモリがないため, 多層パーセプトロン（MLP）は, 現在のトレンドに基づいて価格を予測するケースには不向きであることがわかりました. Long Short-Term Memory (LSTM)は, 過去の記憶とGated Recurrent Net- work (GRU)がモデルに含まれている場合, 比較的最良の予測を提供します. キーワード 深層学習モデル，多層パーセプトロン，Gated Recurrent Network，Long Short-Term Memory，クロスバリデーション，正規化 1. はじめに ビットコインは, 暗号通貨であり, 電子現金の一形態です. ビットコインは, ピア・ツー・ピアのビットコイン・ネットワーク上で, 仲介者を介さずに, ユーザーからユーザーへと送信できるデジタル通貨です. Bitcoinは, ピア間の取引記録を保持し, すべての記録は暗号化されています. 新しい記録は, 以前のブロックの暗号ハッシュを含んでいます. 各レコードには, タイムスタンプと, 送信者, 受信者, 金額のデータが含まれています. ビットコインはまだ生まれたばかりの技術であるため, ビットコインの将来的な価値を予測することはほとんどできません. GreavesとAuは, 線形回帰, ロジスティック回帰, サポートベクターマシンを使用してビットコインの将来価格を予測しましたが, その精度は低いものでした[1]. Indiraらは, 翌日のビットコイン価格を予測するために, Multi-layer Perceptron based non-linear autoregressive with External Inputs (NARX)モデルを提案しました[2]. また，Jakob Aungiersは，S&P500の株価を予測するために，長期・短期記憶型のディープニューラルネットワークを提案しています[3]. 彼の研究は, 株価に類似したビットコインの予測に光を当てています. Madanらは, ビットコイン予測問題に取り組むために, 一般化線形モデルやランダムフォレストなどのより多くの機械学習アプローチを使用しました[4]. 上記の研究は, 翌日のビットコイン価格を予測することに焦点を当てています. しかし, ビットコインははるかに小さい間隔で頻繁に取引されています. 本研究では, 過去のデータを利用して, 翌日の価格ではなく, 1時間後の価格を予測することを試みています. これは, 実世界でのより良いアプリケーションになる可能性があります. まず，min-max正規化や窓付き正規化[5]などのデータ正規化を行い，窓の初期値と変化率に基づいてデータを正規化しました．MLP（Multiple Layer Perceptron），LSTM（Long-Short-Term-Memory），GRU（Gated recurrent units）の各モデルを，クロスバリデーションを用いてテストデータセットで比較した． 2. データセットの探索 本研究で使用するデータは, Kaggle[6]から収集しています. 2012年1月から2018年7月までのビットコインデータを収集しています. タイムスタンプ, Open, High, Low, Closeでの値, BitcoinとUSDでの取引量, 加重価格, 日付があります. 本研究では, 過去24時間の価格を用いて未来1時間のビットコイン価格を予測することに焦点を当てているため, モデルにはタイムスタンプと加重価格のみを使用しています. 加重価格(price weighted)：加重平均をかけた価格 3. 前処理 図1に示すように, データセットは分単位で, 約3,409,920ポイントを含んでいます. 今回は時間単位で価格を予測したため, 1,409,920/60の56,832点のデータが含まれています. このデータセットは, さらにトレーニングセット, 検証セット, テストセットに分けられます. 図2に示すように, トレーニングデータはデータセット全体の80%を占め, バリデーションとテストはそれぞれ10%です. 時系列データであるため，サンプルはランダムではない．最初の24時間のビットコイン価格を入力として, 次の時間のビットコイン価格を予測した. データ処理とモデル収束の効率を高めるために, 他にもいくつかの前処理方法を導入しています. ミニバッチは, 大きなデータを小さなバッチに分割するために使用され, メモリ効率を向上させます. Minimum-Maximum正規化とWindow-based正規化は, トレーニングデータセット全体を(-1, 1)スケールに設定するために使用されます. ウィンドウベースの正規化は, 株式市場の基準に基づいています. この正規化手法は, 各サイズのウィンドウを取得し, ウィンドウの開始時間からの変化率を反映して各ウィンドウを正規化する[3]. 4. モデル 深層学習ネットワークは, 与えられたデータセットの中からパターンを見つけ出し, それに応じて入力を分類するコンピュータモデリングの一種です. 深層学習ネットワークには, 線形活性化関数を持つMLP（Multiple Layer Percep-tron）や, 次の計算に影響を与える隠れユニットを個別に記録するRNN（Recurrent Neural Network）など, さまざまな構造があります. RNNの例としては, LSTM（Long Short-Term Memory）やGRU（Gated Recur-rent Model）などが挙げられます. MLPは, 予測の基本となる手法です. MLPは, すべての入力を順不同で読み込み, 独立変数と従属変数の関係を判断します. 入力層と出力層の間に隠れ層を追加し, 活性化関数を追加することで, 非線形関係をよりよく表現することができます. RNNは, 以前のモデルの結果と新しい入力データから製品を計算する方法のグループです. 実際, MLPにとっては, 前回の計算結果の「経験」が計算に影響を与えることが良いとされています. この「経験」は, モデルから得られたもので, 秘密にされていますが, 次のモデルに引き継ぐことができます. この秘密の変数は「隠れた状態」と呼ばれ, 現在の計算から未来の計算へと引き継がれます. アルゴリズムとは別に, モデルの出力を独立して決定します. しかし, RNNモデルは, 学習用のデータを入力するために, 時系列のように連続した流れに依存します. 長期的にしか繰り返されないパターンの場合, 前回の繰り返しが次の繰り返しに影響するほどの影響力を持たない可能性があります. また, データを時間順に並べる必要があります. そのため, RNNはMLPと異なり, ランダムなサンプルを与えることができません. 長短記憶は, 遠くの出来事がRNNのネットワークに与える影響が小さくなるという問題を解決します. これは, 記憶する特定のイベントを選択できるスイッチを持っています. また, 長期依存ではないので, 学習量も少なくて済みます. 4つの層で出力を決定し, その結果を持つ隠れた状態を次のサイクルに渡します. "また, 4つの層に加えて「忘却ゲート」が存在し, 経験をカウントすべきでないかどうかを判断します. 4つの層と忘却ゲートには, 短期記憶と長期記憶のどちらに焦点を当てるか, 異なる情報を与えることができます. GRU（Gated Recurrent Model）は, LSTMに比べてシンプルなモデルの一つとされており, 「忘れる」ステップと「入力」ステップを1つにまとめ, 結果として隠れユニットが1つで済むようになっています. この3つの手法のうち, MLPはそのシンプルさと少ない計算量で済むことが評価されています. これらは, 入力として同じ量の情報を持っています. しかし, 隠れた層の数と隠れたユニットは, より多くの魔法の数字です. 隠れ層の数や隠れユニットの数はマジックのようなもので, うまくいく場合もあれば, 逆にうまくいかない場合もあります. RNNは, 隠れユニットを介して前のモデルを計算します. その値は計算に使われますが, 介入する必要はありません. このモデルは, 大量のトレーニングセットを持っていることから, 非常に正確になります. しかし, 長期的なパターンを記憶することができないため, 特に近年のような急激な変化が起こった場合には陳腐化してしまう可能性があります. LSTMは, 以前の状態を「忘れる」かどうかを選択することができます. そのため, 長い時間をかけて繰り返される傾向のあるデータを扱うのに適しています. GRUモデルは, 過去の経験を思い出すかどうかを選択することができますが, より高速に学習することができ, 必要なリソースも少し少なくなります. 本研究では, 6つのモデルを比較しました. モデルの設定を以下の表1に示し, 学習結果については次のパートで説明します. 5. 結果 図3に示すように, MLPとRNNのフレームワークでは, ウィンドウベースの正規化が全データセットベースの正規化よりもはるかに優れているという同様の結論が得られました. 時系列データの特徴から, RNNフレームワークはMLP手法よりも早く収束します. 本研究におけるモデルの性能は, データセットの予測価格と真の価格のRMSE（Root Mean Square Error）によって評価されます. その結果を以下の表に示します. 表2に示すように, 窓法による正規化の方がはるかに良い結果が得られています. テストデータセットの真値に対する予測価格を図4に可視化し, 図5では予測価格をより詳しく見るためにズームインしています. LSTMと窓による正規化の組み合わせが最適であることがわかります. すべてのモデルに対して, 10倍のクロスバリデーションを行いました. 図6に示すように, 学習セットを大きくすると, 誤差が小さくなることがわかります. しかし，時系列データの終盤で揺らぎが大きくなると，再び誤差が少し大きくなる. クロスバリデーションの結果, 表3にまとめたように, 2層のGRUが最も優れており, 2層のLSTMが非常に近い性能を示しています. 6. 結論と考察 クロスバリデーションの結果によると, オリジナルのテストデータセットでは, 2層のLSTMが最も良い性能を示し, 2層のGRUが最も良い性能を示しました. 6つのモデルはどれも近い性能を持っているため, 異なるシナリオでは異なるモデルが好まれる可能性があります. MLPモデルは，RNNモデルよりもわずかに性能が低いものの，より少ない計算能力で済む．本研究では, 時間ベースの予測, 過去24時間のデータの使用, ウィンドウによる正規化, 異なる層数のモデルの比較など, いくつかのユニークな特徴を組み合わせています. 今回の研究をもとに, 今後は, より一般的なビットコイン取引のシナリオに適用できるよう, 一連の推定値を予測する研究を行うことができます. おわりに 時系列解析の知見がほとんど，どんなネットワーク構造になるのかあまり想像できてない． 次回はその辺りを探っていく．

はじめに PRIME STUDY (プライム・スタディ) python3エンジニア認定試験 模擬試験 第一回 問31の「エラーと例外ーユーザー定義例外」 で勉強になったことをメモ。 目次 構文エラーと例外 31問目 コード 解説 参考文献 構文エラーと例外 構文エラー:構文の解釈上のエラーで、エラーは検知されたときに出る^マーク前にあることを示している。 例外:構文エラーではないが、エラーとして検知されたもの。必ずしも致命的ではない。 31問目 コード Python.ipynb class OurException(Exception): #(1)独自の例外 pass def raise_her_exception(a): #(2) print(a, "is a") raise OurException #raise:指定した例外(ここでは(1))を発生させる※1 print("easygoing person.") def func(key: int): #(3)key(int型)を引数とする関数 try: if key == 0: raise_her_exception("Saya") #(2)へ→※1のエラーがでる→exceptへ(例外処理) except OurException as e: print("intelligent") raise Exception #raise:指定した例外(ここではException)を発生させる※2 key = 0 try: func(key) #(3)へ→※2のエラーがでる→exceptへ(例外処理) except Exception as f: print("speedster.") Saya is a intelligent speedster 解説 (1)独自の例外 まず、OurExceptionで独自の例外を書きます (3)key(int型)を引数とする関数 key :intで、型をintに指定(引数が実際にintじゃなくても動く) 処理は以下のようになされる。 1. クラスと関数の定義(1),(2),(3) 2. keyの定義 3. func(key)→(3)へ ここからは上のコードの数を追っていけば流れがわかる。 Key コードの(2)のraise OurExceptionをコメントアウトすると、以下のように出力される Saya is a easygoing person. ここではraiseがコメントアウトされたため、(3)で記載している「※1のエラーがでる→exceptへ(例外処理)」の流れがなかった。 したがってその次のexcept OurException as e:に行かず、 結果的にexcept Exception as f:print("speedster.") も処理されなかったわけである。 参考文献 Pythonチュートリアル>>エラーと例外 Pythonチュートリアル 第3版 PRIME STUDY

はじめに 3次元以上のベクトルに対して直交座標から極座標へ変換するメソッドが見つからなかったので自分で作成しました。せっかくなので公開します。pythonのライブラリで直交座標から極座標へ変換するメソッドが既にある場合には紹介していただけると嬉しいです。 なお、このメソッドを用いて生じた問題については一切責任は負いません。 直交座標と極座標の関係 まずは直交座標と極座標の関係を数式で表現します。多次元の極座標変換を参考にしました。 $n \in \mathbb{N}$として、$(x_1, \cdots, x_n) \in \mathbb{R}^n, (r, \theta_1, \cdots, \theta_{n-1}) \in \mathbb{R}_{\geq 0}\times [-\pi, \pi)^{n-1}$とします。このとき、直交座標と極座標の変換公式は、 \begin{align} x_1 &= r{\rm cos}\theta_1 {\rm cos}\theta_2 \cdots {\rm cos}\theta_{n-2}{\rm cos}\theta_{n-1} \\ x_2 &= r{\rm cos}\theta_1 {\rm cos}\theta_2 \cdots {\rm cos}\theta_{n-2}{\rm sin}\theta_{n-1} \\ x_3 &= r{\rm cos}\theta_1 {\rm cos}\theta_2 \cdots {\rm sin}\theta_{n-2} \\ \vdots \\ x_n &= r{\rm sin}\theta_1 \end{align} です。これを$(r, \theta_1, \cdots, \theta_{n-1})$について解くと、 \begin{align} r &= \sqrt{\sum_{i=1}^n x_i^2} \\ \theta_1 &= {\rm arctan}\left(\frac{x_{n}{\rm sin}\theta_2}{x_{n-1}}\right) \\ \theta_2 &= {\rm arctan}\left(\frac{x_{n-1}{\rm sin}\theta_3}{x_{n-2}}\right) \\ \vdots \\ \theta_{n-1} &= {\rm arctan}\left(\frac{x_2}{x_1} \right) \end{align} となります。この変換公式を使って、直交座標から極座標に変換するメソッドをpythonで作成します。 実装 方針 引数が直交座表現のベクトル、出力が極座標表現のベクトルとなるメソッドを作成します。なお、引数、返り値ともにNumPy.arrayとします。 座標系の変換をするにあたり、以下の場合分けを行います。 ノルムが0のとき 入力されたベクトルの次元数が1以下のとき それ以外 入力されたベクトルの次元数が1のとき 入力されたベクトルの次元数が2のとき それ以外 上記のそれぞれの場合について、返り値を以下のように定めます。 1-1. 文字列「2次元以上のベクトルを入力してください」を表示し、Noneを返す。 1-2. 入力されたベクトルの同じ大きさのゼロベクトルを返す。 2. 文字列「2次元以上のベクトルを入力してください」を表示し、Noneを返す。 3. 極座標ベクトル $\left(\sqrt{x_1^2+x_2^2},\ {\rm arctan}\left(\frac{x_2}{x_1}\right) \right)$ を返す。 4. 極座標ベクトル $\left(\sqrt{x_1^2+x_2^2},\ {\rm arctan}\left(\frac{x_{n}{\rm sin}\theta_2}{x_{n-1}}\right),\cdots, \ {\rm arctan}\left(\frac{x_2}{x_1}\right) \right)$ を返す。 実装例 上記の場合分けに基づいて直交座標から極座標に変換するメソッドは下記です。 get_pc.py #import import numpy as np from math import atan2, sin #直交座標⇒極座標変換メソッド def get_pc(cartesian_coordinate_vec): cartesian_coordinate_vec = np.array(cartesian_coordinate_vec) #NumPy arrayへの変換 dim = len(cartesian_coordinate_vec) #入力ベクトルの次元数を取得 r = np.linalg.norm(cartesian_coordinate_vec) #入力ベクトルのノルムを取得 if r == 0: #上記1に該当 if dim <= 1: #上記1-1に該当 print("2次元以上のベクトルを入力してください") #2次元以上のベクトルしか極座標変換できない else: #上記1-2に該当 return(np.zeros(shape=dim)) #引数と同じ大きさのゼロベクトルを返す elif dim == 1: #上記2に該当 print("2次元以上のベクトルを入力してください") #2次元以上のベクトルしか極座標変換できない elif dim == 2: #上記3に該当 polar_coordinate_vec = np.array([]) #出力ベクトルを入れるベクトルを用意 polar_coordinate_vec = np.append(polar_coordinate_vec, r) #出力ベクトルにノルムを格納 polar_coordinate_vec = np.append(polar_coordinate_vec, atan2(cartesian_coordinate_vec[1], cartesian_coordinate_vec[0])) return(polar_coordinate_vec) else: #上記4に該当 polar_coordinate_vec = np.array([]) #出力ベクトルを入れるベクトルを用意 polar_coordinate_vec = np.append(polar_coordinate_vec, atan2(cartesian_coordinate_vec[1], cartesian_coordinate_vec[0])) for i in range(2, len(cartesian_coordinate_vec)): tmp = atan2(cartesian_coordinate_vec[i]*sin(polar_coordinate_vec[i-2]), cartesian_coordinate_vec[i-1]) polar_coordinate_vec = np.append(polar_coordinate_vec, tmp) polar_coordinate_vec = np.append(polar_coordinate_vec, r) polar_coordinate_vec = polar_coordinate_vec[::-1] return(polar_coordinate_vec) 実行例 実際に上記のメソッドを実行してみました。 test0 = [1, 1] print(get_pc(test0)) このとき、表示結果は(1.414, 0.7853)となり、おおよそ$(\sqrt{2}, \pi/4)$となったので問題ないと判断しました。他の実行例も下表にまとめておきます。 引数 表示 返り値 [] 2次元以上のベクトルを入力してください None [0] 2次元以上のベクトルを入力してください None [1] 2次元以上のベクトルを入力してください None [0, 0] - [0, 0] [1, 1] - [1.4142, 0.7853] [1, 2, 3] - [3.7417, 0.9303, 1.1071] [1, -1, 1, -1] - [ 2.000, 0.5236, -2.5261, -0.7854] まとめ 今回は多次元ベクトルに対して直交座標から極座標に変換するメソッドを作成しました。しかし、自作のメソッドは不安なので、どなたか3次元以上のベクトルに対しても直交座標から極座標に変換するメソッドをご存知の方がいらっしゃたら教えてください。 参考文献 多次元の極座標変換

GitHubのエラー Githubの操作に慣れるためにファイルに追加変更を記入したところで、基本操作である add, commit, pushを行ったところエラーが出てきた。 エラーコード error: failed to push some refs to 'https://github.com/~~~~~' hint: Updates were rejected because the tip of your current branch is behind hint: its remote counterpart. Integrate the remote changes (e.g. hint: 'git pull ...') before pushing again. hint: See the 'Note about fast-forwards' in 'git push --help' for details. 自分の中では特に変な操作をした覚えがないのでなんでエラーが出るかわからず、とりあえず 「error: failed to push some refs to」の部分をググった。 検索先のトップにこちらの記事が出てきた参考にさせてもらったところ、確かに自分もREADMEを操作していたのでエラーの原因はREADMEを操作したことにあると仮定した。 その上でこの記事を参考にして解決してみた。 解決できた流れとしては git branch --contains #今自分がどのブランチにいるのか確認 ↓ master #masterと出たのでとりあえずマスターブランチにいるのかなと考えた git pull origin master #次にこのコードを記入してみたところ From https://github.com/~~~~~~~~~ * branch master -> FETCH_HEAD Merge made by the 'recursive' strategy. README.md | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) create mode 100644 README.md #上記のコードが出てきたのでとりあえずうまくいったと推測 #ここで１回GitHubブラウザを確認したがファイルの更新はできていなかった。 git push -u origin master #最後にpushしてみたところターミナルにpushが成功したと出てきた。 以上の流れで無事にファイルを変更することができたので一安心できた。 感想 自分の中ではエラーになるようなことをしたつもりがなくても、想定外のところでエラーが出てきてしまうこともあるので、小さなことでもエラーが出てしまうんだなと認識できた。

本書は筆者たちが勉強した際のメモを、後に学習する方の一助となるようにまとめたものです。誤りや不足、加筆修正すべきところがありましたらぜひご指摘ください。継続してブラッシュアップしていきます。 © 2021 NPO法人AI開発推進協会 本書は物体検出の代表モデルであるSSDについて説明します。 （CNNの基礎を理解している前提で記載しています。まだ理解していない方は別冊のCNNの基礎を先に読んでください。） 【参考文献、サイト】 ・ SSD原論文 、日本語翻訳サイト ・ ニューラルネットワーク/ ディープラーニング １．物体検出（Object Detection） はじめに、物体検出とは何かについて説明します。物体検出とは、画像内の物体領域の検出と、物体領域のクラス判定（Classification ex.猫、犬）を同時に行うものです。クラス判定をClassfication、物体領領域の検出をLocalizationと呼びます。 図 1　 ClassficationとLocalization 物体検出では主にバウンディングボックス（以降、BBoxと表記）と呼ばれる矩形で物体を検出します。具体的には、以下の５つの変数を予測します。 class_name（クラス名） bounding_box_top_left_x_coordinate（矩形左上のＸ座標） bounding_box_top_left_y_coordinate（矩形左上のY座標） bounding_box_width（矩形の幅） bounding_box_height（矩形の高さ 　※SSDではX座標、Y座標は中心座標を用います ２．物体検出での評価指数 次に、物体検出を評価するための代表的な評価指数について説明します。物体検出では、検出精度の評価はAP/mAPやIoU、処理速度についてはfpsが主に使われます。 AP(Average Precision) ： クラスごとの平均適合率 mAP(mean Average Precison) ： 平均適合率APの全クラス平均 IoU(Intersection over Union) ： BBoxの一致度（オーバーラップ率、Jaccard係数とも呼ばれています） fps(frames per second) ： 1秒間に何フレーム（画像）処理できるか (1) 二値分類の評価指標について 物体検出の評価指標を説明する前に、基本の二値分類の評価指標について説明します。 二値分類では、予測結果と実際の値との間で以下の４つのパターン（①TP（真陽性）、②TN（真陰性）、③FP（偽陽性）、④FN（偽陰性））が存在し、混同行列として表現されます。 ※ここで１文字目は予測が正解(True)か不正解(False)か、２文字目は予測が正(Positive)か予測が負(Negative)か、を表しています。 この混同行列を使って、以下のような指標を算出し、モデルを評価します。 正解率：すべての事象の中で予測が正解していた割合 　　　　　正解率 ＝ (TP + TN) / (TP + FP + FN + TN) 適合率（Precision）：正と予測した中で実際に正であった割合（上記表の１列目） 　　　　　Precision = TP / (TP + FP) 再現率（Recall）：　実際に正であるものの中で、正と予測できた割合（上記表の1行目） 　　　　　Recall = TP / (TP + FN) Ｆ値： 適合率と再現率を調和平均 　　　　　Ｆ値　= 2 × [ (適合率 × 再現率) / ( 適合率 + 再現率)] (2) 物体検出の評価指標について 次に一般的な物体検出の評価指標について説明します。物体検出では、クラスの分類結果に加えて物体領域の検出結果を含めた評価を行う必要があります。 ⅰ) 検出したBBoxの評価 検出したBBoxと正解のBBoxとの一致度をIoU(Intersection over Union)と呼ばれる指標を用いて評価します。IoUは以下のように算出し、重なり度合いを0～1の間の数値で表現します。閾値を設定し、検出が有効か（TP）か、そうでないか（FP）を判断します（閾値以上が有効）。 図 2　 IoU(Intersection over Union) ⅱ) AP・mAPの算出 APは、PR曲線で描かれた下側の面積として算出するのが一般的な方法です。PR曲線は、信頼度が高い順にサンプル数を増やしながら、複数の適合率と再現率を算出し、縦軸を適合率、横軸を再現率としてプロットすることで作成します。 また、mAPはAPの全クラス平均として算出します。 　　　　　出展元）【物体検出】mAP ( mean Average Precision ) の算出方法 ３．物体検出モデルの種類 物体検出モデルは分類モデルとして検出を行うモデルと、回帰モデルとして検出を行うモデルとがあります。分類モデルとしては、R-CNN、SPP-net、FAST R-CNN、Faster R-CNNがありますが、これらの手法は「最初に物体候補が生成」され、次に「これらの候補を分類／回帰に送る」というパイプラインを構築することで分類問題として扱っています（詳細はそれぞれの論文や関連記事を参照ください）。　 一方回帰問題として扱うモデルはいくつかありますが、最も代表的なものに「YOLO」と今回説明する「SSD」があります。 以降、SSD（Single Shot multibox Detector）について説明をします。 ４．SSDについて SSDは、様々な階層の出力層からマルチスケールな検出枠を出力できるよう設計されています。　　 モデルアーキテクチャは以下の通りです。 図 4　SSDのアーキテクチャ SSD の基本的なアーキテクチャは、上図のように、フィードフォワード型（順方向）の多層 CNN をベースに構成されています。 ネットワークの最初の部分のレイヤーは、画像分類に使用される標準的なアーキテクチャ（上図では VGG-16）に基づいて構成されており、これをベースネットワークといいます。 このベースネットワークで、特徴量を検出します（特徴量はより広い領域から抽出）。 その後のレイヤーは、マルチスケール特徴マップによる多様な物体検出のための補助的な構造です。 [補足] 畳み込み層の出力「Detections 8732 per Class」について 「8732」はデフォルトモデルでのバウンディングボックスの数を示しています。上図の通り、300×300の元画像を38×38、19×19、10×10、5×5、3×3、1×1に均等に分け、それぞれ4個または6個のデフォルトボックスを生成します。具体的には以下の計算式で算出できます（’_’はその層で生成するバウンディングボックスの数）。1枚の画像から8,732のバウンディングボックスでクラス判定を行うことになります。 （38×38×4）+（19×19×6）+（10×10×6）+（5×5×6）+（3×3×4）+（1×1×4）＝ 8,732 ５．SSDの物体検出の概要について 図 5　特徴マップとデフォルトボックス ①それぞれの層のデフォルトボックス（４～6個）ごとに、「Offset情報（自身が物体からどのくらい離れていてどのくらい大きさが異なるのか）、および、クラス情報」を予測します。 ②①でクラス毎にスコアの高いデフォルトボックスを選択します。 ③②で選択したデフォルトボックスを予測値をもとにバウンディングボックスを確定します。 ④クラス毎に重なり率（IoU）を算出し、IoUが高い場合はスコアが低いデフォルトボックスを除します。 top-k filtering：クラス所属の確信度が上位 k 個のもののみを抽出します。 non-maximum suppression アルゴリズム：推論されたデータに対し、バウンディングボックスの重複防止のために non-maximum suppression アルゴリズムを適用します。 図 6　 non-maximum suppression アルゴリズム 以降、特徴について詳細に説明します。 (A)基本コンセプト 「多層CNNでは、conv層やpooling層で、特徴マップがダウンサンプリングされて、後段に行くほど、特徴マップのグリッドサイズが小さくなります。これにより、各々の層の特徴マップで、色々なサイズの物体を検出出来る情報が含まれていることを意味しています。 従って、SSDモデルの後段の特徴マップ（特徴レイヤー）の各グリッドでは大きな物体の情報を、前段では、小さな物体の情報を取得することが出来ます。そして、各グリッドにおける特徴量を使用して、バウンディングボックス（BBox）のアスペクト比、所属クラス、座標のオフセットを学習させる。」というのが、基本的なコンセプトです。 (B)検出のためのマルチスケール特徴マップ 畳み込み特徴レイヤーを、（途中で打ち切られている）ベースネットワークの最後尾に追加しています。これらのレイヤーは、特徴マップのサイズを小さくさせ（上図５）、マルチスケールでの検出の予想を可能にします。（下図７） 図 7　マルチスケール特徴マップ 　 検出を予想するための畳み込みモデルは、各特徴レイヤーにおいて異なっています。（モデル図の青線元のレイヤー） （※ YOLO では、対照的に１つのスケールの特徴マップを扱っています） (C)検出のための、畳み込み予想器 ベースネットワークの後尾に追加された各特徴レイヤーは、畳み込みフィルタ（フィルタ行列）の集合を使用して、固定の検出予想の集合を生成可能です。 p 個のチャンネルをもつサイズ m×n の特徴レイヤーに対しての、潜在的なパラメータ予想のための基本要素は、3×3×p の小さなカーネル（カーネル行列）であり、このカーネルは、予想カテゴリ（所属クラス）の 検出スコア（物体の中心座標が含まれているなら１となるようなスコア）、または、デフォルトボックスの座標に関してのオフセット値を生成（算出）します。 また、バウンディングボックスの座標オフセットの出力値は、各特徴マップの位置に対するデフォルトボックスの位置に対して測定されます。 SSD モデルは、（前述のように）いくつかの特徴レイヤーを、ベースネットワークの最後に追加しますが、これらの特徴レイヤーは、異なるスケールとアスペクト比のデフォルトボックスに対するオフセット値と、それに付随する確信度を予測します。 (D)デフォルトボックスとアスペクト比（デフォルトボックスと回帰によるオフセット予想） ネットワークのトップで、複数の特徴マップに関して、各特徴マップのセル（グリッド）とデフォルトボックスの集合を関連付けています。 その対応するセル（グリッド）に対して、各デフォルトボックスの位置が固定されるように、デフォルトボックスは、畳み込みのやり方で、特徴マップを隙間なく敷き詰め [tile] ています。（下図８参照） 図 8　デフォルトボックスの敷き詰め 　 各特徴マップのセル（グリッド）において、検出したい物体が、各BBox内に存在することを示す指標である検出スコア（物体の中心座標が含まれているなら１となる）と、セル中のデフォルトボックスのオフセット値を予想します。 このオフセット値は、具体的には、ある座標位置で k 個のボックスそれぞれに対し、c 個のクラスの検出スコア（所属クラスの確率）と、元のデフォルトボックスに対するオフセット値 (x, y, w, h) の４つ（x, yの中心座標２つ、幅w、高さh）を計算します。 その結果、特徴マップに適用されるフィルタは、合計 (c + 4)×k 個となり、m × n の特徴マップに対して、(c + 4)×k×mn 個の出力が生成されることになります。 幾つかの特徴マップで、異なるデフォルトボックスの形状を利用することにより、出力されるボックス（出力特徴マップ）の形状を、効率よく離散化 [discretize] することが出来ます。 (E)デフォルトボックスとアスペクト比の選択（デフォルトボックスと回帰によるオフセット予想）の詳細 異なるスケールの物体検出を取り扱うために、元画像を異なるサイズで処理し、それらの結果を結合する手法を提案している研究も存在します。 しかしながら、このSSDのように、単一のネットワークの中の幾つかの異なるレイヤーの特徴マップを利用することで、全ての物体のスケールについて、同じパラメータを共有しながら、同様の効果を得ることが出来ます。（更に、こちらの手法では同一の単一のネットワークの使用するため、処理が軽くなるメリットがあります。） また、先述のように特定の特徴マップが、特定のスケールの物体に対応するように学習させるために、デフォルトボックスを“敷き詰めて”設計されています。デフォルトボックスを中心とする提案領域は、物体のスケール値、中心座標、高さ、幅が合っていないことがあるため、スケール値、幅、高さ、中心座標に回帰する畳み込み層を追加しています（BBoxの形状回帰）。具体的には、今、m 個の特徴マップを予想に使用するケースにおいて、各特徴マップ k についてのデフォルトボックスのスケール は、以下のようにして計算されます。 この式より、最下位のレイヤーのスケール 0.2 と最上位のレイヤーのスケールは 0.9 となり、中間レイヤーのスケールは、上式に従って規則的な間隔で設定されます。 また、デフォルトボックスのアスペクト比に関しては、ar = {1, 2, 3, 1/2, 1/3 } の異なるアスペクト比を設定します。 このスケールとアスペクト比により、各デフォルトボックスに対して、幅 Wk(a) = sk√arと、高さ hk(a) = sk / √ar が設定されます。 (F)特徴マップとデフォルトボックスについて SSD が訓練中に必要とするのは、入力画像（図５の(a)）と、各物体それぞれの正解ボックス（デフォルトボックスの内、各物体が収まるボックス。上図の赤枠と青枠）のみです。 各層での畳み込み処理のやり方 [in a convolutional fashion] において、いくつかの特徴マップでの各位置（中心座標）において、異なるアスペクト比デフォルトボックスの少数の集合（上記例では４個）を、異なるスケールの特徴マップ内（例えば、上記の (b) の 8×8 の特徴マップ内、(c) の４×４の特徴マップ内）で評価します。 そして、これらのデフォルトボックスそれぞれにおいて、形状のオフセットloc :(cx, cy, w, h) と、全ての物体カテゴリー (c1, c2, ..., cp) に関する確信度 confを予測します。 ※w : デフォルトボックスの幅、h : デフォルトの高さ、cx,cy : デフォルトボックスの中心座標 訓練時には、最初にこれらのデフォルトボックスと正解 [ground truth] ボックスのマッチ度を図ります。上図の例では、２つのデフォルトボックスの内、１つ目はネコ、２つ目はイヌとマッチさせていますが、この組み合わせは正 [positive] として扱われ、残りは負 [negative] として扱われます。 モデルの誤差（損失関数）は、 位置特定誤差 [localization loss] （例えば、Smooth L1）と、確信度誤差 [confidence loss] （例えば、softmax）との間の重み付き和 [weighted sum] です。 ６．訓練（学習）アルゴリズム (a)マッチング戦略 訓練では、どのデフォルトボックスが正解ボックスとなるのか決定する必要があり、その結果を元にネットワークを学習させます。 各正解ボックスは、座標位置、アスペクト比、スケール値が異なる幾つかのデフォルトボックスから選択しますが、これらデフォルトボックスに対して、jaccard overlap の最良値（最もエリアが正解と重複している）で、各正解ボックスのマッチ度（エリアの重複度）を算出することになります。 この際、ベストマッチした（最もエリアが重複している）デフォルトボックスだけでなく、jaccard overlap が 0.5 の値よりも大きいデフォルトボックスを正解ボックスと判定させ、学習させます。 これにより、正解ボックスに複数に重なり合っているデフォルトボックスについて、高いスコア予想が可能になります。 (b)「損失関数」 SSD の損失関数は、位置特定誤差（loc）と確信度誤差（conf）の重み付き和であり、（SSD の学習は、複数の物体カテゴリーを扱うことを考慮して行われるため２つの線形和をとる。）以下の式で与えられます。 ・N：正解ボックスとマッチ（一致）したデフォルトボックスの数 　　　ただし、N = 0 の場合（一致するデフォルト）は、loss値を0とする ・l：予想あれたボックス ・g：正解ボックス 　 ここで、上式の位置特定誤差 Lloc は、予想されたボックス（l）と正解ボックス（g）の間の Smooth L1 誤差（関数）であり、以下の式で与えられます。 ※Smooth L1 誤差（関数） 又、確信度誤差 Lconf は、所属クラスのカテゴリ（c）に対する softmax cross entropy 誤差（関数）であり、以下の式で与えられます。 (c)ハードネガティブマイニング [hard negative mining] マッチング工程の後、有効なデフォルトボックスの数が多い場合多くのデフォルトボックスは、負 [negative] に判定され、正と負の訓練データの比率が不均衡になってしまいます。 （＝典型的な画像の多くの面積は、背景によって占められますが、これを有効なデフォルトボックスとして採用すると、背景しか出力しないネットワークでも、ある程度 loss 値を下げることが出来てしまい、結果として検出能力の低いモデルになってしまいます。） この問題に対する対策として、負に判定される全訓練データを使用する代わりに、これらの（訓練データとしての） デフォルトボックスに対しての誤差関数が高い順（降順）にソートし、負と正の比率が、最大でも３：１になるように、誤差関数の値が上位のもののみを選択します。 これにより、より速くモデルが最適化され、又、安定した学習に繋がります。 (d)データ拡張 [data augumentation] モデルを様々な物体の大きさと形状に対して、よりロバスト（堅牢）にするために、各訓練画像は、以下のオプションをランダムに選択し、（画像中の領域の）サンプリングを行います。 元の入力画像全体を使用する。 物体画像との最小の jaccard overlap が、0.1 , 0.3 , 0.5 , 0.7 , 0.9 となるように、画像中の領域（サンプルパッチ）をサンプリングする。 画像中の領域（サンプルパッチ）をランダムにサンプリングする。 この各画像中の領域（サンプルパッチ）のサイズは、元の画像サイズの 0.1倍 ~ 1.0倍で、アスペクト比は 1/2 ~ 1.0 倍です。 但し、サンプルパッチの中に正解ボックスの中心座標が存在する場合は、正解ボックスの重複部分はそのままにします（サイズやアスペクト比を元の画像から変更しない）。 ７．おわりに 以上が物体検出・SSDの基礎になります。今後、分類モデルとして物体検出を行うFaster R-CNNや、最近話題になっているTransformerを使った物体検出モデルであるDETR（DEtection TRnsformer)について執筆予定です。