目次 ・藤崎モデルのパラメータ逆推定をやってみる ・パラメータ逆推定 ・実装 　- ①前処理 　- ②アクセント成分の近似 　- ③フレーズ成分の近似 　- ④AbS最適化 　- ⑤結果 　- ⑥事後処理も ・おまけ 藤崎モデルのパラメータ逆推定をやってみる \begin{align} \log F_0[t] &= F_b + F_p[t] + F_a[t] \\ &= F_b + \sum_i Ap_iG_p(t - T_{0i}) + \sum_j Aa_j\{G_a(t - T_{1j}) - G_a(t - T_{2j})\} \\ G_p(t) &= \alpha^2 t \exp(-\alpha t) \hspace{10pt} (t \geqq 0,\hspace{5pt} else\ 0) \\ G_a(t) &= \rm{min}[1 - (1 + \beta t)\exp(-\beta t),\ \gamma] \hspace{10pt} (t \geqq 0,\hspace{5pt} else\ 0) \end{align} のように、基底成分$F_b$、フレーズ成分$F_p$、アクセント成分$F_a$の3つの和でモデル化可能。これが藤崎モデル（超ザックリ）。 フレーズ成分 こいつは発声してから緩やかに減衰していく成分。 $G_p(t)$はインパルス応答で、$T_{0i}$は$i$番目のフレーズ成分の生起時刻、$A_{pi}$は$i$番目のフレーズ成分の振幅を表す。 $\alpha$はフレーズ成分の固有各周波数で、話者や言語の違いによる差が小さいとされており$\alpha = 3\rm{rad/s}$で定義される。 アクセント成分 こいつは急峻な増減成分を担う。。 $G_a(t)$はステップ応答で、$T_{1j}$と$T_{2j}$、$A_{aj}$はそれぞれ$j$番目のアクセント成分の生起時刻、終わり時刻、振幅に対応する。 $\beta$はアクセント成分の固有各周波数でこちらも話者・言語によらず$\beta = 20\rm{rad/s}$に。 $\gamma$は同様に$\gamma = 0.9$で。 パラメータ逆推定 藤崎モデルがパラメータから音声の基本周波数を生成する順問題なので、音声の基本周波数からパラメータを逆推定する逆問題を解く。 これができると・・・ - 基本周波数を少数のパラメータで表現できる - 基本周波数の加工・操作が容易になる というメリットがある。 ただ、当然だけど逆問題は簡単じゃない。 しかも詳しく提案法が書かれてる論文を最近はあまり見かけないな〜という感じ（たぶん）。 今回は20年くらい前の論文だけど、この辺を再現できればな〜っていう感じ。 （A METHOD FOR AUTOMATIC EXTRACTION OF MODEL PARAMETERS FROM FUNDAMENTAL FREQUENCY CONTOURS OF SPEECH） 大まかな流れ ※発話区間毎に ① 前処理：抽出したF0のグロスエラーを直し、マイクロプロソディを除去し、無声区間を補間して、平滑化する。 ② アクセント成分の近似 ③ フレーズ成分の近似 ④ Analysis-by-Synthesisで各パラメータの最適化 っていう感じ。 ①前処理 世間の提案法とは違って - F0抽出は、WORLDのHARVESTから抽出後、STONE MASKを適用して得たものとする - マイクロプロソディはローパスフィルタによって除去する ことにしてる。 WORLD信者なのでF0はWORLDで抽出したいし、そうすると抽出過程の情報を用いてマイクロプロソディとかを除去するという論文の手法を使えないからっていうだけの理由。 ちなみにWORLD信者なので（２回目）、抽出したF0は高精度だと信用しきってるため、「グロスエラー？何それ？」の如く 前処理から省いてまーす。 前処理[1] マイクロプロソディの除去 マイクロプロソディは急峻なF0変化のことで、F0軌跡の高周波成分に該当する。 なのでカットオフ周波数10Hzのローパスフィルタで除去する。 （「マイクロプロソディ ローパスフィルタ」でググればいくつか論文出てきます。はい。） def remove_mp(f0, period, n_taps=51, cutoff=10): padd1, padd2 = np.zeros(n_taps) + f0[0], np.zeros(n_taps) + f0[-1] f0_ = np.hstack((np.hstack((padd1, f0)), padd2)) fil = sl.firwin(numtaps=n_taps, fs=int(1./period), cutoff=cutoff, pass_zero=True) f0_ = sl.lfilter(fil, 1, f0_) f0_ = f0_[int(n_taps/2)+n_taps:len(f0)+n_taps+int(n_taps/2)] return f0_ periodはF0のサンプリング周期（WORLDのデフォルトは0.005）。 発散しないまでもフィルタに通すと両端が大振りになるので 最初にf0の両端を伸ばしておく。どれくらい伸ばすかはご自由に（？） また、FIRフィルタはタップ数の半分だけ遅延するので その遅延分を取り戻してあげる（最後の行）。 FIRフィルタじゃなくても自分でFFTしても良いかと。 FFTで振幅スペクトルと位相を求めて、10Hz以上の振幅を1e-20くらいにして iFFTで戻すだけ。 タップ数やパディング数は実際にやってみて良さそうな数字にしてる。 「ローパスフィルタでマイクロプロソディを除去しましょう」っていう話は細かいフィルタ設計までは説明見当たらないので 自分でいじってみるしかなさそう。。。 前処理[2] 夢精区間の補間 F0は無声区間は0で不連続なのでそこの補間が必要だが、HARVESTはそもそも無声区間でも補間してあるのでここはスルー（ハッピー！！！） 前処理[3] 平滑化 分析区間毎に3次関数で平滑化する。このとき分析区間の繋ぎ目でも連続かつ微分可能なように平滑化する。 ここは初学習者でも論文の数式を見れば実装できるかと。英語も優しいし。 ①前処理の結果 ②アクセント成分の近似 先ほど、F0は区間毎に3次関数で近似した。区間の境目でも微分可能にしてあるので、この平滑化F0の"１次導関数の"極大値と極小値を求め（つまり変曲点ですね）、それぞれからアクセント指令を求める。 このとき、極大値が0.6未満、極小値の絶対値が0.4未満の極値は無視する。 参考文献ではそう書かれてないけど載せてる図は明らか閾値処理してるし和文の方にそう書いてあるので 1次導関数の極値を求める def inflection_point(coef, ix1, ix2, timeaxis): dy = coef[:,1] + 2*coef[:,2] * timeaxis + 3*coef[:,3] * (timeaxis**2) dy_ = dy[ix1:ix2+1] # inflection point ix_ = np.hstack((maxid, minid)) ix = np.sort(ix_) label = np.ones(len(ix)) tmp = np.asarray([np.where(ix == i)[0] for i in minid]).reshape(len(minid)) label[tmp] = -1 candidate = ix+ix1 return dy, candidate, label.astype(int) ・coefは3次近似したときの係数で、T-by-4の次元のarray配列。 ・ix1, ix2は発話区間のオンセットとオフセットインデックス ・timeaxisは時間軸の配列で、WORLDで得られるやつ アクセント指令の時刻を求める 極値を順に見ていって、極大値をとる時刻$t_{1j}$から$1/\beta$だけ遡った時刻をオンセット($T_{1j}$と呼ぶ)、極小値をとる時刻$t_{2j}$から$1/\beta$だけ遡った時刻をオフセット($T_{2j}$と呼ぶ)として採用していく。 ただし極大値（または極小値）が連続する場合は最も大きい極大値（または最も小さい極小値）を採用する。 ほいでもし最初が極小値から始まる場合、その極小値をとる時刻から平均1モーラ長遡った時刻を仮想的にオンセットとする。 和文の論文から引用した。 ちなみに今回、最後が極大値なら発話区間のオフセットを仮想的にアクセント指令のオフセットとした。 アクセント指令の振幅を求める まず結論 $j$番目のアクセント指令の振幅$A_{aj}$は \begin{align} A_{aj} &= 0.5(A_{a1j} + A_{a2j}) \\ A_{a1j} &= \frac{\rm{e}}{\beta}G_0(t_{1j}) \\ A_{a2j} &= \frac{\rm{e}}{\beta}\big|G_0(t_{2j})\big| \\ G_0(t) &= a_1 + 2a_2 t + 3a_3 t^2 \end{align} ただし$G_0(t)$の $a_k$ は①前処理で対数F0を3次関数近似した際の近似係数で、 $t_{1j} = T_{1j} + 1/\beta$で、$t_{2j} = T_{2j} + 1/\beta$ 解説 英文ではきちんと記述されてないので和文から数式を持ってくる。 $t_{1j} = T_{1j} + \frac{1}{\beta}$ は対数F0の変曲点を取る時刻であり、1次導関数における極大値を取る時刻になっている。 \begin{align} \frac{\partial}{\partial t} \log F_0[t] &\simeq \frac{\partial}{\partial t} A_{a1j}\big(G_a(t - T_{1j}) - G_a(t - T_{2j})\big) \\ &= X[t] \\ X[t_{1j}] &= \frac{\beta}{\rm{e}}A_{a1j} \end{align} 対数F0の$t = t_{1j}$における1次微分係数がこれ。 なぜかフレーズ成分の項を消し飛ばしてるけど。 一方で、①前処理で対数F0は3次関数近似してるから、その近似係数を用いて \begin{align} \Big(\frac{\partial}{\partial t} \log F_0[t]\Big)[t_{1j}] &= a_1 + 2a_2 t_{1j} + 3a_3 t_{1j}^2 \end{align} と表せるから \begin{align} \frac{\beta}{\rm{e}}A_{a1j} &= a_1 + 2a_2 t_{1j} + 3a_3 t_{1j}^2 \\ &= G_0(t_{1j}) \\ A_{a1j} &= \frac{\rm{e}}{\beta}G_0(t_{1j}) \end{align} となる。 同様に$A_{a2j} = \frac{\rm{e}}{\beta}\big|G_0(t_{2j})\big|$ 負にならないよう絶対値を取る 以上から、$j$番目のアクセント指令の振幅$A_{aj}$は $A_{aj} = 0.5(A_{a1j} + A_{a2j})$ とする。 ③フレーズ成分の近似 アクセント成分を減算した残差成分を積分して、ある値を超え始めた時点にフレーズ指令を立てる 発話区間最初のフレーズ指令 さっき求めたアクセント指令パラメータで合成したアクセント成分$y_a(t)$を$\log F0$から引く（→残差成分$R_0(t)$をゲット)。 R_0(t) = \log {\rm F_0}(t) - y_a(t) 次に、区間 $[T_S, T_{21}]$ において$R_0(t)$の最大値$R_{0max}$を取る位置$T_{R_{0max}}$を求めてる！→第一フレーズ指令の位置$T_{01}$をゲット！ T_{01} = T_{R_{0max}} - \frac{1}{\alpha} 残差$R_0$は$y_p$に依存するんだから$y_p$が最大なら$R_0$も最大でしょ！ってことで・・・ $R_0[t]$の最大値は$R_0[t_{0max}]$だから \begin{align} R_0[t] &= y_p[t] + y_b \\ R_0[T_{R_{0max}}] &= y_p[T_{R_{0max}}] + y_b \\ R_0[T_{R_{0max}}] - y_b & = y_p[T_{R_{0max}}] \end{align} なので、$y_p[T_{0max}]$が区間$[T_S, T_{21}]$での$y_p$の最大値になるはず。 だから数式をごにょごにょして、第一フレーズ指令の振幅$A_{p1}$は A_{p1} = \frac{\rm{e}}{\alpha} \big(R_{0max} - y_b\big) になる。 以降のフレーズ指令 今、$i$個のフレーズ指令が見付かったとする。 現状求めた全てのフレーズ指令を用いて$y_p(t)$を求め、 R_i = \log {\rm F_0} - ya - yp と、残差を求める。 $R_i$が大きかったらまだまだフレーズ指令が存在するということ。 なので、$R$を積分していって$\theta_1$を超えたら$\theta_2 ( < \theta_1)$ を超えた時点に指令を立てる( $T_{0i+1}$の誕生)。 もし$T_{0i+1}$がアクセント成分内だったら当該アクセント指令のオンセットと、一個前のアクセント指令のオフセットの間に立て直す。 本当はもう少し細かい処理をするけど面倒で...ごめんなさい 上図例のように、生起位置候補$T_{0i+1}$が $j+1$番目のアクセント指令内にあった場合（青矢印）、 ``あらゆる指令は同居しない''というルールを破るので、移動させる（赤矢印）。 $T_{0i+1}$は先ほど同様、候補時刻から$1/\alpha$遡った時刻だと仮定して、 $R_{imax}$は$y_p(T_{0i+1}+1/\alpha)$の最大値と一致すると仮定すれば またまた数式をごにょごにょして A_{pi+1} = \frac{\rm{e}}{\alpha} \big(R_{imax} - y_b\big) はい。「以降のフレーズ指令」を繰り返す。 和文では$i+1$番目のフレーズ指令を求めたら、その都度最適化して微調整してるけど面倒なので今回はしてません 積分の閾値は$\theta_1 = 7.0, \theta_2 = 3.0$にしてます。 　閾値は記載ないのでやってみて良さげなもので。 ④AbS最適化 求めた全ての指令を時系列順に１個ずつ最適化する。Analysis-by-Synthesisで、F0を再合成し、3次関数近似したF0との最少二乗誤差を最小化させる。 ただし、$k$番目のある指令$q_k$の最適化は q_k^{\rm{new}} = \rm{argmin}_{q_k} \Big(\log {\rm F_0}(t) - \big\{y_a(t) + y_p(t) + y_b\big\}\Big)^2 \hspace{10pt} (t = [t_1, t_2]) $t_1$は指令の生起位置、$t_2$は指令の影響が無視できる時刻として、この区間内の誤差が小さくなればOK。 t2を探す：アクセント指令の場合 $t$が発散すると$y_a$は0に収束するので、"$y_a$の影響が無視できる時刻"を考えるのが面倒なので、 次の指令のオンセットの手前にしてみた。とりあえず。はい。 もし発話内最後の指令なら次の指令は存在しないので発話自体のオフセットにしておく。 あと！もう一つ！ アクセント指令を更新した時に隣の指令に侵食すると良くないので、 どんなに更新しても現段階の両隣の指令に侵食しないように！ →「もし更新して隣の指令と重なるor跨ぐなら、更新しない」 t2を探す：フレーズ指令の場合 本当はこれもちゃんと考えないといけないけどさっきと同じ理由でサボった。 さっきと同じで。はい。すみません。 最適化する 最適化は「山登り法」で行った。 勾配法でも良かったけど指令の時刻を１インデックス分増やす・減らす・留めるの３パターンしか無いんだから計算時間もかからないし山登り法の法が確実じゃないですか。 ただ、振幅については実数値だから大人しく勾配法を使った。 エポック数は20でいい感じになったよ〜 ⑤結果 推定結果。アクセント指令が「ほんまに？」っていう感じもするけど 音源的にはフレーズ成分は２発しかなさそうな感じなので・・・まぁ・・・ ちなみに最適化せず、初期推定のままだとこんな感じ。 ⑥事後処理も オリジナルのlogF0と3次近似したlogF0の残差を求めておいて、 再構築したlogF0にその残差を加えてあげましょう。 zansa = numpy.log(f0_orijinal) - f0_approximated f0_recon = f0_recon + zansa f0_recon_ = numpy.exp(f0_recon) 本当はこのままやるとlog(0)の警告、log(0)が生んだ-infの演算による警告が出るので ちゃんと0の取り扱いには配慮しましょう。 オリジナルのF0（黒）と再合成したF0（赤）が概形一致してる。 各発話区間の頭とケツがピッタリ合うのは おまけに書きまーす おまけ WORLDの基本周波数推定アルゴリズムを使ってF0を抽出してます。 アルゴリズムは - DIO - HARVEST の２つあります。 HARVESTはDIOと違い連続値でF0を得られますが、 その推定されたF0は、候補生達の中から信頼度が高いためF0として選ばれてます。 HARVESTは頭とケツのF0よりちょい前後の候補生達も、「滑らかに接続できそうやん！」ってなれば もれなく敗者復活でF0として選択するアルゴリズムになってます（多分そんな感じだったかと）。 今回、3次近似するときにそれが邪魔な場合があるので F0の有声区間についてのみDIOのものを使いました。 なので今回の藤崎モデルパラメータ推定はF0よりも少し狭い範囲での物語になります。 再構築するときは残差を用いて微細変動を復元してるため、 その都合で推定時に無視した区間のF0が元のF0と完全一致します。 これが良いのか悪いのかは知りませんが、 その端っこのF0の挙動によって3次近似の概形が変わるのと、 そうすると先々のパラメータ推定にも影響してくるわけで 僕のチープな肌感だとHARVESTのF0延長する仕様は邪魔だと思ったのでこうしました。 おしまい 間違いや改善点などご指摘いただけると幸いです。 個人的に、DeepLearningよりこういうアルゴリズムというか処理の方が好きなので もっともっと勉強していきたい所存です。

1 < x < 10 などの記述ができるんだよ，という記事は多いけど，別に2つの比較演算を単純に記述できるだけじゃないんだね。 1 < 2 <= 2 < 3 == 3 > 2 >= 1 == 1 < 3 != 5 の結果はどうなる？ 順番にやっていけばわかるけど，True だ。

PyTorchにはPyTorch Image Modelsなど学習済モデルがたくさん公開されていて、これを転移学習に使うことも多いです。その際、学習済モデルを少し改造して試したい場合、どうすればいいのか。直接編集するわけではありませんが、同等の効果がある方法をご紹介します。 要点 1. 挿入するにはtorch.nn.Sequentialで置き換える 2. 削除するにはtorch.nn.Indentityで置き換える まずはtorch.hubからお馴染みのresnet18を取得して、その構造を表示してみましょう。 import torch model = torch.hub.load('rwightman/pytorch-image-models:master', 'resnet18', pretrained=True) print(model) Using cache found in /root/.cache/torch/hub/rwightman_pytorch-image-models_master ResNet( (conv1): Conv2d(3, 64, kernel_size=(7, 7), stride=(2, 2), padding=(3, 3), bias=False) (bn1): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act1): ReLU(inplace=True) (maxpool): MaxPool2d(kernel_size=3, stride=2, padding=1, dilation=1, ceil_mode=False) (layer1): Sequential( (0): BasicBlock( (conv1): Conv2d(64, 64, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn1): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act1): ReLU(inplace=True) (conv2): Conv2d(64, 64, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn2): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act2): ReLU(inplace=True) ) (1): BasicBlock( (conv1): Conv2d(64, 64, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn1): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act1): ReLU(inplace=True) (conv2): Conv2d(64, 64, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn2): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act2): ReLU(inplace=True) ) ) (layer2): Sequential( (0): BasicBlock( (conv1): Conv2d(64, 128, kernel_size=(3, 3), stride=(2, 2), padding=(1, 1), bias=False) (bn1): BatchNorm2d(128, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act1): ReLU(inplace=True) (conv2): Conv2d(128, 128, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn2): BatchNorm2d(128, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act2): ReLU(inplace=True) (downsample): Sequential( (0): Conv2d(64, 128, kernel_size=(1, 1), stride=(2, 2), bias=False) (1): BatchNorm2d(128, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) ) ) (1): BasicBlock( (conv1): Conv2d(128, 128, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn1): BatchNorm2d(128, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act1): ReLU(inplace=True) (conv2): Conv2d(128, 128, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn2): BatchNorm2d(128, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act2): ReLU(inplace=True) ) ) (layer3): Sequential( (0): BasicBlock( (conv1): Conv2d(128, 256, kernel_size=(3, 3), stride=(2, 2), padding=(1, 1), bias=False) (bn1): BatchNorm2d(256, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act1): ReLU(inplace=True) (conv2): Conv2d(256, 256, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn2): BatchNorm2d(256, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act2): ReLU(inplace=True) (downsample): Sequential( (0): Conv2d(128, 256, kernel_size=(1, 1), stride=(2, 2), bias=False) (1): BatchNorm2d(256, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) ) ) (1): BasicBlock( (conv1): Conv2d(256, 256, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn1): BatchNorm2d(256, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act1): ReLU(inplace=True) (conv2): Conv2d(256, 256, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn2): BatchNorm2d(256, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act2): ReLU(inplace=True) ) ) (layer4): Sequential( (0): BasicBlock( (conv1): Conv2d(256, 512, kernel_size=(3, 3), stride=(2, 2), padding=(1, 1), bias=False) (bn1): BatchNorm2d(512, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act1): ReLU(inplace=True) (conv2): Conv2d(512, 512, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn2): BatchNorm2d(512, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act2): ReLU(inplace=True) (downsample): Sequential( (0): Conv2d(256, 512, kernel_size=(1, 1), stride=(2, 2), bias=False) (1): BatchNorm2d(512, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) ) ) (1): BasicBlock( (conv1): Conv2d(512, 512, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn1): BatchNorm2d(512, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act1): ReLU(inplace=True) (conv2): Conv2d(512, 512, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn2): BatchNorm2d(512, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act2): ReLU(inplace=True) ) ) (global_pool): SelectAdaptivePool2d (pool_type=avg, flatten=Flatten(start_dim=1, end_dim=-1)) (fc): Linear(in_features=512, out_features=1000, bias=True) ) 括弧で括られた部分がそのレイヤーの名前です。例えば(fc)となっている最終レイヤーにアクセスするには以下のようにします。 print(model.fc) Linear(in_features=512, out_features=1000, bias=True) ここで、出力の特徴数を変更するには以下です。 model.fc.out_features = 100 print(model.fc) Linear(in_features=512, out_features=100, bias=True) 同様にレイヤーを置き換えるのは簡単です。 print('before:', model.act1) model.act1 = torch.nn.LeakyReLU(inplace=True) print('after:', model.act1) before: ReLU(inplace=True) after: LeakyReLU(negative_slope=0.01, inplace=True) それではあるレイヤーを挿入したい場合はどうすればいいでしょうか。レイヤーを挿入するメソッドを見つけられなかったので、torch.nn.Sequentialを使って、以下のようにしました。以下はfc層の直前にDropoutを挿入する例です。 from torch.nn import Sequential, Dropout model.fc = Sequential(Dropout(0.5), model.fc) print(model) ResNet( (conv1): Conv2d(3, 64, kernel_size=(7, 7), stride=(2, 2), padding=(3, 3), bias=False) (bn1): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act1): LeakyReLU(negative_slope=0.01, inplace=True) (maxpool): MaxPool2d(kernel_size=3, stride=2, padding=1, dilation=1, ceil_mode=False) (layer1): Sequential( (0): BasicBlock( (conv1): Conv2d(64, 64, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn1): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act1): ReLU(inplace=True) (conv2): Conv2d(64, 64, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn2): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act2): ReLU(inplace=True) ) (1): BasicBlock( (conv1): Conv2d(64, 64, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn1): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act1): ReLU(inplace=True) (conv2): Conv2d(64, 64, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn2): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act2): ReLU(inplace=True) ) ) (layer2): Sequential( (0): BasicBlock( (conv1): Conv2d(64, 128, kernel_size=(3, 3), stride=(2, 2), padding=(1, 1), bias=False) (bn1): BatchNorm2d(128, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act1): ReLU(inplace=True) (conv2): Conv2d(128, 128, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn2): BatchNorm2d(128, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act2): ReLU(inplace=True) (downsample): Sequential( (0): Conv2d(64, 128, kernel_size=(1, 1), stride=(2, 2), bias=False) (1): BatchNorm2d(128, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) ) ) (1): BasicBlock( (conv1): Conv2d(128, 128, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn1): BatchNorm2d(128, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act1): ReLU(inplace=True) (conv2): Conv2d(128, 128, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn2): BatchNorm2d(128, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act2): ReLU(inplace=True) ) ) (layer3): Sequential( (0): BasicBlock( (conv1): Conv2d(128, 256, kernel_size=(3, 3), stride=(2, 2), padding=(1, 1), bias=False) (bn1): BatchNorm2d(256, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act1): ReLU(inplace=True) (conv2): Conv2d(256, 256, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn2): BatchNorm2d(256, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act2): ReLU(inplace=True) (downsample): Sequential( (0): Conv2d(128, 256, kernel_size=(1, 1), stride=(2, 2), bias=False) (1): BatchNorm2d(256, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) ) ) (1): BasicBlock( (conv1): Conv2d(256, 256, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn1): BatchNorm2d(256, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act1): ReLU(inplace=True) (conv2): Conv2d(256, 256, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn2): BatchNorm2d(256, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act2): ReLU(inplace=True) ) ) (layer4): Sequential( (0): BasicBlock( (conv1): Conv2d(256, 512, kernel_size=(3, 3), stride=(2, 2), padding=(1, 1), bias=False) (bn1): BatchNorm2d(512, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act1): ReLU(inplace=True) (conv2): Conv2d(512, 512, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn2): BatchNorm2d(512, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act2): ReLU(inplace=True) (downsample): Sequential( (0): Conv2d(256, 512, kernel_size=(1, 1), stride=(2, 2), bias=False) (1): BatchNorm2d(512, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) ) ) (1): BasicBlock( (conv1): Conv2d(512, 512, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn1): BatchNorm2d(512, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act1): ReLU(inplace=True) (conv2): Conv2d(512, 512, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn2): BatchNorm2d(512, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act2): ReLU(inplace=True) ) ) (global_pool): SelectAdaptivePool2d (pool_type=avg, flatten=Flatten(start_dim=1, end_dim=-1)) (fc): Sequential( (0): Dropout(p=0.5, inplace=False) (1): Linear(in_features=512, out_features=100, bias=True) ) ) 上手く挿入されました。 それではレイヤを削除したい場合はどうすればいいでしょうか。例えば最終層のクラス分類を外して、直前のレイヤーの出力を利用したいような場合です。この場合も、レイヤーを削除するメソッドが見つけられなかったので、以下のようにnn.Identityで置き換えます。 model.fc = torch.nn.Identity() print(model) ResNet( (conv1): Conv2d(3, 64, kernel_size=(7, 7), stride=(2, 2), padding=(3, 3), bias=False) (bn1): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act1): LeakyReLU(negative_slope=0.01, inplace=True) (maxpool): MaxPool2d(kernel_size=3, stride=2, padding=1, dilation=1, ceil_mode=False) (layer1): Sequential( (0): BasicBlock( (conv1): Conv2d(64, 64, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn1): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act1): ReLU(inplace=True) (conv2): Conv2d(64, 64, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn2): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act2): ReLU(inplace=True) ) (1): BasicBlock( (conv1): Conv2d(64, 64, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn1): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act1): ReLU(inplace=True) (conv2): Conv2d(64, 64, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn2): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act2): ReLU(inplace=True) ) ) (layer2): Sequential( (0): BasicBlock( (conv1): Conv2d(64, 128, kernel_size=(3, 3), stride=(2, 2), padding=(1, 1), bias=False) (bn1): BatchNorm2d(128, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act1): ReLU(inplace=True) (conv2): Conv2d(128, 128, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn2): BatchNorm2d(128, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act2): ReLU(inplace=True) (downsample): Sequential( (0): Conv2d(64, 128, kernel_size=(1, 1), stride=(2, 2), bias=False) (1): BatchNorm2d(128, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) ) ) (1): BasicBlock( (conv1): Conv2d(128, 128, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn1): BatchNorm2d(128, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act1): ReLU(inplace=True) (conv2): Conv2d(128, 128, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn2): BatchNorm2d(128, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act2): ReLU(inplace=True) ) ) (layer3): Sequential( (0): BasicBlock( (conv1): Conv2d(128, 256, kernel_size=(3, 3), stride=(2, 2), padding=(1, 1), bias=False) (bn1): BatchNorm2d(256, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act1): ReLU(inplace=True) (conv2): Conv2d(256, 256, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn2): BatchNorm2d(256, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act2): ReLU(inplace=True) (downsample): Sequential( (0): Conv2d(128, 256, kernel_size=(1, 1), stride=(2, 2), bias=False) (1): BatchNorm2d(256, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) ) ) (1): BasicBlock( (conv1): Conv2d(256, 256, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn1): BatchNorm2d(256, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act1): ReLU(inplace=True) (conv2): Conv2d(256, 256, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn2): BatchNorm2d(256, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act2): ReLU(inplace=True) ) ) (layer4): Sequential( (0): BasicBlock( (conv1): Conv2d(256, 512, kernel_size=(3, 3), stride=(2, 2), padding=(1, 1), bias=False) (bn1): BatchNorm2d(512, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act1): ReLU(inplace=True) (conv2): Conv2d(512, 512, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn2): BatchNorm2d(512, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act2): ReLU(inplace=True) (downsample): Sequential( (0): Conv2d(256, 512, kernel_size=(1, 1), stride=(2, 2), bias=False) (1): BatchNorm2d(512, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) ) ) (1): BasicBlock( (conv1): Conv2d(512, 512, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn1): BatchNorm2d(512, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act1): ReLU(inplace=True) (conv2): Conv2d(512, 512, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn2): BatchNorm2d(512, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act2): ReLU(inplace=True) ) ) (global_pool): SelectAdaptivePool2d (pool_type=avg, flatten=Flatten(start_dim=1, end_dim=-1)) (fc): Identity() ) これで一つ前の(global_pool)の出力をそのまま出力するようにできたので、クラス分類層を削除したのと同様になりました。 ところで、モデルの構造をよく見ると、(0)のように示されているものがあります。これはmodel.layer1.0などではアクセスできませんので、以下のようにします。 model.layer1[0] BasicBlock( (conv1): Conv2d(64, 64, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn1): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act1): ReLU(inplace=True) (conv2): Conv2d(64, 64, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (bn2): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (act2): ReLU(inplace=True) ) 無事にアクセスできました。これでどんな層も自由に編集できます。 以上です。

目的 指定のサイトで使用頻度が高い単語をワードクラウドを使って視覚化したい コード import requests from bs4 import BeautifulSoup import MeCab from wordcloud import WordCloud # ワードクラウドの保存ファイル名 file_name = "hage.png" # ワードクラウド化の対象サイトのURL url = 'https://ja.wikipedia.org/wiki/%E3%83%8F%E3%82%B2' # 対象サイトのURLのページを取得して解析 html = requests.get(url) soup = BeautifulSoup(html.content, 'html.parser') # scriptとstyleを要素から除去 for element in soup(["script", "style"]): element.decompose() # テキストを抽出 text = soup.get_text() # 行ごとに分けて両端の空白文字を削除 lines = [line.strip() for line in text.splitlines()] # 改行文字を結合 text = "\n".join(line for line in lines if line) # 形態素解析する mecab = MeCab.Tagger("-Ochasen") node = mecab.parseToNode(text) word_list = [] # 品詞を使って単語を分ける while node: # 品詞 part = node.feature.split(",")[0] # 配列の追加する単語を品詞で指定する if part in ["動詞","副詞","形容詞","名詞"]: word = node.surface # 単語 word_list.append(word) node = node.next words = ' '.join(word_list) # ワードクラウドで使用するフォント（メイリオ）のパス fpath = "C:/Windows/Fonts/meiryo.ttc" # ワードクラウドの作成 wordcloud = WordCloud(background_color="white", font_path=fpath, width=800, height=600).generate(words) # 指定の名前で保存する wordcloud.to_file("./"+file_name)

djangoのアクセス制限に関してです。 ログインしているユーザーのみといったアクセス制限の方法は日本語資料でも確認できましたが、 superuserのみといったアクセス制限の方法については、あまり見かけなかったので解説します。 誰かの役に立てば、幸いです。 環境 Django 2.2.4 関数で開発している場合 デコレータをつけることでいけます。 views.py from django.contrib.auth.decorators import permission_required from django.http import HttpResponse @permission_required('admin.can_add_log_entry') def test(request): return HttpResponse("Only access superuser!") classで開発している場合 せっかく、djangoを使っているなら、classベースで開発したいですよね。 そんな時のためのTips。 1. mixins.pyを定義 djangoのアプリ内で、mixins.pyを作成します。 mixins.py from django.contrib.auth.mixins import UserPassesTestMixin class SuperuserRequiredMixin(UserPassesTestMixin): def test_func(self): return self.request.user.is_superuser 2.views.pyで継承 ここで、定義したclassをviews.pyで使用します。 views.py from .mixins import SuperuserRequiredMixin class TestView(SuperuserRequiredMixin,TemplateView): template_name = "test/home.html" アクセス制限ができれば、 こんな感じになるはず。 参考になったら、幸いです。 参考文献 How to limit a view to superuser only?

はじめに Webアプリケーションでブラウザから直接S3にマルチパートファイルアップロードをする方法を備忘録的に残します。苦しんでる人の助けになればと思います。 動作する全てのコードはここにあります。 https://github.com/inunekousapon/s3_sts_upload_sample やること STSを利用したS3への直接アップロード マルチパートアップロード やらないこと マルチパートアップロードの自前実装 署名付きURLを利用したオブジェクトのアップロード ブラウザからS3に直接アップロードする方法 2つの方法があります。 署名付きURLアップロード AWS Security Token Serviceを使った方法 通常、署名付きアップロードを使うのが最適な選択肢になると思います。しかし、マルチパートアップロードを選択しようと思うと実装が難しいことがわかりました。 下記は署名付きURLを利用したマルチパートアップロードの記事です。 https://zenn.dev/kujime/articles/484b6e25d058fb このSDKのissueでは署名付きURLを利用したマルチパートアップロードが使いづらいという点について話し合われています。（2021.08.06現在解決していません） https://github.com/aws/aws-sdk-js/issues/1603 一方、AWSSDK for JavascriptにはManagedUploadというメソッドがあり、完全にマルチパートアップロードに対応しています。なぜ、署名付きURLアップロードに使えないのか理解に苦しみます。 The managed uploader allows for easy and efficient uploading of buffers, blobs, or streams, using a configurable amount of concurrency to perform multipart uploads where possible. This abstraction also enables uploading streams of unknown size due to the use of multipart uploads. マルチパートアップロードなどという面倒な実装はAWS公式のライブラリに任せたいので、「AWS Security Token Serviceを使った方法」を紹介したいと思います。 コード サーバー側 app.py import os import json import boto3 from flask import Flask, render_template app = Flask(__name__) @app.route('/', methods=['GET']) def upload_file(): sts_connection = boto3.client('sts', aws_access_key_id=os.getenv('AWS_ACCESS_KEY'), aws_secret_access_key=os.getenv('AWS_SECRET_KEY'), region_name='ap-northeast-1', ) b = sts_connection.get_federation_token( Name="Bob", DurationSeconds=900, Policy='{"Version":"2012-10-17","Statement":[{"Sid":"Stmt1","Effect":"Allow","Action":"s3:PutObject","Resource":"*"}]}' ) access_key_id = b['Credentials']['AccessKeyId'] secret_access_key = b['Credentials']['SecretAccessKey'] session_token = b['Credentials']['SessionToken'] return render_template('index.html', access_key_id=access_key_id, secret_access_key=secret_access_key, session_token=session_token, bucket_name='your bucket name', object_key='your object key', ) flaskによるサンプルです。 ブラウザで利用できる一時的な認証情報を取得するにはAsumeRoleまたはGetFederationTokenのどちらかを利用するといいと思います。今回はGetFederationTokenを利用しています。 （問題がある場合のご指摘お待ちしております） get_federation_tokenのPolicyをより厳密に書くことで、アップロードしたいキーのみ許可するなど署名付きURLによるアップロードと同等のアクセス制御を作ることができます。 FederationTokenを発行するには下記のアクセス権限が発行元のユーザーやロールに必要なので事前に設定する必要があります。 { "Version": "2012-10-17", "Statement": [ { "Sid": "VisualEditor0", "Effect": "Allow", "Action": "sts:GetFederationToken", "Resource": "*" } ] } クライアント側 index.html <!doctype html> <html> <head></head> <body> <title>Upload new File</title> <h1>Upload new File</h1> <input id="upload" type="file"> <button id="addFile" onclick="addFile()"> アップロード </button> <div id="progress"></div> <script src="https://sdk.amazonaws.com/js/aws-sdk-2.961.0.min.js"></script> <script type="text/javascript"> function addFile() { var files = document.getElementById("upload").files; if(!files.length) { return alert("No file."); } var bucket = "{{ bucket_name }}"; var file = files[0]; var key = "{{ object_key }}"; AWS.config.update({region: 'ap-northeast-1'}); AWS.config.credentials = new AWS.Credentials( accessKeyId="{{access_key_id}}", secretAccessKey="{{ secret_access_key }}", sessionToken="{{ session_token　}}", ) var promise = new AWS.S3.ManagedUpload({ params: { Bucket: bucket, Key: key, Body: file } }).on('httpUploadProgress', function(evt) { var disp = "Uploaded :: " + parseInt((evt.loaded * 100) / evt.total)+'%'; console.log(disp); document.getElementById('progress').innerText = disp; }) .promise(); promise.then( function(data) { alert("Successfully uploaded file."); }, function(err) { return alert("There was an error uploading your file: ", err.message); } ); } </script> </body> </html> 自前でマルチパートアップロードを書くよりかなり記述量が少なくなったかと思います。 'httpUploadProgress'のイベントでアップロードの進捗状況を取得することができます。 まとめ 公式のライブラリにマルチパートアップロードする仕組みがあるから、それをなんとか使う方法を模索したほうが楽です。 セキュリティ的に盲点があるかもしれないため、問題点があるのならコメントがほしいです。 一緒にAWSで苦しんでいける仲間を募集しています。

動機 普段はunityでアプリを作っていますが、アプリをインストールしてくれない時代になってきたので webアプリも習得したいが、個人的な勉強なので安価な環境が欲しい 最終的には、mediapipe、opencv 使いたい　 GCP、IBMcloudでインスタンス停止し忘れた経験あるので、高額請求されるくらいなら止めて欲しい 目指す環境 スマホカメラで撮影し、サーバーで処理し画像を返す ubuntu18.14 flask SSL コンテンツ app.py from flask import Flask app = Flask(__name__) @app.route("/") def hello(): return "Hello World!" if __name__ == "__main__": app.run( port=8000, debug=True) 自宅のPCに環境を構築 ngrokで公開してみる install.sh ./ngrok authtoken <your_auth_token> ./ngrok http 80 run.sh cd app.pyがある場所 python3 app.py 遅すぎるのと、アクセス回数が限られていたので断念 vpsを借りてみる プリペイドが使え残金がなくなると停止するconohaを使う VPSレンタル代 530 円/月 ドメイン代　￥2000/年 SSL設定 caddy２という、SSLが設定できるサーバーを見つけたので 設定してみる　ネットで探すとほとんどV1で、情報がまざっていて大変でした install.sh sudo apt install -y debian-keyring debian-archive-keyring apt-transport-https curl -1sLf 'https://dl.cloudsmith.io/public/caddy/stable/gpg.key' | sudo tee /etc/apt/trusted.gpg.d/caddy-stable.asc curl -1sLf 'https://dl.cloudsmith.io/public/caddy/stable/debian.deb.txt' | sudo tee /etc/apt/sources.list.d/caddy-stable.list sudo apt update sudo apt install caddy sudo nano /etc/caddy/Caddyfile install.sh www.sample.com:443 tls mail@gmail.com reverse_proxy localhost:8000 encode zstd gzip run.sh sudo ufw allow 443 sudo ufw allow 80 sudo caddy run --config /etc/caddy/Caddyfile 確認 install.sh cd app.pyがある場所 python3 app.py ブラウザから　https://www.sample.com　で Hello World!　と表示され なんとか、できたので次回、画像認識をやっていきたい

目的 都内在住・在勤の者として、細かく区市町村ごとの陽性者数推移が知りたかったので開発しました。 同じ都内といっても、区市町村ごとに傾向が違っています。 【今回作ったもの】 https://nobukuni-hyakutake.github.io/covid19tokyo/ の下の方の区市町村名をクリックするとグラフが表示されます。例： 原則毎日更新して公開しています。 なお、同サイトの上方にある分布マップの解説はこちらの記事をご覧ください。 環境 言語：Python3.9 使用ライブラリ： データ前処理はPandasとNumpy、グラフ作成はPlotly 使用サービス： GitHub pages コード解説 東京都の公開しているCSVは、各区市町村の各日付に対する累計陽性者数です。ここからPandasを使って各区市町村の1日ごとや、7日間ごとの陽性者数を出しました。また、各区市町村の人口を付けています。全体をCSV出力したものはこちらです。 covid19tokyo_preprocessed.csv(抜粋) ,group_code,date,count_sum,count_1day,count_7days,count_14days,last7days_ratio,last_day,Unnamed: 0,area,label,ruby,population,pref,en,number,jp_order 0,131016,2020-04-22,22.0,3.0,6.0,13.0,0.8571428571428571,2021-08-05,0,特別区,千代田区,ちよだく,65942.0,東京都,Chiyoda,0,35 1,131016,2020-04-23,22.0,0.0,4.0,13.0,0.4444444444444444,2021-08-05,0,特別区,千代田区,ちよだく,65942.0,東京都,Chiyoda,0,35 2,131016,2020-04-24,23.0,1.0,4.0,13.0,0.4444444444444444,2021-08-05,0,特別区,千代田区,ちよだく,65942.0,東京都,Chiyoda,0,35 3,131016,2020-04-25,23.0,0.0,4.0,13.0,0.4444444444444444,2021-08-05,0,特別区,千代田区,ちよだく,65942.0,東京都,Chiyoda,0,35 その後の、推移グラフ作成についてを詳しく解説します。 62区市町村+東京都全体(以下、各自治体)で、合計63回グラフ作成をfor文で繰り返します。各自治体には、あらかじめ0~62のnumberを振ってあります。 main.py for i in range (63): iに対応した各自治体のデータだけを、データフレームから取り出します。人口で割って10万人あたりの数を出します。 なおステージ3基準、4基準の数字はそれぞれ１週間あたり15人、25人を７で割った数字です。 main.py dfgraph00201=out.loc[(out['number']==i),['group_code','label','ruby','date','count_1day','count_7days','population','en','number']].copy().reset_index() label00201=dfgraph00201['label'][0] en00201=dfgraph00201['en'][0] title00201=dfgraph00201['label'][0]+' 新型コロナウイルス陽性者数 スマートフォンは横向きにして下さい' dfgraph00201['sevendays_ave']=round((dfgraph00201['count_7days']/7)/dfgraph00201['population']*100000,1) dfgraph00201['count_1day_p']=round(dfgraph00201['count_1day']/dfgraph00201['population']*100000,1) dfgraph00201['stage4']=3.6 dfgraph00201['stage3']=2.1 １日ごとの陽性者数の棒グラフ部分を作ります。hover(マウスオーバーさせると出てくるテキスト)を表示させたくないのでhoverinfo = "none"とします。 main.py fig00202=go.Bar( x=dfgraph00201["date"], y=dfgraph00201["count_1day_p"], name='10万人あたり', marker={"color": "#99cc66"}, hoverinfo = "none", ) 7日間平均の陽性者数の折れ線グラフ部分を作ります。 main.py fig00203=go.Scatter( x=dfgraph00201["date"], y=dfgraph00201["sevendays_ave"], name='10万人あたり7日間平均', line={"color": "#cc6600","width":2}, ) ステージ３、４の破線を作ります。 main.py figstage4=go.Scatter( x=dfgraph00201["date"], y=dfgraph00201["stage4"], name='ステージ4基準', line={"width":2, "color": "red", "dash":"dash"}, hoverinfo = "none" ) figstage3=go.Scatter( x=dfgraph00201["date"], y=dfgraph00201["stage3"], name='ステージ3基準', line={"width":2, "color": "#ffcc00", "dash":"dash"}, hoverinfo = "none" ) フォントサイズ、グラフタイトルの設定をします。 main.py layout=go.Layout( font=dict(size=20), title={"text":title00201, }, このrangeselectorの設定によって、表示期間を全期間または直近３ヶ月間か、グラフを見る人がクリックで選べます。plotlyは、こういうインタラクティブなグラフが作れて便利です。 main.py xaxis={ "linecolor": "black", "rangeselector":{ "buttons":[ {"label":"全期間","step":"all"}, {"label":"直近3ヵ月間","step":"month","count":3,"stepmode":"backward"} ], }, "type":"date", }, y軸のタイトルと色、プロット部分全体の背景色を設定 main.py yaxis={ "title":{ "text": '10万人あたり', }, "linecolor": "black", }, hovermode='x', plot_bgcolor="#ffffff", ) ここまでで作ったグラフオブジェクトを全て、Figureとして登録します。 main.py fig04=go.Figure(data=[fig00202, fig00203, figstage4, figstage3], layout=layout) 軸、凡例の調整 main.py fig04.update_yaxes( rangemode="nonnegative" ) fig04.update_layout(legend_orientation="h") fig04.update_layout(legend={"x":0,"y":-0.18}) 最後の日にannotation(注釈のこと。上で出した例だと、一番直近の日に「21/8/4: 47.8」と書いてあるテキストボックスと矢印の部分)をつけます。 デフォルトだと見づらいので、ax=-80, ay=-80としてプロットの左上側に表示されるように設定しました。 main.py sevendays_ave_lastday=dfgraph00201.loc[(dfgraph00201['date']==last_day1),['sevendays_ave']].mean()[0] fig04.update_layout( annotations=[ go.layout.Annotation( x=last_day1, y=sevendays_ave_lastday, xref="x", yref="y", text=str(last_day1.year)[2:4]+"/"+str(last_day1.month)+"/"+str(last_day1.day)+": "+str(sevendays_ave_lastday), showarrow=True, arrowhead=1, bgcolor="#cc6600", font={"size":15,"color":"black"}, ax=-80, ay=-80, opacity=0.7, ) ] ) 余白の調整。また、x軸の表示形式が、デフォルトだと"Mar 2021"など英語なので日本式の21/3/1という形にしました。 main.py fig04.update_layout( margin={"t":120} ) fig04.update_yaxes(automargin=False) fig04.update_xaxes(type='date', tickformat="%y/%-m/%-d", tick0='2020-05-01', dtick="M2") 最後に、"[各自治体名]_g.html"というファイル名で出力します。例えば新宿区だったら"Shinjuku_g.html"です。 main.py fig04.write_html("docs/"+en00201+"_g.html") 以上を前述のfor文で繰り返すことで、各自治体63箇所のグラフの出力が一気にできます。 私のMacBook Air (CPU:1.1 GHz クアッドコアIntel Core i5) ではデータ取り込み・前処理・全てのグラフ出力までを15秒で処理できました。 コード全文 データ取り込み、前処理も含めた全文はこちらをご参照ください。 https://github.com/Nobukuni-Hyakutake/covid19tokyo 参考サイト 見栄えの調整で参考にさせていただきました。 軸の日付表示形式で参考にさせていただきました。 参考文献 Pythonインタラクティブ・データビジュアライゼーション入門 終わりに MITライセンスとして公開しておりますので、東京都以外の道府県で同様のグラフを作られる際はどうぞ参考にしてみて下さい。 【今回作ったもの】 https://nobukuni-hyakutake.github.io/covid19tokyo/ 下の方の区市町村名のリンクをクリックすると、グラフが表示されます。 【例：新宿区はこちらです】 https://nobukuni-hyakutake.github.io/covid19tokyo/Shinjuku_g.html 原則毎日、グラフを更新しています。 最後になりましたが感染された方へお見舞い申し上げます。

PyinstallerでWindowsで実行可能なexe化を行った際にPydriveがインポートされず、 GoogleDriveへのファイルのアップロードが失敗してしまったので、そのエラーの対処方法の覚書です。 環境 ・OS: Windows 10 64bit ・Python 3.8.10 ・PyInstaller 4.3 エラー箇所 エラーの原因となっていたのは、下記のライブラリをインポートした部分です。 やりたかったこととしては、ローカルファイルを共有ドライブにアップロードするでしたが、Uploadの際に、Anaconda上では問題なく動作していましたが、Pyinstallerでexe化するとアップロードの際にエラーが発生して、GoogleAPIErrorが発生してしまった。 GoogleDriveup.py from pydrive.auth import GoogleAuth from pydrive.drive import GoogleDrive ローカルファイルを共有ドライブにアップロードする方法は以下の記事を参照しました。 pydrive を使って 共有ドライブにアップロードする方法 エラー内容 コンソール上で確認できたエラーはgoogleapiclient.errors.UnknownApiNameOrVersionというものでした。もう一つのモジュールでGoogleDrive上のspreadsheetの書き込みが成功していたので、PyDriveの問題ということが分かりました。 対処方法 現状のPyinstallerはexe化実行時に、PyDriveを読み込み出来てずそのせいでエラーが起きていました。 対処方法としては、Pyinstaller実行後に出力されるspecファイルの中身を書き換える必要があります。 1.ファイルの先頭にfrom PyInstaller.utils.hooks import collect_data_filesを追加する。 2.data=[],の部分をdatas=collect_data_files("googleapiclient"),に書き換える。 GoogleDriveup.spec # -*- mode: python ; coding: utf-8 -*- from PyInstaller.utils.hooks import collect_data_files block_cipher = None a = Analysis(['EXE化したいファイル名'], pathex=['ファイルへのパス'], binaries=[], datas=collect_data_files("googleapiclient"), hiddenimports=[], hookspath=[], runtime_hooks=[], excludes=[], win_no_prefer_redirects=False, win_private_assemblies=False, cipher=block_cipher, noarchive=False) pyz = PYZ(a.pure, a.zipped_data, cipher=block_cipher) exe = EXE(pyz, a.scripts, [], exclude_binaries=True, name='Enabler_registor', debug=False, bootloader_ignore_signals=False, strip=False, upx=True, console=False ) coll = COLLECT(exe, a.binaries, a.zipfiles, a.datas, strip=False, upx=True, upx_exclude=[], name='EXE化したいファイル') これで問題なくファイルのアップロードが出来ました。

自動でMicrosoft teamsで個人あてにチャットを送る方法 　相手からの電子メッセージにできるだけ早く、かつ漏れなく確認することは、業務上大事である。 メールは全体へ連絡したり情報共有するのは使いやすいが、自身とあまり関係のないメッセージも大量に来るので漏れなく確認するのは難しい。チャットにはプッシュ通知機能があり、その連絡に気付きやすく効率的な会話ができる。今回は自動でteamsにメッセージを送る方法を紹介する。 install python selenium by pip pip installでseleniumをinstallする。 pip install selenium Web版teamsの操作 　今回はChrome web driverを使ってWeb版TeamsをChromeで開いて操作する。 　開くページのURLにチャットの相手のユーザーネームが書いてるので、複数送りたい相手がいる場合はその箇所を送りたい相手に変更する。 driver_path = r'/path/chromedriver.exe' driver = webdriver.Chrome(executable_path = driver_path) driver.get('https://teams.microsoft.com/dl/launcher/***' + Name +'**') driver.maximize_window() デスクトップアプリ版を使うかどうかアラートが出る場合があるので、例外処理でアラートへswichして"いいえ"を押す。 try: alert = driver.switch_to.alert text = alert.text alert.dismiss() except: pass メッセージボックスをクリックして入力可能状態にし、過去の未送信のメッセージがある場合があるのでBack Spaceで消して、送りたいメッセージを送る。 driver.find_element_by_xpath('message boxのxpath').click() for i in range(30): driver.find_element_by_xpath('message boxのxpath').send_keys(Keys.BACK_SPACE) message = """ Hello! jon How are you doing? Would you like to go fishing with me this weekend? """ driver.find_element_by_xpath('message boxのxpath').send_keys(message) driver.find_element_by_id('send-message-button').click() まとめ Microsoft teamsで個人あてにチャットを送る方法を紹介した。定期業務で忘れやすい人はこれでreminderを送るのもいいかもしれません。

アラビア数字から漢数字に変換する処理　備忘録 無量大数って有限だったんだ、、知らんかった d1 = {u'1': u'一', u'2': u'二', u'3': u'三', u'4': u'四', u'5': u'五', u'6': u'六', u'7': u'七', u'8': u'八', u'9': u'九'} d3 = {1: u'十', 2: u'百', 3: u'千'} d4 = {4: u'万', 8: u'億', 12: u'兆', 16: u'京', 20: u'垓', 24: u'予', 28: u'穣', 32: u'溝', 36: u'潤', 40: u'正', 44: u'載', 48: u'極', 52: u'恒河沙', 56: u'阿僧祇', 60: u'那由他', 64: u'不可思議', 68: u'無量大数'} def int2kanji(integer): sio = str(integer) i = len(sio) - 1 ret = u"" for s in sio: mod4 = i % 4 if (mod4 == 0): ret += d1.get(s, u"") ret += d4.get(i, u"") elif (s != u'0'): if s != u'1': ret += d1.get(s, s) ret += d3.get(mod4, u"") i -= 1 return ret print(int2kanji(110304567)) # -> 一億千三十万四千五百六十七

Logicaとは LogicaとはGoogleが発表した論理プログラミング言語で、BigQueryやPostgreSQLのクエリにコンパイルして実行できます。 日本語ではこれらの記事で紹介されており、だいたいの概要は分かると思います。 logica で readable, testable な SQL の夢を見る 論理プログラミング言語Logicaでデータサイエンス100本ノック 引用すると次のような特徴があり、私の所属するチームではBigQueryでそれなりに複雑な処理が実装されているので、特に「モジュール機能があること」「テストが実施できる?(ドキュメントにも明確に記載がないのでまだ詳細不明)」ことに期待しています。また、GitHubのDiscussions（PageRank Alogrithm）を見る限り再帰も（回数は限定されてそうですが）使えるらしく、もしかするとBigQuery側で柔軟にアルゴリズムを実装していくことも可能かもしれません。 最初に、プログラミング言語Logicaの特徴を纏めておく。 論理型プログラミング言語: このカテゴリではPrologが有名 SQLにコンパイルされる: 現状BigQueryとPostgreSQLに対応 モジュール機構がある: SQLと比較した強み コンパイラはPythonで書かれている: Jupyter NotebookやGoogle Colabですぐ始められる Colabでチュートリアルが用意されているので、まずこちらからやると良いと思う。 ただ、チュートリアルやドキュメントを見る限りlogicaのコマンドで実行するか、Collabで実行するかという方法しか紹介されていません。そのため、私は「複雑なクエリをPythonのバッチ処理で実行したい」と考え、Python自体でlogicaをimportして実行してみました。 （※ただ、後から思いついたのですが、Python内で利用できなくても、Dockerのビルド時にlogicaをSQLに変換した結果を保存して利用すればいいだけかもしれません） LogicaのプログラムをSQLに変換するPythonスクリプト チュートリアルにあるこちらのLogicaプログラムをPython内で変換します。 mentions.l @OrderBy(Mentions, "mentions desc"); @Limit(Mentions, 10); Mentions(person:, mentions? += 1) distinct :- gdelt-bq.gdeltv2.gkg(persons:, date:), Substr(ToString(date), 0, 4) == "2020", the_persons == Split(persons, ";"), person in the_persons; 以下のようなスクリプトで実行できます。logica.pyの内容を追って、不要な処理を削除したものです。 main.py from logica.parser_py import parse from logica.compiler import universe with open("./mentions.l") as f: program = f.read() # 構文木(?)に変換する parsed_rules = parse.ParseFile(program)["rule"] # SQLに変換する program = universe.LogicaProgram(parsed_rules) sql = program.FormattedPredicateSql("Mentions") # 呼び出したい変数(Mentions)を指定する print(sql) 結果、以下のSQLが出力されます。このSQL文字列をBigQueryのクライアントで利用すればokです。 SELECT x_3 AS person, SUM(1) AS mentions FROM gdelt-bq.gdeltv2.gkg AS gdeltbq_gdeltv2_gkg, UNNEST(SPLIT(gdeltbq_gdeltv2_gkg.persons, ";")) as x_3 WHERE (SUBSTR(CAST(gdeltbq_gdeltv2_gkg.date AS STRING), 0, 4) = "2020") GROUP BY person ORDER BY mentions desc LIMIT 10; 参考資料 Logica Tutorial GitHubリポジトリ

前回の記事では非適応戦略では28bitあたり14回のチェックで同定する手法を示した。 これ以上は何故か一意に解くことが出来なくなったため、以降詳しくは書いていない。 とはいえ問題が発生する確率は当初の見積もりよりも非常に低い確率であったので「高確率で同定可能な手法」という但し書きが可能ならもっと28bitあたり14回よりも高効率で同定が可能である。 一意に解ける係数は求められていないが、高確率では解けるのでその係数を参考にメモしておく。 解11 15/32*N回 4bitのチェック関数である$bit4(a,b,c,d)$を仮定して立ててみた。今までの傾向から1個のbitは値が1回だけ、残りの3bitは15回中8回が値が1にならなければならない。この傾向は以下の手順から推測される。 例えば12bitを7回のチェックで同定する手法では7回のチェックの合計値(赤枠の総和)は4,4,2,4,4,2,4,2,4,1,4,4となっている。 これを2で割った商が2,2,1,2,2,1,2,1,2,0,2,2、2で割った余りが0,0,0,0,0,0,0,0,0,1,0,0である。ここで10bit目が同定される。 次に緑枠の合計値は2,2,1,2,2,1,2,2,0,1,2,2なので前述の商と余りの和との差は0,0,0,0,0,0,0,-1,2,0,0,0となる。この値の2で割った時の商と余りから8bit目と9bit目が同定できる。 同様に青枠の合計値との計算で6bit目と7bit目が同定できる。 また同様に4bit目の合計が0になるように合計枠を選べば3bit目と4bit目が同定できる。 次に紫枠(5bit目の合計が0になるように合計枠)の合計値を考えれば今、6bit目と8bit目は同定されているから差分は0,0,0,0,2,0,0,0,0,0,1,0である。これの2で割った時の商と余りから5bit目と11bit目が同定される。 また同様に1bit目の合計が0になるように合計枠を選べば1bit目が同定できる。 また同様に2bit目の合計が0になるように合計枠を選べば2bit目が同定できる。 このように見ていけば$bit2(a,b)$の3回チェック時点の合計値は$(1,2)$ $bit3(a,b,c)$の7回チェック時点の合計値は$(1,4,4)$だから $bit4(a,b,c,d)$の15回チェック時点の合計値は$(1,8,8,8)$となることが推測される。 import numpy as np from sympy import * def test(y_pred, y_true): return len(y_pred) - np.sum(np.abs(y_pred - y_true)) def bit2(a,b): value = int(a + b) r = np.array(list(format(value, '02b'))).astype('int64') return r def bit3(a,b,c): value = int(a + b + c) if (a,b,c)==(1,1,1): value += 4 r = np.array(list(format(value, '03b'))).astype('int64') return r def bit4(a,b,c,d): r1,r2,r3 = [3, 5, 6, 8, 9,10,12,15], [1, 2, 3, 7,12,13,14,15], [1, 2, 4, 7,11,12,14,15] value = 0 r0 = int(a + 2 * b + 4 * c + 8 * d) if r0==15: value += 8 if r0 in r1: value += 4 if r0 in r2: value += 2 if r0 in r3: value += 1 r = np.array(list(format(value, '04b'))).astype('int64') return r def Solve11(length): global iter m = 15 Am = 32 M = np.zeros((Am,m)) M_list = [[1], [2], [1,2], [4], [1,4], [2,4], [1,2,4], [8], [1,8], [2,8], [1,2,8], [4,8], [1,4,8], [2,4,8], [1,2,4,8]] index_list = [] k = 0 for i in range(m): index_list.append(k) k += len(M_list[i]) k = 0 for i in range(m): index = index_list[i] if len(M_list[i])==1: for j in range(m): M[index,j] = np.array(list(format(j+1, '06b'))).astype('int64')[-1-k] k += 1 if len(M_list[i])==2: index2 = index_list[M_list[i][0]-1] index3 = index_list[M_list[i][1]-1] for j in range(m): M[index:index+2,j] = bit2(M[index2,j], M[index3,j]) if len(M_list[i])==3: index2 = index_list[M_list[i][0]-1] index3 = index_list[M_list[i][1]-1] index4 = index_list[M_list[i][2]-1] for j in range(m): M[index:index+3,j] = bit3(M[index2,j], M[index3,j], M[index4,j]) if len(M_list[i])==4: index2 = index_list[M_list[i][0]-1] index3 = index_list[M_list[i][1]-1] index4 = index_list[M_list[i][2]-1] index5 = index_list[M_list[i][3]-1] for j in range(m): M[index:index+4,j] = bit4(M[index2,j], M[index3,j], M[index4,j], M[index5,j]) if len(M_list[i])==5: index2 = index_list[M_list[i][0]-1] index3 = index_list[M_list[i][1]-1] index4 = index_list[M_list[i][2]-1] index5 = index_list[M_list[i][3]-1] index6 = index_list[M_list[i][4]-1] for j in range(m): M[index:index+5,j] = bit5(M[index2,j], M[index3,j], M[index4,j], M[index5,j], M[index6,j]) y_check = np.array([1] * int(Am)) result = np.array([0]*m) x = [Symbol('x%d' % (i)) for i in range(Am)] for index in range(0,length//Am*Am,Am): y_true2 = y_true[index:index+Am] for j in range(m): result[j] = test(y_check[:np.count_nonzero(M[:,j])], y_true2[M[:,j]==1]) iter += 1 eq = [] for i in range(Am): eq.append(x[i] * x[i] - x[i]) for j in range(m): eq0 = 0 for i in range(Am): eq0 += M[i,j] * x[i] eq.append(eq0 - result[j]) y_pred[index:index+Am] = np.array(solve(eq, x)[0]) print(index, test(y_pred[index:index+Am], y_true[index:index+Am])) n = 1024*16 iter = 0 y_pred = np.array([0] * n) y_true = np.random.randint(0, 2, (n)) Solve11(n) print(iter, test(y_pred, y_true)) ------------------------------ ... 6688 32 6720 32 6752 32 6784 18 <= failure!!! 6816 32 6848 32 適当に$bit4(a,b,c,d)$関数を立てて計算させてみるが、これを実行した場合に低確率で失敗する。とはいえ16384bitで1回現れる程度なので確率的には低い。勿論、出ない時もある。 これは$bit4(a,b,c,d)$関数が一意に求めるには正しくないのかもしれない。(とはいえ自分の方でも適当な乱数で少しは探したのでこの値は適当に決めた係数の中ではましな方である) 解12 30/75*N回 その次に$bit5(a,b,c,d,e)$関数が出てくるのは$m=31,Am=80$である。ということは$bit4(a,b,c,d)$関数で書けるのはその一つ前までだから$m=30,Am=75$で求める。これは75bitあたり30回のチェックで同定できるから効率的には0.4である。 def Solve12(length): global iter m = 30 Am = 75 M = np.zeros((Am,m)) M_list = [[1], [2], [1,2], [4], [1,4], [2,4], [1,2,4], [8], [1,8], [2,8], [1,2,8], [4,8], [1,4,8], [2,4,8], [1,2,4,8], [16], [1,16], [2,16], [1,2,16], [4,16], [1,4,16], [2,4,16], [1,2,4,16], [8,16], [1,8,16], [2,8,16], [1,2,8,16], [4,8,16], [1,4,8,16], [2,4,8,16]] ------------------------------ ... 1500 75 1575 75 1650 62 <= failure!!! 1725 75 ... 10500 75 10575 62 <= failure!!! 10650 75 10725 75 10800 55 <= failure!!! 10875 75 ... 13875 75 13950 60 <= failure!!! 14025 75 ... 15150 75 15225 60 <= failure!!! 15300 75 失敗の件数は16384bit(75bitは218回)で5回程度で解11よりは高い。 とはいえ98%は成功するのでテストデータが1000個程度ならほとんど失敗しないだろう。 解13 31/80*N回 $bit5(a,b,c,d,e)$の31回チェック時点の合計値は$(1,16,16,16,16)$となることが推測される。 この数字を合計16となるように決めて実行した所、以下の様な$bit5$関数を定義した。 def bit5(a,b,c,d,e): r01 = [2, 3, 5,10,13,14,15,16,19,20,21,22,27,29,30,31] r02 = [1, 2, 4, 6, 9,10,11,13,15,17,19,21,26,29,30,31] r03 = [2, 3, 4, 7, 8,11,12,13,14,16,18,22,24,28,30,31] r04 = [2, 5, 6, 7, 9,10,13,14,16,19,23,24,25,26,28,31] value = 0 r0 = int(a+2*b+4*c+8*d+16*e) if r0==31: value += 16 if r0 in r01: value += 8 if r0 in r02: value += 4 if r0 in r03: value += 2 if r0 in r04: value += 1 r = np.array(list(format(value, '05b'))).astype('int64') return r def Solve13(length): global iter m = 31 Am = 80 M = np.zeros((Am,m)) M_list = [[1], [2], [1,2], [4], [1,4], [2,4], [1,2,4], [8], [1,8], [2,8], [1,2,8], [4,8], [1,4,8], [2,4,8], [1,2,4,8], [16], [1,16], [2,16], [1,2,16], [4,16], [1,4,16], [2,4,16], [1,2,4,16], [8,16], [1,8,16], [2,8,16], [1,2,8,16], [4,8,16], [1,4,8,16], [2,4,8,16], [1,2,4,8,16]] 失敗の件数は1024*16bit(80bitは204回)で3回程度で解13とさほど変わらない。というか解13よりも少ない。 解14 62/186*N回 $bit5(a,b,c,d,e)$関数で作れるのは$bit6(a,b,c,d,e,f)$関数が出てくる$m=63,Am=192$の一歩前、つまり$m=62,Am=186$までである。効率は0.333でテストデータの$1/3$回のチェックで同定が出来る。これは理論上限が$α*N/log_{2}N$で$2≦α≦log_{2}9(3.17)$より$N=1024$の場合、$200～320$と見積もられるからその値に近い。 しかし、問題は計算に非常に時間が掛かることである。 def Solve14(length): global iter m = 62 Am = 186 M = np.zeros((Am,m)) M_list = [[1], [2], [1,2], [4], [1,4], [2,4], [1,2,4], [8], [1,8], [2,8], [1,2,8], [4,8], [1,4,8], [2,4,8], [1,2,4,8], [16], [1,16], [2,16], [1,2,16], [4,16], [1,4,16], [2,4,16], [1,2,4,16], [8,16], [1,8,16], [2,8,16], [1,2,8,16], [4,8,16], [1,4,8,16], [2,4,8,16], [1,2,4,8,16], [32], [1,32], [2,32], [1,2,32], [4,32], [1,4,32], [2,4,32], [1,2,4,32], [8,32], [1,8,32], [2,8,32], [1,2,8,32], [4,8,32], [1,4,8,32], [2,4,8,32], [1,2,4,8,32], [16,32],[1,16,32],[2,16,32],[1,2,16,32],[4,16,32],[1,4,16,32],[2,4,16,32],[1,2,4,16,32],[8,16,32],[1,8,16,32],[2,8,16,32],[1,2,8,16,32],[4,8,16,32],[1,4,8,16,32],[2,4,8,16,32]] 時間比較 実行時間の比較を行ってみる。16384bitの長さのチェックを$m=14,15,30,31,47,62$の時間と失敗数をメモする。 m = 14 Am = 28 M = np.zeros((Am,m)) M_list = [[1], [2], [1,2], [4], [1,4], [2,4], [1,2,4], [8], [1,8], [2,8], [1,2,8], [4,8], [1,4,8], [2,4,8]] # m = 15 # Am = 32 # M = np.zeros((Am,m)) # M_list = [[1], [2], [1,2], [4], [1,4], [2,4], [1,2,4], [8], [1,8], [2,8], [1,2,8], [4,8], [1,4,8], [2,4,8], [1,2,4,8]] # m = 30 # Am = 75 # M = np.zeros((Am,m)) # M_list = [[1], [2], [1,2], [4], [1,4], [2,4], [1,2,4], [8], [1,8], [2,8], [1,2,8], [4,8], [1,4,8], [2,4,8], [1,2,4,8], # [16], [1,16], [2,16], [1,2,16], [4,16], [1,4,16], [2,4,16], [1,2,4,16], [8,16], [1,8,16], [2,8,16], [1,2,8,16], [4,8,16], [1,4,8,16], [2,4,8,16]] # m = 31 # Am = 80 # M = np.zeros((Am,m)) # M_list = [[1], [2], [1,2], [4], [1,4], [2,4], [1,2,4], [8], [1,8], [2,8], [1,2,8], [4,8], [1,4,8], [2,4,8], [1,2,4,8], # [16], [1,16], [2,16], [1,2,16], [4,16], [1,4,16], [2,4,16], [1,2,4,16], [8,16], [1,8,16], [2,8,16], [1,2,8,16], [4,8,16], [1,4,8,16], [2,4,8,16], [1,2,4,8,16]] # m = 47 # Am = 128 # M = np.zeros((Am,m)) # M_list = [[1], [2], [1,2], [4], [1,4], [2,4], [1,2,4], [8], [1,8], [2,8], [1,2,8], [4,8], [1,4,8], [2,4,8], [1,2,4,8], # [16], [1,16], [2,16], [1,2,16], [4,16], [1,4,16], [2,4,16], [1,2,4,16], [8,16], [1,8,16], [2,8,16], [1,2,8,16], [4,8,16], [1,4,8,16], [2,4,8,16], [1,2,4,8,16], # [32], [1,32], [2,32], [1,2,32], [4,32], [1,4,32], [2,4,32], [1,2,4,32], [8,32], [1,8,32], [2,8,32], [1,2,8,32], [4,8,32], [1,4,8,32], [2,4,8,32], [1,2,4,8,32]] # m = 62 # Am = 186 # M = np.zeros((Am,m)) # M_list = [[1], [2], [1,2], [4], [1,4], [2,4], [1,2,4], [8], [1,8], [2,8], [1,2,8], [4,8], [1,4,8], [2,4,8], [1,2,4,8], # [16], [1,16], [2,16], [1,2,16], [4,16], [1,4,16], [2,4,16], [1,2,4,16], [8,16], [1,8,16], [2,8,16], [1,2,8,16], [4,8,16], [1,4,8,16], [2,4,8,16], [1,2,4,8,16], # [32], [1,32], [2,32], [1,2,32], [4,32], [1,4,32], [2,4,32], [1,2,4,32], [8,32], [1,8,32], [2,8,32], [1,2,8,32], [4,8,32], [1,4,8,32], [2,4,8,32], [1,2,4,8,32], # [16,32],[1,16,32],[2,16,32],[1,2,16,32],[4,16,32],[1,4,16,32],[2,4,16,32],[1,2,4,16,32],[8,16,32],[1,8,16,32],[2,8,16,32],[1,2,8,16,32],[4,8,16,32],[1,4,8,16,32],[2,4,8,16,32]] m 14 15 30 31 47 62 Am 28 32 75 80 128 186 効率(m/Am) 0.500 0.469 0.400 0.388 0.367 0.333 実行数(16384bit/Am) 585 512 218 204 128 88 失敗数 0 1 5 3 5 - 実行時間 730[s] 986[s] 4883[s] 6255[s] 18709[s] - 実行時間/実行数 1.3[s] 1.9[s] 22.4[s] 30.7[s] 146[s] 1117[s] $m=62$では時間が掛かりすぎて終わらなかった。何故か途中で止まってしまっている。 演算時間は$m=62$以外では$Am^3$にほぼ比例している。その領域の関係式でいえば$m=62$の実行時間/実行数は$400[s]$程度が妥当である。 $m=15$以降ではやはり低確率で失敗してしまうので$bit4,bit5$関数の係数の決め方が良くないと思われる。 とはいえTaitanic問題はテストデータが418件と数が多くないのでこの二値分類コンペではこの不完全な手法でも確率的には数回の試行中で失敗する確率は低いと思われる。 まとめ 低確率で失敗してしまうがサブミット回数はテストデータの約$1/3$(186bitあたり62回のチェック)で済む手法を示した。また、参考までに同様に次の$bit6$関数を定義出来た場合、理論的には$m=126,Am=441$まで定義できる。 Taitanic問題のテストデータ418件を探すとちょうど$m=122,Am=418$であるから418bitの結果は122回のサブミット結果のスコアから答えのラベルを同定できる。効率は0.292で3割を下回る。 しかし、現実的には$m=122$でこの連立方程式を解くのには演算時間が$Am^3$に比例すると仮定しても$4500[s]$かかり、また、失敗する確率も考えると今回の手法では現実的ではない。 また、Kaggleのスコアは少数5桁まで出るのでテストデータが1000件未満のサブミット時の正解bitの数を求めるのは容易である。例えば上記なら正解数は$158/418$である。

問題 5つの品物の価格と価値が以下の表に示されるようになっているとする． 価格の合計が 15以下で，価値の合計が 19以上になる場合はあるか? 品物 価格 価値 品物0 3 4 品物1 4 5 品物2 5 6 品物3 7 8 品物4 9 10 回答 example_knapsack.py from z3 import * price = [3, 4, 5, 7, 9] value = [4, 5, 6, 8, 10] X = [Int("X%s" % i) for i in range(5)] s = Solver() s.add([And(0 <= X[i], X[i] <= 1) for i in range(5)]) s.add(Sum([price[i]*X[i] for i in range(5)]) <= 15) s.add(Sum([value[i]*X[i] for i in range(5)]) >= 19) print(s.check()) if s.check() == sat: print(s.model()) 出力 価値の合計が19の場合、 unsat 価値の合計を18の場合、 sat [X1 = 0, X0 = 1, X4 = 0, X3 = 1, X2 = 1] 解説 s.add(Sum([value[i]*X[i] for i in range(5)]) >= 19) の19が価値の合計に対する制約 ここを変えることで、価値の制限による違いを見ることができる。 他の例題 Z3Py個人的ポータル へ

アーティストのアルバム一覧取得 SpotipyGetAlbumList.py # アーティストのアルバム一覧取得 import spotipy from spotipy.oauth2 import SpotifyClientCredentials client_id = "クライアントID" # Spotify for Devで作成したものを貼り付け client_secret = "シークレットID" # Spotify for Devで作成したものを貼り付け client_credentials_manager = spotipy.oauth2.SpotifyClientCredentials( client_id, client_secret) spotify = spotipy.Spotify( client_credentials_manager=client_credentials_manager) # アーティストのID artist = "spotify:artist:4rd6QvHo3CdZNmZrPJs8yz" # 例：人間椅子 # アルバム情報を取得して、アルバム名、リリース日、URLなどを表示する results = spotify.artist_albums(artist) name = [] release_date = [] # for item in results['items']: # if item['available_markets'] == ['JP']: # print(item['name'] + "/" + item['release_date']) # marketsに対するifの有無による結果の変化を見てみる for item in results['items']: print(item['name'] + "/" + item['release_date']) `

はじめに 自分自身の知識の整理も兼ねて、需要予測の手法についてPythonで実装・まとめていきたいと思います。 業務上、よく需要予測をしている職場もあると思いますが、ずっと昔から使われていてエクセルで実装できる移動平均のようなクラシックな手法以外にも、最近では機械学習を使った、より高度な手法もあります。 以下ではPythonを使って各種需要予測手法を実装・比較していきます。 需要予測はなぜ重要なのか サプライチェーンマネジメント（SCM）を考える上で、需要予測は欠かせません。 というのも、多期間における在庫モデルや生産計画を考える上で、需要予測がないと、いつにどれだけ作って、どれだけ在庫を持っておくべきかというような最適化問題を解くことができないからです。 もう少しデメリットについて具体的に言うと、需要予測ができないと在庫の適正化ができず、担当者の感覚に頼って余裕を持って作りすぎて、結果としてCCC（Cash Conversion Cycle：原料仕入れに現金を投入してから最終的に現金化されるまでの日数）が悪化するなどの弊害が生じます。 また、需要予測をしておくことによって、たとえ予測が外れたとしてもハズレた外れた場合のオペレーション計画を事前に立てることもできます。 なので、精度の良し悪しはありますが、需要予測をしっかりしていくことが非常に重要といえます。 今回使用するサンプルデータセット概要 概要 サンプルのデータセットとして、Kaggle　のStore Item Demand Forecasting Challengeのデータを使って販売実績を予測（需要予測）していきます。 このデータでは10店舗ある小売店のそれぞれで50種類の商品を売っており、各商品について2013/1/1〜2027/12/31までの日付ごとの販売実績(913000行)がトランザクションデータとしてあります。以下にサンプルを示します。 Salesデータの確認 簡単に中身を確認してみます。 今回はある1つの店舗、ある１つの商品に絞って２０１３年から２０１７年の販売実績の推移を確認してみます。 データ確認その1 #データセットの中身を確認。すべてを見るのは大変なので、Store==1、item==1のデータに絞って見てみる。 small_df = df.query('store == "1" & item == "1"') #matplotlibで可視化する import matplotlib.pyplot as plt plt.ylabel("Sales") plt.plot(small_df["date"], small_df["sales"], label="daily_sales") plt.legend(loc = 'upper left') plt.show() 可視化結果は以下です。 １つの年の中で周期性がありそうです。 更に区間を狭めて２０１３年の１年間の販売実績を見てみようと思います。 データ確認その２ #２０１３年に絞ってデータを確認 small_df_1 = small_df.query('year == "2013"') #matplotlibで可視化する plt.plot(small_df_1["date"],small_df_1["sales"], color ="red",label = '2013') plt.ylabel("daily_Sales") plt.xlabel("date") plt.legend(loc = 'upper right') plt.show() 可視化結果は以下です。 どうやらこの商品は季節性があり、６〜８月によく売れているようです。 今回は５０製品のうち１製品だけに着目しましたが、他の製品についても同様に確認することができます。 需要予測の各種手法 以下では代表的な手法を用いてサンプルデータセットについて需要予測をしていきます。 移動平均法 移動平均法は過去の売上から算出される移動平均をもとに需要予測をする単純な方法です。 この手法ではスパイクになっている急激な需要変動を予測することはできませんが、季節ごとの大まかな需要変動のトレンドを捉えることができます。 一般的には移動平均を使った需要予測は昨年の販売実績データの移動平均値を使って今年の需要を予測します。 まずは２０１３年の販売実績について移動平均を求めてみます。 移動平均と可視化 #Store==1、item==1、year==2013のデータに絞る。 small_df = df.query('store == "1" & item == "1" & year =="2013') #今回は前３日間の移動平均をとったデータをデータセットに追加して生データと比較する。 small_df["rolling_mean_3"] = small_df["sales"].rolling(window = 3).mean() #可視化 plt.plot(small_df["date"],small_df["sales"], color ="red",label = 'row_data') plt.plot(small_df["date"],small_df["rolling_mean_3"], color ="blue",label = 'rolling_mean') plt.ylabel("daily_Sales") plt.xlabel("date") plt.legend(loc = 'upper right') plt.show() 求めた2013年のStore1, item1の移動平均と販売実績（Sales）を比較すると以下のようになります。 ギザギザになっている販売実績の変動が均されるような形に移動平均がなっていることがわかります（青線）。 では次に２０１３年の移動平均結果を使って、２０１４の販売実績と比較して見ます。 ２０１４年の実績と２０１３年の移動平均値の比較 #２０１３年と２０１４年のstore=1, item=1に絞る small_df = df.query('store == "1" & item == "1"& year =="2013"') small_df2 = df.query('store == "1" & item == "1"& year =="2014"') #日付列の追加 days = [i for i in range(1,366)] small_df["days"]=days small_df2["days"]=days #移動平均を計算・列に追加 small_df["rolling_mean_3"] = small_df["sales"].rolling(window = 3).mean() small_df2["rolling_mean_3"] = small_df2["sales"].rolling(window = 3).mean() #可視化 plt.plot(small_df2["days"],small_df2["sales"], color ="red",label = '2014_raw') plt.plot(small_df["days"],small_df["rolling_mean_3"], color ="blue",label = '2013_rolling') plt.ylabel("daily_Sales") plt.xlabel("date") plt.legend(loc = 'upper right') plt.show() #MAEで精度の確認 small_df = small_df.dropna()#window=3の移動平均なので最初２つの行はNaNになっている from sklearn.metrics import mean_absolute_error print(mean_absolute_error(small_df2['sales'][2:], small_df['rolling_mean_3'])) #4.992653810835629 以下に可視化の結果を示します。 2013年の移動平均結果と２０１４年の販売実績を比較すると多少の需要の変化はありますが、概ね季節変動を捉えてられており、２０１３年実績の移動平均をもとに２０１４年の需要を想定することができることがわかります。 加重移動平均法 加重移動平均法は最新の需要変動の影響に重み付けする手法で、移動平均を取る際に荷重係数を掛け合わせることによって求められます。 直前の変動の影響が大きい（例えば連休やイベント）ような製品の需要予測ではこちらのほうが移動平均よりも精度がよくなることがあります。 加重移動平均（WMA：Weighted Moving Average）を求めるのにはTA-Libというチャート分析用ライブラリを使います。 先程定義したsmall_dfとsmall_df2を使ってやっていきます。 加重移動平均（WMA)と可視化 #ライブラリのインポート import talib #データの前準備 date = small_df["date"] small_df.drop('date', axis=1, inplace=True) small_df = small_df.astype("float64") #加重移動平均の算出 small_df["wma_3"]=talib.WMA(small_df["sales"].values, timeperiod=3) #可視化 plt.plot(small_df2["days"],small_df2["sales"], color ="red",label = '2014_raw') plt.plot(small_df["days"],small_df["wma_3"], color ="blue",label = '2013_wma_3') plt.ylabel("daily_Sales") plt.xlabel("date") plt.legend(loc = 'upper right') plt.show() #MAEで精度を確認 from sklearn.metrics import mean_absolute_error print('MAE(wma):') print(mean_absolute_error(small_df2['sales'][2:], small_df['wma_3'])) #4.952249770431588 以下に可視化の結果を示します。 先程の移動平均よりもグラフの振動が大きくなっている一方で、MAEも若干ですが加重移動平均のほうが低くなっており、２０１４年の販売実績と比較するとよりトレンドが一致しているように見えます。 移動平均と加重移動平均は何をしたいかの目的に応じて使い分けると良さそうです。 指数平滑法（EWM：Exponentially weighted method） 重み付けをするという意味では加重移動平均法と近いですが、こちらは指数関数的な重み付けをする手法です。 加重移動平均法のときと同様のデータで比較してみます。 指数平滑と可視化 #指数平滑移動平均の算出 small_df["ewm"]=small_df["sales"].ewm(alpha=0.5).mean() #可視化 plt.plot(small_df2["days"],small_df2["sales"], color ="red",label = '2014_raw') plt.plot(small_df["days"],small_df["ewm"], color ="blue",label = '2013_ewm') plt.ylabel("daily_Sales") plt.xlabel("date") plt.legend(loc = 'upper right') plt.show() #MAEで精度を確認 from sklearn.metrics import mean_absolute_error print('MAE(ewm):') print(mean_absolute_error(small_df2['sales'], small_df['ewm'])) #4.811589224904538 以下に可視化の結果を示します。 ぱっと見る限り加重移動平均で求めたものとあまり変わりませんが、大きな違いとして加重移動平均では移動平均を取る区間を自分で設定するのに対して、指数平滑法では平滑パラメーターalphaを調節して目的に合う予測値を出して行きます。 今回のデータでは指数平滑法のほうが移動平均や加重移動平均よりも誤差が小さくなっていますが、扱うデータによって一長一短がありそうですので、目的にあった手法を選ぶことがいいかなと思います。 時系列分析法 （機械学習：NeuralProphet） これまではクラシックな需要予測の手法を使って直感的にわかりやすい需要予測をしてきました。 一方でこういった手法は人間の直感・経験に頼った処理（例えば移動平均の区間をどうやって取るか）が必要で、使い手によって精度はまちまちになりそうです。 そこでここではFacebookから公開されているNeuralProphetというライブラリを使った機械学習による需要予測を実装してみようと思います。 NeuralProphetは、ProphetとAR-Netという自己回帰型のニューラルネットワークモデルを（Deep Learning）元にした時系列データ予測モデルです。 Deep Learningと聞くと、ブラックボックスで解釈性に乏しいという印象があると思いますが、NeuralProphetは学習後にモデルの捉えたトレンドや周期性の可視化ができ、何をもとに予測されているのかを確認できるため、比較的結果の解釈がしやすいというのが特徴のライブラリになっています。 学習データとしてstore1,item1について、2013~2016年のSalesデータを用い、２０１７年の需要予測（Salesの予測）を行います。 NeuralProphetで需要予測 #NeuralProphetをインポート from neuralprophet import NeuralProphet #データセットの準備 df1 = pd.read_csv("train.csv") df1 = df1.query('store == "1" & item == "1"') df1 = df1.drop(["store", "item"], axis=1) df1 = df1.rename(columns = {"date": "ds", "sales": "y"}) #データセットの分割 train = df1[:-365]#2013-2016年の販売実績 test = df1[-365:]#2017年の販売実績 #機械学習モデルのインストタンス作成 m = NeuralProphet(seasonality_mode='multiplicative') #validationの設定(3:1)&モデルのトレーニング metrics = m.fit(train,freq="D") #作ったモデルで２０１７年のSalesを予測(trainデータの先３６５日分を予測) future = m.make_future_dataframe(train, periods=365) forecast = m.predict(future) #予測結果をtestのデータフレームに追加 test["pred"] = forecast["yhat1"].to_list() #MAEで精度を確認 from sklearn.metrics import mean_absolute_error print('MAE(NeuralProphet):') print(mean_absolute_error(test['y'], test['pred']))#3.9889207134508107 #可視化 test.plot(title='Forecast evaluation',ylim=[0,50]) 以下に可視化の結果を示します。 得られた予測結果は先程までのクラシックな手法に比べてかなり異なる挙動に見えますが、MAEはかなり改善されました。 過去のデータの単純な需要予測手法に比べてDeep Learningを使ったNeuralProphetが優れていることがわかります。 この予測根拠になったトレンド分析結果をplot_components()とplot_parameters()を使って確認してみます。 NeuralProphetのトレンド分析結果を確認 fig_comp = m.plot_components(forecast) fig_param = m.plot_parameters() こちらの分析結果は予測対象の２０１７年のトレンドと季節変動をモデルの中で分解した結果です。 なだらかな上昇傾向と年間の季節変動、１週間の間の季節変動が自動的に算出されていることがわかります。 以下の２つ目のグラフは学習データの分析データです。 数年間のトレンドの上昇傾向、各年の時期ごとでの変動率等の情報が出力されています。 これらのデータを見ることで、おおよその過去の傾向のsummaryがわかり、かつ予測結果が何から導かれているのかが理解できます。 おわりに 需要予測について様々な手法を見てきました。 SCMにおける需要予測の重要性については最初にお伝えしたとおりですが、一方で需要予測は万能ではありません。 COVID-19のパンデミックによって既存の需要予測はほとんど当てにならないものになったという話を聞きました。 それは上記で紹介した技術すべてに共通しているのですが、パンデミックや戦争のような人生で1度きりみたいな大きく不連続な変化を過去のデータから予測することは難しいからです。 ですので、需要予測するにあたっては、その要素が何で、どんな前提に基づいた予測なのかを知っておくことが重要です。 そうすれば、需要予測が外れ始めたときに何が原因なのか、どうすればいいのか(場合によっては需要予測のモデルを使わないという選択をする)について考える指針になると思いますので、本記事がその参考になれば幸いです。 参考文献 需要予測で在庫管理を効率化！計算式や精度を上げる方法を紹介 Python言語によるサプライ・チェイン・アナリティクス 【需要予測の手法】製造業における在庫管理の改善・在庫削減に活かす Machine Learning for Retail Demand Forecasting 【FX,Python】PythonとTA-Libで加重移動平均(WMA)を引く Mac OSでTA-Libをpipで入れられなかった時の対処法 Pythonでサクッと作れる時系列の予測モデルNeuralProphet（≒FacebookのProphet × Deep Learning） NeuralProphetで時系列データ予測 NeuralProphet公式

こんにちは！データサイエンティストのじゃむです。 本日は、時系列データで使用する移動平均と移動分散を逐次的に求めるアルゴリズムの紹介です。 問題設定 長さTの時系列データX=(X_t)_{t=1,・・・,T}に対して、k項移動平均と移動分散を逐次的に求めることを考えます。 たとえば、X=[1,2,3,4,5]、k=2として、 移動平均 = [1.5, 2.5, 3.5, 4.5] 移動分散 = [0.25, 0.25, 0.25, 0.25] をどうやったら効率よく算出できるかを考えます。 単純な方法 たとえば先ほどの例でX=[1,2,3,4,5]の2項移動平均を愚直に求めるとすると (1+2)/2 = 1.5, (2+3)/2 =2.5 ・・・ 同じように、2項移動分散を求めると、 ((1-1.5)^2 + (2-1.5)^2)/2 = 0.25, ((2-2.5)^2 + (3-2.5)^2)/2 = 0.25, ・・・ となります。 一般的には、j番目のk項移動平均 $m_j$ と移動分散 $v_j$ は次のようにかけます。(j=0,1,2,・・・, T-1) $$ m_j = \frac{1}{k}\sum_{i=j}^{j+k-1} X_i $$ $$ v_j = \frac{1}{k}\sum_{i=j}^{j+k-1} (X_i-m_i)^2 $$ このままでは、計算量が高々O(T*k)になるので、今回のアルゴリズムでは、毎回Σ計算をせずに計算量がO(T)となるように工夫します。 なにをするか ０番目のk項移動平均の計算は、 $$ m_0 = \frac{1}{k}\sum_{i=0}^{k-1} X_i $$ ですが、1番目の移動平均を計算するときに、0番目の移動平均を使うことを考えます。1番目の移動平均を計算するためには、$X_1$から$X_k$までの合計値が必要です。これは$X_0$から$X_{k-1}$までの合計値から、$X_k$を足して、$X_0$を引いたものと考えれば、 $X_1$ から $X_k$ までの合計値は $m_0k - X_0 + X_{k}$ と算出することができます。よってΣの計算をせずに各移動平均は、O(1)で求められ、$i$番目のk項移動平均は次のように逐次的に求められます。 $$m_i = m_{i-1} - \frac{X_{i-1} - X_{i+k-1}}{k}$$ 次に、移動分散ですがこれは少々厄介で、次の分散公式が必要です。 $$\frac{1}{n}\sum_{i=1}^{n} (X[i]- m)^2 = \frac{1}{n}\sum_{i=1}^{n} X^2[i] - m^2 $$ ここで、$m$ は$X$ の平均です。これを使って、$v_j$ を書き換えると、 $$ v_j = \frac{1}{k}\sum_{i=j}^{j+k-1} (X[i]-m_i)^2 = \frac{1}{k}\sum_{i=j}^{j+k-1} X^2[i] -m^2_i$$ になります。移動平均と同じように$v_0$を求めると、上記の式から $$ v_0 = \frac{1}{k}\sum_{i=0}^{k-1} X^2[i] -m^2_0$$ $v_1$を求めるには、$X_1$から$X_k$までの二乗和が必要ですが、移動平均と同じようにすれば、$X_0$から$X_{k-1}$までの二乗和から、$X_k$の二乗を足して、$X_0$の二乗を引けばよいです。よって、$X_1$から$X_k$までの二乗和は$k(v_0+m^2_0) + X^2_k - X^2_0$ になります。よって、分散公式から、$v_1$は、 $$ v_1 = (v_0+m^2_0) + \frac{X^2_k - X^2_0}{k} - m^2_{1} $$ と算出できます。よってΣの計算をせずに各移動分散は、O(1)で求められ、$i$番目のk項移動分散は次のように逐次的に求められます。 $$ v_i = (v_{i-1}+m^2_{i-1}) + \frac{X^2_{i+k-1} - X^2_{i-1}}{k} - m^2_{i} $$ ソースコード pythonで実装すると、次のようになります。 import numpy as np def calc_moving_kth_average_and_variance(x_array, k): m_init = np.mean(x_array[:k]); mean_list = [m_init] #平均の初期値 v_init = np.var(x_array[:k]); var_list = [v_init] #分散の初期値 for i in range(len(x_array)-k): # 更新式 mean_list.append(mean_list[-1] + (-x_array[i] +x_array[i+k])/k) var_list.append(var_list[-1]+mean_list[-2]**2 +(x_array[i+k]**2 - x_array[i]**2)/k -mean_list[-1]**2) return mean_list, var_list

はじめに 競艇のルール，予想の仕方など私自身何もその辺の知識がないことをご了承ください． それらすべてをデータ分析，機械学習に任せることによって知識のない人間でも利益をプラスにできるのではないか，そういう魂胆でこの記事を書いています． 環境 MacOS BigSur 11.4 Jupyter Notebook(Colaboratoryでも可) データについて 競艇公式ホームページ BOAT RACEからスクレイピングによってデータを取得しました． 会場：住之江会場 データ数：43560件 期間：2019/01 ~ 2020/04 今回の目的 競艇はスタート位置でほぼ決まる．みたいな話を聞いたことがあります． データ分析によってこれが本当なのかを調べるのが今回の目標です． 本編 ライブラリのimport import numpy as np import pandas as pd import matplotlib.pyplot as plt import seaborn as sns from itertools import product 　前処理 データはsuminoe_40000_program.csvとsuminoe_40000_answer.csvに格納されているため，PandasによってDataFrame形式でデータを読み込みます． 説明変数 names = ["Place", "Number", "Name", "Age", "Live", "Weight", "Rank", "A_1st", "A_2nd", "B_1st", "B_2nd", "Moter_No", "Moter_2nd", "Bote_No", "Bote_2nd"] df_program = pd.read_csv("../../Boat-Race-Prediction/suminoe_40000_program.csv", names=names) print(f"shape: {df_program.shape}") df_program.head() 目的変数 df_answer = pd.read_csv("../../Boat-Race-Prediction/suminoe_40000_answer.csv", names=["idx", "target"]) print(f"shape: {df_answer.shape}") df_answer OneHotEncoding 説明変数を見ると，Name, Live, Rankはカテゴリ変数となっていることがわかります． 名前Nameと地域LiveについてはOneHotEncodingしてしまうと　，列数がとんでもないことになってしまうため今回はデータから弾くことにします． よってRankをOneHorEncodingによって，A1, A2, B1, B2の列を新たに作成します． (目的変数のidxも必要ないため弾いています．) df_data = df_program.merge(df_answer, how="inner", right_index=True, left_index=True) df_data = df_data.drop(["idx", "Name", "Live"], axis=1) df_dummies = pd.get_dummies(df_data, columns=["Rank"]) df_dummies.head() 目的変数の前処理 競艇には転覆やスタートミスによって失格判定になる場合があります． 生データではtargetの部分がK, Fのようなカテゴリ変数となっています． できるだけ数値のみのデータとして扱いため，10, 11, 12の値に変換してあります． 今回の分析でこれらのデータは使用しないため，np.nanによって欠損値扱いとします． df_dummies.target[df_dummies.target == 10] = np.nan df_dummies.target[df_dummies.target == 11] = np.nan df_dummies.target[df_dummies.target == 12] = np.nan 各列とtargetの相関を見る targetと各変数の関係を相関係数を計算することによって調べます． わかりやすさのために，Seaboanのヒートマップを使用しています． # 相関係数を表示 plt.rcParams['figure.figsize'] = (20.0, 10.0) plt.rcParams['font.family'] = 'serif' sns.heatmap(df_data.corr(), annot=True, square=True, cmap='coolwarm'); このヒートマップから，targetとPlaceの間には，0.38の正の相関があることがわかります． 0.38という値は強い相関とは言えませんが，やや相関があると言えます． なので，targetとPlaceをヒストグラムで示してみましょう． targetとPlaceの関係をヒストグラムで可視化する 各スタート位置に対する順位をヒストグラムによって可視化してみます． これによって位置がどのくらい順位に影響するのかがわかるはずです． fig, ax = plt.subplots(2, 3, figsize=(15, 10)) for i, tu in zip(range(1, 7), product(range(2), range(3))): df_place = df_dummies[df_dummies.Place == i] ax[tu[0], tu[1]].hist(df_place.target, bins='auto') ax[tu[0], tu[1]].title.set_text(f'Place {i}') ax[tu[0], tu[1]].set_xlabel('Ranking') ax[tu[0], tu[1]].set_ylabel('Count') plt.tight_layout() スタート位置が1である場合に1位になる確率はとてつもなく高そうです． スタート位置が6に近づくにつれて1位になる確率は下がっていっていることがわかります． では，位置に対する順位の確率を算出してみましょう． for i in range(1, 7): df_place = df_dummies[df_dummies.Place == i] print(f"Place: {i}") for j in range(1, 7): df_target = df_place[df_place.target == j] probability = df_target.shape[0] / df_place.shape[0] * 100 print("\t{}: {:.1f}%".format(j, probability)) Place: 1 1: 55.9% 2: 16.9% 3: 9.1% 4: 6.8% 5: 5.5% 6: 4.6% Place: 2 1: 15.7% 2: 25.4% 3: 18.3% 4: 15.1% 5: 13.0% 6: 11.3% Place: 3 1: 11.3% 2: 21.5% 3: 20.3% 4: 18.2% 5: 15.1% 6: 12.2% Place: 4 1: 8.7% 2: 16.1% 3: 19.5% 4: 19.0% 5: 17.8% 6: 17.3% Place: 5 1: 5.4% 2: 11.7% 3: 17.0% 4: 20.0% 5: 23.0% 6: 21.4% Place: 6 1: 3.0% 2: 8.4% 3: 15.7% 4: 20.6% 5: 24.4% 6: 26.4% スタート位置が1であったときに1位になる確率は，55.9%であることがわかりました！ よって，1レース6人であるため，データでは1レース6行で表現されます．なので，今回のデータは43560件であるため， 43560 / 6 = 7260レース分のデータが存在しています． この7260レース中のうち，55.9%が1位になるため，スタート位置が1である選手は7360中4058回１位になっていることが今回のデータからわかりました． まとめ 絶対というわけではないですが，スタート位置だけでも順位にそれなりに影響するようです． しかしスタート位置が2, 3, 4の場合は順位にあまり差がないため，3連単などを当てるためには他の点で分析していく必要があると思います． ここまでのデータ分析で使えるのは単勝予想くらいだと思います．

はじめに Python Tweepyを使ってツイートの送信・取得・検索やユーザーのフォロワーの確認、DMの確認・送信をしてみた。 TweepyのinstallとTwitter Developerの登録 　tweepyを使うにはTwitter APIが必要なので、下記の公式tweepyドキュメンテーションを見て、twitter developersへ登録する。もちろんクライアント側はpython moduleも必要なので、tweepyをinstallする。tweepyの使い方も以下のリンク資料に書いてある。 pip install tweepy 認証キーの設定とOAuth認証、APIのインスタンス生成 　IDやパスワードが不要でアプリケーション間の連動ができるOAuth認証を採用しており、各ユーザーがアクセスキーを生成してそのアクセスキーで認証を行う。Webシステムとの対話は、URL形式でHTTP Request/Responseをするtwitter REST APIでアクセスすることができ、得られたレスポンスはModel class instanceとして、INSTANCE.CLASSで各classにアクセスすることができる。 consumer_key = "API key" consumer_secret = "API secret key" access_token = "Access token" access_token_secret = "Access token secret" auth = tweepy.OAuthHandler(consumer_key,consumer_secret) auth.set_access_token(access_token,access_token_secret) api = tweepy.API(auth) 1. ツイートする user = 'my_userid' message = "Hello!" api.update_status(message) 2. 特定のユーザーのtimeline取得 user = 'qiita_python' for i, status in enumerate(tweepy.Cursor(api.user_timeline,id=user,tweet_mode="extended").items(10)): print(i," : ",status.full_text) 3. 特定のユーザーをフォローしている人フォローされている人を取得 user = "qiita_python" print("Name : @Name") for status in tweepy.Cursor(api.friends,id=user).items(10): print("%s : %s "%(status.name,status.screen_name) ) print("-"*100) for status in tweepy.Cursor(api.followers,id=user).items(10): print("%s : %s "%(status.name,status.screen_name) ) 4. 特定のユーザーのプロフィールを取得 user_id = "qiita_python" user = api.get_user(user_id) print(user.description) print(user.profile_image_url_https) 5. 特定のユーザーのメディアファイルをダウンロードする key_account = 'qiita_python' search_results = tweepy.Cursor(api.user_timeline, screen_name=key_account, include_rts=False).items(10) for result in search_results: if hasattr(result, 'extended_entities'): ex_media = result.extended_entities['media'] for image in ex_media: image_url = image['media_url'] image_name = image_url.split("/")[len(image_url.split("/"))-1] urllib.request.urlretrieve(image_url + ':orig', 'img/' + image_name) else: print('No media') 6. 特定ユーザーへDMの送信 ID = 'qiita_python' text = 'Hello!' event = { "event": { "type": "message_create", "message_create": { "target": { "recipient_id": ID }, "message_data": { "text": text } } } } api.send_direct_message_new(event) 7. 受け取ったDMの確認 direct_messages = api.list_direct_messages() for direct_message in direct_messages: print(direct_message) 8. 特定箇所から半径何km圏内のツイートを検索して、ユーザーと内容を取得 for status in tweepy.Cursor(api.search,q="東京駅",geocode="35.681698005508686,39.7671958524044,10km").items(10): if "location" in status._json["user"].keys(): print("%s : %s : %s" %(status.created_at,status.user.location,status.text)) else: print("%s : nan : %s" %(status.created_at,status.text)) まとめ Python Tweepyを使ってみた。 できることは以下の公式ページに書いてるので、それを参考にすればしたいことができると思います。 Twitter API v1.1 Reference

起きた問題 VSCode上でPylanceという拡張機能を使っていて sample.py import numpy as np import matplotlib.pyplot as plt などと書くとエラーを示す波線が出るのですが実際実行してみると問題なく実行できるという状態になりました。しかし、この「問題ない」というのはあくまでコードが正しく書けていればの話でもしnumpyのメソッドにtypoがあったときなどにVSCodeが「間違っているよ！そんなメソッドは存在しないよ！」と言ってくれなくなります。それにimportする度にエラーの波線が出るのは鬱陶しいので解決法を調べました。 環境 MacOS BigSur, python 3.6.6 64-bit やったこと 初めにsettings.jsonというファイルをいじりました。このファイルはshift+command+pを押して、settings.jsonと打ち込めば設定画面に移ることができます。ここで設定画面には「ユーザー」と「ワークスペース」というものがあると思うのですが個人的なイメージでは「ユーザー」はVSCode全体の設定、「ワークスペース」は今開いているフォルダに対する設定のようなものだと認識しています。 まず「ユーザー」を選択してからVSCodeの右上の「ファイルに矢印がついているアイコン」をクリックします。すると「ユーザー」のsettings.jsonが開きます。そこにどの行でもいいので(正確には他の行のネストに入らない行)にターミナルで terminal which python として出力されたpythonのパスをsettings.jsonに加えます。私の場合は下のようなコードを加えました。 User/settings.json "python.pythonPath": "/Users/ユーザー名/.pyenv/shims/python" これによってVSCodeで動作するpythonとlocal環境で動作するpythonを合わせることができるはずです。私の場合はこれでimport errorの波線は消えました。 その後「ワークスペース」を選択して同じくVSCodeの右上の「ファイルに矢印がついているアイコン」をクリックすると「ワークスペース」のsettengs.jsonが開きます。それか今のフォルダに新しく .vscode というフォルダを作ってそのフォルダの中に settings.json を作っても構いません。どちらでも同じです。このファイルには .vscode/settings.json { "python.languageServer": "Pylance" } と書きます。これで今のワークスペースでPylanceを使うことができるようになります。 その後にVSCodeの左下のPythonのバージョンを選択できるところがあるのでそれをはじめのsettings.jsonに書いたPythonのバージョンと同じものを指定します。 最後に再びcommand+shift+pを押して、settings.jsonと打って「ユーザー」を選択します。すると左側のタブに拡張機能という欄があるのでスクロールしていくとPylanceがあり、そこで Diagnostic Severity Overrides という欄を変更します。ここでは不具合、バグの種類によって波線の色を変えることができます。ここの裁量は自分で決めてしまっていいと思います。コードを書いていってみてこの種類の不具合はwarningでいいかなとかこの不具合はerrorだなとか決めていけばいいと思います。 終わりに 以上が自分が調べてやってみたことです。ここ間違ってるよとかここをこう変えるとより良いよなどございましたらぜひコメントお願いします！

はじめに 　調査観察データの統計科学に記載の傾向スコアを使った因果効果の推定手法とpythonコードをまとめました。本記事では最初に傾向スコアの定義をして、次に共変量調整の方法を記載します。共変量調整の方法としては、傾向スコアマッチング、層別解析、IPW、カーネルマッチングを記載しました。最後にこれらの手法をpythonで実装し、シミュレーションを実施した結果を記載します。 \def\middlemid{\;\middle|\;} \def\set#1{\left\{{#1}\right\}} \def\setin#1#2{\left\{{#1} \middlemid {#2}\right\}} \def\sub#1#2{{#1}_{#2}} 傾向スコア 　傾向スコアの定義は下記のようになります。第$i$被験者の介入変数を Z_i = \left\{ \begin{array}{ll} 1 & (介入があった場合) \\ 0 & (介入がなかった場合) \end{array} \right. とします。また、第$i$被験者の共変量を$\boldsymbol{x}_i$とします。このとき、第$i$被験者の傾向スコアは $$e_i = P(Z_i=1|\boldsymbol{x}_i)$$ と定義されます。 　傾向スコアを使って共変量を調整することで、セレクションバイアスを排除した介入の効果を推定することが今回の目的です。傾向スコアを使ってセレクションバイアスを排除する必要性については、効果検証入門にめちゃくちゃわかりやすく書かれています。 共変量調整の方法(理論) 傾向スコアマッチング 　最も一般的(?)な方法。介入群と対照群で傾向スコアが同じ被験者をペアにしてその結果変数の差を取り、その差の平均を因果効果の推定値とする手法です。数式で表すと下記のようになります。 \tau = E_e[E[Y_1|e, Z=1]-E[Y_0|e, Z=0]].　\tag{1} ただし、$Y_1, Y_0$はそれぞれ介入があった場合の結果変数と、介入がなかった場合の結果変数です。また、$\tau$は推定対象の因果効果で$\tau = E[Y_1-Y_0]$です。1この数式をつぶさに見てみると、最初に$E[Y_1|e, Z=1]-E[Y_0|e, Z=0]$を計算し、それを傾向スコア$e$について平均を取っています。感覚的には$E[Y_1|e, Z=1]-E[Y_0|e, Z=0]$の計算が「介入群と対照群で傾向スコアが同じ被験者をペアにしてその結果変数の差を取り」に対応し、$e$による平均の計算が、「その差の平均」の計算に対応しています。 　傾向スコアマッチングは直感的で非常にわかりやすい一方で、傾向スコアが完全に同じ介入群と対照群の被験者をペアにできるのかという問題をはらんでいます。実際、傾向スコアは0から1の間の実数なので、傾向スコアが完全に同じペアは存在しない場合が多いです。そこで、傾向スコアが最も近い被験者同士をマッチングさせる最近傍マッチングや、介入群の各被験者の傾向スコアに最も近い$M$人の対照群の被験者をマッチングさせる1対多のマッチングなどが提案されています。後で記載するpythonのコードには最も基本的な最近傍マッチングを実装しています。 層別解析 　傾向スコアマッチングでは、完全に同じ傾向スコアでマッチングさせることは不可能であることが多いと書きました。そこで、マッチングさせるのではなく傾向スコアの値によって層を複数作り、各層で因果効果を推定し、それらの平均をとることで因果効果を推定するという方法が提案されました。この手法を層別解析と呼びます。 　具体的に数式で表すと因果効果の推定量は、下記となります。 \hat{\tau} = \sum_{k=1}^K \frac{N_k}{N} (\bar{y_1}_k-\bar{y_0}_k) \tag{2} ただし、被験者全体のサンプルサイズを$N$、第$k$層のサンプルサイズを$N_k$、第$k$層の介入群の結果変数の平均を$\bar{y_1}_{k}$、対照群の結果変数の平均を$\bar{y_0}_k$としています。 　層の数については分析者が決めるパラメータであり、調査観察データの統計科学によれば、5層が使われることが多いそうです。(理論的な根拠はないそうですが...。 IPW推定量 　傾向スコアマッチングと層別解析の基本的な考え方は、介入群と対照群の被験者をペアにすることで因果効果を推定します。しかし、これらの方法はマッチングの仕方や層数の決め方が分析者任せであり、恣意的であるとの批判がありました。そこで、傾向スコアを使って結果変数を重みづけることで因果効果の推定量を導く方法としてIPW推定量が提案されています。 　IPW推定量とは、 \hat{\mu_1} = \sum_{i}^N \frac{Z_i Y_i}{e_i}/\sum_{i=1}^N\frac{Z_i}{e_i}, \ \ \hat{\mu_0} = \sum_{i=1}^N \frac{(1-Z_i) Y_i}{1-e_i}/\sum_{i=1}^N \frac{1-Z_i}{1-e_i} \tag{3} です。ただし、$\hat{\mu_1},\hat{\mu_0}$は $\mu_1 = E[Y_1], \mu_0 = E[Y_0]$の推定量です。2 　IPW推定量の良さの一つは、漸近正規性を持つことです。すなわち、 $$\sqrt{N}(\hat{\mu_1}-\hat{\mu_0}-\mu_1-\mu_0) \stackrel{d}{\longrightarrow } N(0, V(\theta_0))$$ が成り立ちます。3　したがって、被験者数が十分に大きいときは、IPW推定量を使うことを検討すると良いと考えています。(実装も簡単！) カーネルマッチング 　傾向スコアマッチングの課題は、傾向スコアが同じ被験者が介入群と対照群に存在しないことでした。それに対してカーネルマッチングの基本的な考え方は傾向スコアが同じ被験者が介入群と対照群に存在しないなら作ってしまえです。 　ということで、傾向スコアから結果変数を予測するモデルを作成し、介入群の介入を受けていない場合の結果変数を傾向スコアから予測し、観測されている介入があったときの結果変数から、予測された結果変数を引くことで因果効果を推定する手法が提案されています。文章で書くと意味不明なので、手順を追って説明します。 対照群の傾向スコア$e_0$から結果変数$Y_0$を予測するモデルを作成する。これを$Y_0 = f(e_0)$とする。 介入群の各被験者の傾向スコア$e_{1i}$から、介入を受けていない場合の結果変数を1のモデルで予測する。これを$\hat{Y}_{0i} = f(\sub{e}{1i})$とする。 観測されている介入群の結果変数から、2で予測した介入群の介入を受けていない場合の結果変数を引く。これを$\sub{Y}{1i}-\sub{\hat{Y}}{0i}$とする。 3で求めた差分の平均を取り、因果効果を推定する。すなわち、$\hat{\tau} = \frac{1}{\sub{N}{1}}\sum_{i=1}^{N_1}(\sub{Y}{1i}-\sub{\hat{Y}}{0i})$とする。($N_1$は介入群のサンプルサイズ) この方法で大事なのは、モデル$f$の精度です。傾向スコアから結果変数を予測するモデルは無数に存在しますが、ここでは局所回帰モデルすなわち、カーネルで重みづけられた回帰モデルを使います。そのため、カーネルマッチングという名前がついています。局所回帰モデルについては、局所線形回帰を完全に理解したいが超絶分かりやすかったです。また、調査観察データの統計科学の付録でも局所回帰の雰囲気は理解できます。 手法の比較 　ここまで、傾向スコアマッチング・層別解析・IPW推定量・カーネルマッチングについて述べてきました。手法がいくつか出てくると、それらを比較してより良い手法を選択したくなりますよね。効果検証入門によれば、平均のつり合いを見ることである程度手法の精度を評価することができます。 　具体的には、各手法ごとに下図を作成して比較します。 上図の縦軸は共変量、オレンジ色、青色の棒グラフはそれぞれ調整前と調整後の介入群と対照群の共変量の差の絶対値を表しています。これらは各手法によって算出の仕方が異なりますが、基本的には介入群と対照群の共変量の平均を出すという考え方で問題ないです。介入群の第$i$被験者の共変量を$X_{m1i}$、対照群の第$j$被験者の共変量を$X_{m0j}$とします。($m$は上図では1,2,3,4のいずれか) 調整前の共変量の差の絶対値 介入群と対照群の各共変量の平均を算出する。つまり、$\sub{\bar{X}}{m1} = \frac{1}{N_1}\sum_{i=1}^{N_1}X_{m1i}$と$\sub{\bar{X}}{m0} = \frac{1}{N_0}\sum_{j=1}^{N_0}X_{m0j}$を計算する。 共変量の平均の差の絶対値を計算する。つまり、$|\sub{\bar{X}}{m1}-\sub{\bar{X}}{m0}|$を計算する。これが、オレンジ色の棒グラフである。 調整後の共変量の差の絶対値 傾向スコアマッチングの場合 介入群の第$i$被験者と対照群の第$j$被験者がマッチングされたとする。このとき、それぞれの共変量の差分を計算する。つまり、$X_{m1i}-X_{m0j}$を計算する。 共変量の差分の平均をとり、絶対値を計算する。つまり、$|\frac{1}{N_1} \sum(X_{m1i}-X_{m0j})|$を計算する。これが、傾向スコアマッチングの場合の青色の棒グラフである。 層別解析の場合 (2)式の$\bar{y_1}_k, \bar{y_0}_k$をそれぞれ$\sub{\bar{X}}{m1}, \sub{\bar{X}}{m0}$に入れ替える。 入れ替えた後の式に絶対値をつけた値を計算する。これが層別解析の場合の青色の棒グラフである。 IPW推定量の場合 (3)式の$Y_i$をそれぞれ$X_{mi}$に置き換える。 置き換えた$\hat{\mu}_0, \hat{\mu}_1$の差分の絶対値を計算する。これがIPW推定量の場合の青色の棒グラフである。 カーネルマッチングについては、上記のような評価をすることができません。したがって、傾向スコアマッチングとカーネルマッチング間で比較をすることは難しいです。しかし、カーネルマッチングで用いられるモデル間での評価をすることはできます。カーネルマッチングにおいては、重みづけに使うカーネル、カーネルのバンド幅、回帰の次元数などがハイパーパラメータとなりますが、これらを変えて汎化性能が良いモデルを選択することで、ある程度の精度は担保できると考えられます。 シミュレーション シミュレーションデータ 下記の条件でシミュレーションデータを生成しました。 サンプルサイズ：$N=10^3$ 共変量の次元数：$M = 4$ 共変量の分布：全て$[-1, 0]$の一様分布 真の傾向スコアのモデル：$P(Z=1|X_1, X_2, X_3, X_4) = \frac{1}{1+{\rm exp}(-X_1-X_2-X_3-X_4)}$ 真の結果変数と共変量・介入の関係：$Y = X_1+2X_2+3X_3+4X_4+Z$ この条件に合わせてシミュレーションデータを発生させます。 シミュレーションデータの発生 #import import numpy as np import pandas as pd import matplotlib.pyplot as plt #シミュレーション用のパラメータの設定 ##サンプルサイズ N = 10**3 ##共変量の次元数 K = 4 cov_list_str = ["X1", "X2", "X3", "X4"] ##一様分布の範囲 a, b = [-1, 0] #汎用関数 ##真の傾向スコアのモデル def propensity_score(vec, coef = np.array([-1, -1, -1, -1])): linear = vec*coef return(1/(1+np.exp(np.sum(linear)))) ##真の結果変数 def result(vec, coef = np.array([1, 2, 3, 4, 1])): linear = vec*coef return(np.sum(linear)) #シミュレーションデータの生成 ##共変量作成 df = pd.DataFrame((b-a)*np.random.rand(N,K)+a, columns = cov_list_str) ##真の傾向スコアの計算 true_propensity_score = pd.DataFrame(df.apply(propensity_score, axis=1), columns = ["true_propensity_score"]) ##真の傾向スコアをもとのデータフレームにマージ df_true_ps = pd.merge(df, true_propensity_score, how="left", left_index=True, right_index=True) ##真の傾向スコアに基づき介入変数を乱数で発生させる treatment = pd.DataFrame([np.random.binomial(n=1, p=p) for p in true_propensity_score["true_propensity_score"]], columns = ["treatment_flg"]) ##介入変数をもとのデータフレームにマージ df_treatment = pd.merge(df_true_ps, treatment, how="left", left_index=True, right_index=True) ##結果変数を計算 Y = df_treatment[cov_list_str+["treatment_flg"]].apply(result, axis=1) result = pd.DataFrame(Y, columns = ["result"]) ##結果変数をもとのデータフレームにマージ df_result = pd.merge(df_treatment, result, how="left", left_index=True, right_index=True) このようにして発生させたデータについて、単純に介入群と対照群の結果変数を比較すると下記になりました。(乱数で発生させているので、毎回結果は変わります) 介入群の結果変数の平均 対照群の結果変数の平均 介入群ー対照群 -3.30 -5.09 1.79 真の因果効果は、$Y = X_1+2X_2+3X_3+4X_4+Z$の$Z$の係数であることがわかっているので、今回のセレクションバイアスは$1.79-1=0.79$です。 傾向スコアの推定 　今回は傾向スコアの推定はロジスティック回帰モデルを用いて行います。 傾向スコアの推定 #ロジスティック回帰モデルのインスタンス作成 lr = LogisticRegression() #ロジスティック回帰モデルの推定 lr.fit(df_result[cov_list_str], df_result["treatment_flg"]) #推定された傾向スコア est_propensity_score = lr.predict_proba(df_result[cov_list_str])[:,1] #傾向スコアをもとのデータフレームにマージ df_estps = ( pd.merge( df_result, pd.DataFrame(est_propensity_score, columns = ["est_propensity_score"]), how="left", left_index=True, right_index=True) ) 今回は、真の傾向スコアのモデルがロジスティック回帰モデルであることがわかっているので、モデルの評価は省略しますが、実務上では$AUC>0.7$であることが望ましいそうです。(ロジスティック回帰分析と傾向スコア解析) 今回シミュレーションで作ったデータフレームに含まれるカラムは下記です。 X1 X2 X3 X4 true_propensity_score treatment_flg result est_propensity_score 共変量 共変量 共変量 共変量 真の傾向スコア 介入変数 結果変数 傾向スコアの推定値 傾向スコアマッチング 最近傍マッチングで傾向スコアマッチングを実行するためのclass。クラス内のメソッドや引数の意味は後述します。 傾向スコアマッチングをするclass #傾向スコアマッチング(最近傍マッチング)を実行するクラス class propensityScoreMatching: #__init_で傾向スコアマッチングを行い、indexの組み合わせを取得する def __init__(self, df, cov_list, z, ps, rs): self.df = df self.cov_list = cov_list self.z = z self.ps = ps self.rs = rs #介入の有無でデータフレームを分割 self.df0 = self.df[self.df[z] == 0] self.df1 = self.df[self.df[z] == 1] #df0のindex取得 tmp_df0_index = self.df0[self.ps].index.values #マッチング後のindexを入れる変数を定義(介入なし) self.df0_index = np.empty(0, dtype="int64") #マッチング後のindexを定義(介入あり) self.df1_index = self.df1[self.ps].index.values #最近傍マッチングで傾向スコアマッチング for eps1 in self.df1[self.ps]: #介入ありと介入なしの傾向スコアの差の絶対値を計算 abs_eps_list = [abs(eps1-eps0) for eps0 in self.df0[self.ps]] #傾向スコアの差の絶対値が最小のindexを取得 min_index = tmp_df0_index[abs_eps_list.index(min(abs_eps_list))] #マッチング後のindexを入れる self.df0_index = np.append(self.df0_index, min_index) def get_effect(self): #各群の結果変数の平均を取得 self.mean0 = np.average(self.df0.loc[self.df0_index][self.rs]) self.mean1 = np.average(self.df1.loc[self.df1_index][self.rs]) return (self.mean1-self.mean0) def compare_cov_mean_var(self): #共変量調整後の共変量の平均 self.cov_mean0 = np.mean(self.df0[self.cov_list].loc[self.df0_index]) self.cov_mean1 = np.mean(self.df1[self.cov_list].loc[self.df1_index]) #共変量調整前の共変量の平均 self.not_cov_mean0 = np.mean(self.df0[self.cov_list]) self.not_cov_mean1 = np.mean(self.df1[self.cov_list]) #共変量調整後の共変量の平均の差 self.adjusted_abs_diff_mean = abs(self.cov_mean1-self.cov_mean0) #共変量調整前の共変量の平均の差 self.not_adjusted_abs_diff_mean = abs(self.not_cov_mean1-self.not_cov_mean0) #作図 ##パラメータ top = np.arange(len(self.adjusted_abs_diff_mean)) height = 0.3 ##作図関数 #plt.subplot(1, 2, 1) plt.barh(top, self.adjusted_abs_diff_mean, label="adjusted", height=0.3, align="edge", zorder=10) plt.barh(top+height, self.not_adjusted_abs_diff_mean, label="not adjusted", height=0.3, align="edge", zorder=10) plt.yticks(top+height, self.adjusted_abs_diff_mean.index) plt.vlines([0.1], ymin=0, ymax=len(self.adjusted_abs_diff_mean), linestyle = "dashed") plt.ylim(0, len(self.adjusted_abs_diff_mean)) plt.xlabel("Absolute Mean Differences") plt.ylabel("Covariate Balance") plt.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=0, fontsize=8) クラスに含まれるパラメータは下記の通りです。 クラスのパラメータ 意味 df 分析対象のpandas.DataFrame。 cov_list dfに含まれる共変量のカラム名。文字列のリストで与える。 z dfに含まれる介入変数のカラム名。文字列で与える。 ps dfに含まれる傾向スコアの推定値のカラム名。文字列で与える。 rs dfに含まれる結果変数のカラム名。文字列で与える。 また、メソッドの意味は下記の通りです。 クラスのパラメータ 意味 get_effect 傾向スコアマッチングによる因果効果の推定値を返す。マッチング方法は最近傍マッチング。 compare_cov_mean_var 傾向スコアマッチングによる共変量調整前後の共変量の平均のつり合いをみる。棒グラフを返す。 このクラスを実際に使うときのスクリプトは下記のようになります。 傾向スコアマッチングの実行 #インスタンス作成 match_result = ( propensityScoreMatching( df = df_estps, #対象のデータフレーム cov_list = cov_list_str, #共変量のカラム名のリスト。cov_list_str = ["X1", "X2", "X3", "X4"] z = "treatment_flg", #介入変数のカラム名 ps = "est_propensity_score", #傾向スコアの推定値のカラム名 rs = "result" #結果変数のカラム名 ) ) #因果効果の推定量計算 match_result.get_effect() #平均の釣り合いの確認 match_result.compare_cov_mean_var() 結果は「シミュレーション結果」で報告します。 IPW推定量 IPW推定量を計算するためのclass。クラス内のメソッドや引数の意味は傾向スコアマッチングと同じです。 IPW推定量を計算するclass #IPW推定量によって効果を算出するクラス class inverseProbbilityWeight: #__init_でIPWを算出する def __init__(self, df, cov_list, z, ps, rs): self.df = df self.cov_list = cov_list self.z = z self.ps = ps self.rs = rs #介入の有無でデータフレームを分割 self.df0 = self.df[self.df[z] == 0] self.df1 = self.df[self.df[z] == 1] def ipw(self, column): ipw1 = np.sum((self.df1[self.z]*self.df1[column])/self.df1[self.ps])/np.sum(self.df1[self.z]/self.df1[self.ps]) ipw0 = np.sum(((1-self.df0[self.z])*self.df0[column])/(1-self.df0[self.ps]))/np.sum((1-self.df0[self.z])/(1-self.df0[self.ps])) return(ipw1, ipw0) def get_effect(self): #各群のIPW推定量を計算 self.ipw_mean1, self.ipw_mean0 = self.ipw(column=self.rs) return (self.ipw_mean1-self.ipw_mean0) def compare_cov_mean_var(self): #共変量調整後の共変量の平均の差の絶対値 ipw_list = [self.ipw(cov) for cov in self.cov_list] self.adjusted_abs_diff_mean = [abs(ipw1-ipw0) for ipw1, ipw0 in ipw_list] #共変量調整前の共変量の平均の差の絶対値 ##共変量調整前の共変量の平均 self.not_cov_mean0 = np.mean(self.df0[self.cov_list]) self.not_cov_mean1 = np.mean(self.df1[self.cov_list]) ##差の絶対値 self.not_adjusted_abs_diff_mean = abs(self.not_cov_mean1-self.not_cov_mean0) #共変量調整前の共変量の分散の差 self.not_adjusted_abs_diff_mean = abs(self.not_cov_mean1-self.not_cov_mean0) #作図 ##パラメータ top = np.arange(len(self.adjusted_abs_diff_mean)) height = 0.3 ##作図関数 #plt.subplot(1, 2, 1) plt.barh(top, self.adjusted_abs_diff_mean, label="adjusted", height=0.3, align="edge", zorder=10) plt.barh(top+height, self.not_adjusted_abs_diff_mean, label="not adjusted", height=0.3, align="edge", zorder=10) plt.yticks(top+height, self.cov_list) plt.vlines([0.1], ymin=0, ymax=len(self.adjusted_abs_diff_mean), linestyle = "dashed") plt.ylim(0, len(self.adjusted_abs_diff_mean)) plt.xlabel("Absolute Mean Differences") plt.ylabel("Covariate Balance") plt.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=0, fontsize=8) このクラスを実際に使うときのスクリプトは下記のようになります。 IPW推定量の計算の実行 #インスタンス作成 ipw_result = ( inverseProbbilityWeight( df = df_estps, #対象のデータフレーム cov_list = cov_list_str, #共変量のカラム名のリスト。cov_list_str = ["X1", "X2", "X3", "X4"] z = "treatment_flg", #介入変数のカラム名 ps = "est_propensity_score", #傾向スコアの推定値のカラム名 rs = "result" #結果変数のカラム名 ) ) #因果効果の推定量計算 ipw_result.get_effect() #平均の釣り合いの確認 ipw_result.compare_cov_mean_var() 結果は「シミュレーション結果」で報告します。 カーネルマッチング 　classは作らず局所回帰をするためのメソッドを作成しました。まず最初にガウシアンカーネルとイパネクニコフ2次カーネルの定義をします。 カーネルの定義 ## ガウシアンカーネル def gaussianKernel(x, x_prime, width): k = np.exp(- np.linalg.norm(x - x_prime)**2 / (2 * width**2)) return(k) ##イパネクニコフ2次カーネル def epanechinikovKernel(x, x_prime, width): t = abs(x-x_prime)/width D = 3/4*(1-t**2) if abs(t) <= 1 else 0 return(D) 引数の意味は下記です。 引数 意味 x 説明変数 x_prime 説明変数。xとx_primeの間の「距離」が大きいとカーネルは小さくなります width バンド幅 次に局所回帰をするための関数を定義します。 局所回帰の実行 def localRegression(x, resvar, exvar, treatvar, df, kernel = "gaussianKernel", width = 0.1, dim = 1): #説明変数の行列作成 X = np.ones(len(df[exvar])) for p in range(1, dim+1): xp = df[exvar]**p X = np.c_[X, xp] #カーネルが入った対角行列作成 if kernel == "gaussianKernel": kernel_list = [(1-z_ele)*gaussianKernel(x, x_ele, width) for x_ele, z_ele in zip(df[exvar], df[treatvar])] elif kernel == "epanechinikovKernel": kernel_list = [(1-z_ele)*epanechinikovKernel(x, x_ele, width) for x_ele, z_ele in zip(df[exvar], df[treatvar])] else: print("please input kernel") return() W = np.diag(kernel_list) coef = np.linalg.inv((X.T)@W@X)@(X.T)@W@df[resvar] predict = np.sum(coef*[x**p for p in range(0, dim+1)]) return(predict) 引数の意味は下記です。 引数 意味 x 説明変数 resvar 結果変数のカラム名。 exvar 傾向スコアの推定値のカラム名。 treatvar 介入変数のカラム名。 df 対照のpandas.DataFrame kernel カーネル。gaussianKernel(ガウシアンカーネル)とepanechinikovKernel(イパネクニコフ2次カーネル)が選択可能。 width バンド幅 dim 局所回帰の説明変数の最大次数 関数たちは実装できたので、実際に局所回帰してみます！ まずは介入群と対照群にデータフレームを分割します。 介入群と対照群に分割 #介入群と対照群に分割する df_estps1 = df_estps[df_estps["treatment_flg"] == 1] df_estps0 = df_estps[df_estps["treatment_flg"] == 0] 次に、訓練データとテストデータに分割します。 訓練、テストデータに分割 df_estps0_train, df_estps0_test = train_test_split(df_estps0, test_size = size) 分割出来たら汎化性能の測定のために、局所回帰のハイパーパラメータを色々変えて局所回帰を実行してみます。 今回は、 カーネル：(ガウシアンカーネル,　イパネクニコフ2次カーネル)の2パターン バンド幅：(0.1, 0.2, ..., 5.0)の50パターン 局所回帰の説明変数の最大次数：(1, 2, 3, 4, 5)の5パターン の全組み合わせで局所回帰を実行してMSEで評価しました。そのスクリプトがこちら。 汎化性能の測定 #対照群の訓練データを用いたモデルの作成とテストデータを用いた汎化性能の測定 ##パラメータの候補を設定 kernel_list = ["gaussianKernel", "epanechinikovKernel"] width_list = np.arange(0.1, 5.0, 0.1) dim_list = list(range(1, 6)) ##パラメータの全列挙パターン作成 param_list = itertools.product(kernel_list, width_list, dim_list) ##各パラメータについてMSEを計算 MSE_df = pd.DataFrame(columns=["kernel", "width", "dim", "MSE"]) for k, w, d, in param_list: ##テストデータの予測値作成 predict_list = [ localRegression(x = x_ele, resvar = "result", exvar = "est_propensity_score", treatvar = "treatment_flg", df = df_estps0_train, kernel = k, width = w, dim = d) for x_ele in df_estps0_test["est_propensity_score"] ] #テストデータの真の結果 result_list = df_estps0_test["result"] #MSEをパラメータとともに保存 tmp_MSE = ( pd.DataFrame( {"kernel" : [k], "width" : w, "dim" : d, "MSE" : mean_squared_error(predict_list, result_list) } ) ) MSE_df = MSE_df.append(tmp_MSE) 我が家のPCだと12分かかりましたが、なんとか終わりました。次にMSEが最小のパラメータを抽出します。 MSE最小のパラメータ出力 #MSEが最小のパラメータを出力 min_params_df = MSE_df[MSE_df["MSE"] == min(MSE_df["MSE"])] 今回は、 カーネル：ガウシアンカーネル バンド幅：4.9 局所回帰の説明変数の最大次数：5 のパターンのMSEが最小でした。このパラメータを使って、介入群のユーザが介入を受けなかった場合の結果変数を予測します。 予測値の計算 #選択されたモデルで介入群の介入がなかったときの結果変数を予測する predict_list = [ localRegression(x = x_ele, resvar = "result", exvar = "est_propensity_score", treatvar = "treatment_flg", df = df_estps0, kernel = min_params_df["kernel"].item(), width = min_params_df["width"].item(), dim = min_params_df["dim"].item()) for x_ele in df_estps1["est_propensity_score"] ] 最後に因果効果の推定量を計算します。 因果効果の推定値の計算 #因果効果を測定する effect = np.mean(df_estps1["result"]-predict_list) シミュレーション結果 　傾向スコアマッチング・IPW・カーネルマッチングの各手法によって得られた因果効果の推定量は下記のようになりました。 真値 セレクションバイアスが入った見せかけの効果 傾向スコアマッチング IPW カーネルマッチング 1.00 1.79 0.99 1.07 1.01 いずれの手法も真値である$1.00$に近い値となっており、セレクションバイアスを取り除くことができています。 また、傾向スコアマッチングとIPWについては、平均値の釣り合いを見ることもできます。 傾向スコアマッチングで共変量を調整した場合は下図のようになりました。 また、IPWで共変量を調整した場合は下図です。 これらの図を見ると今回の場合は、X2についてはIPWの方が共変量の釣り合いがとれているものの、他の共変量については傾向スコアマッチングの方が釣り合いがとれています。したがって、今回の場合はIPWに比べて傾向スコアマッチングを選択したほうが良いとわかります。実際に、因果効果の推定値も傾向スコアマッチングは0.99であり、IPWの1.07よりも真値の1.00に近い値となっています。 　また、カーネルマッチングについては傾向スコアマッチングとIPWと比較することができないものの、因果効果の推定値は1.01となっており真値の1.00に近い値を取っています。 まとめ 　傾向スコアマッチング・IPW・カーネルマッチングの各手法の考え方と、pythonコードをまとめました。何かのお役に立てれば幸いです。 参考文献 調査観察データの統計科学 効果検証入門 局所線形回帰を完全に理解したい ロジスティック回帰分析と傾向スコア解析 (1)式の右辺は傾向スコアがバランシングスコアであること、期待値繰り返しの公式を用いることで左辺と一致する。 ↩ なんでこんな変てこりんな推定量が出てきたのかと思ったが、IPW推定量はM推定量の考え方から導出されている。また、直感的な説明は効果検証入門に載っていた。 ↩ $V(\theta_0)$の具体的な式は書こうとしたが、煩雑になりすぎるのでやめた。 調査観察データの統計科学に記載がある。 ↩

1. この記事は 複数銘柄の日足データを銘柄ごとに週足データに算出するコードを掲載する。 2. コードの説明 pd.Grouperを使うのがポイントです。Multiindexの場合、groupbyを使用するとresampleのコードでエラーが発生します。よって、pd.Grouperを使用します。 sample.py import pandas as pd import numpy as np dat=[ [1570,"2021/7/5",15650,15670,15490,15520,3315590], [1570,"2021/7/6",15620,15690,15500,15560,3165551], [1570,"2021/7/7",15050,15340,15030,15260,6047534], [1570,"2021/7/8",15200,15280,14990,15020,4927900], [1570,"2021/7/9",14620,14880,14240,14850,15201967], [1570,"2021/7/12",15440,15500,15360,15460,8908945], [1570,"2021/7/13",15630,15780,15610,15640,6024571], [1570,"2021/7/14",15410,15610,15370,15510,4788303], [1570,"2021/7/15",15430,15480,15120,15160,6687459], [1570,"2021/7/16",14840,15090,14700,14890,8021962], [1570,"2021/7/19",14500,14650,14320,14470,8402245], [1570,"2021/7/20",14200,14400,14140,14220,8901150], [1570,"2021/7/21",14610,14730,14260,14380,7473282], [1305,"2021/7/5",2073,2076,2066,2074,155700], [1305,"2021/7/6",2080,2085,2072,2078,188890], [1305,"2021/7/7",2054,2070,2047,2058,740510], [1305,"2021/7/8",2020,2024,2007,2007,216620], [1305,"2021/7/9",1976,2004,1964,1999,403640], [1305,"2021/7/12",2035,2042,2028,2036,208860], [1305,"2021/7/13",2049,2060,2048,2051,655740], [1305,"2021/7/14",2041,2056,2038,2047,83270], [1305,"2021/7/15",2045,2045,2021,2022,132610], [1305,"2021/7/16",2014,2027,2008,2016,77810], [1305,"2021/7/19",1995,2001,1981,1991,162720], [1305,"2021/7/20",1969,1981,1964,1970,476450], [1305,"2021/7/21",2000,2009,1983,1985,91520] ] #datをDataFrame型変数dfに格納する。 df = pd.DataFrame(dat,columns=["Code","Date","Open","High","Low","Close","Volume"]) df['Date'] = pd.to_datetime(df['Date'], format='%Y/%m/%d') #文字列型を日付型に変更する print("### 銘柄番号1305,1570の日足株価データ") display(df) d_ohlcv = {'Open': 'first','High': 'max','Low': 'min','Close': 'last','Volume': 'sum'} df_stock = df.groupby(['Code', pd.Grouper(freq='W', key='Date')]).agg(d_ohlcv) #株価データ(週足)を取得する。 print("#### 銘柄番号1305,1570の週足株価データ") display(df_stock) 実行結果

フォルダーを再帰的に作りたい いままでは新学期になるたびに講義ファイルをまとめるフォルダーを作っておりました。 大した労力じゃないかもしれないけど、地味に大変だから何とかしたい！ コード フォルダーを作りたい場所で動かすと出来上がり。 はじめに一番上になるフォルダーの名前を決めることができます。 曜日や時限数は調整してみてください。 foldermaker.py import os week = ('1_Monday', '2_Tuesday', '3_Wednesday', '4_Thursday', '5_Fliday') period = ('1', '2', '3', '4', '5') title = input('input folder name...') for w in week: for p in period: os.makedirs(title + '/' + w + '/' + p) 参考

前回の記事では、CATRモデルの実装コードを、githubから落として、動作検証を行ってみました。 （ 前回の記事 ） Transformerベースの画像キャプション生成モデル CATRを動かしてみた 今回は、git cloneした推論工程の実行ファイル predict.pyを改修して、出力されたキャプション文を入力画像に文字列埋込して、結果をファイル保存するように書き換えてみました。 なお、キャプション文が40文字を超える場合は、改行して、改行した結果を入力画像に埋め込むようにしました。4行目以降は、埋め込まずに捨てることにしました。 （ 入力ファイル ） （ 出力ファイル ） （ 実行画面 ） Terminal electron@diynoMacBook-Pro catr % python3 predict_output_file_save.py --path ./image_files/test008.jpeg --v v3 Using cache found in /Users/ocean/.cache/torch/hub/saahiluppal_catr_master A man sitting on a beach with a surfboard. 2 electron@diynoMacBook-Pro catr % （ 入力ファイル ） （ 出力ファイル ） （ 実行画面 ） Terminal electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % （ 実行画面 ） Terminal electron@diynoMacBook-Pro catr %  electron@diynoMacBook-Pro catr % python3 predict_output_file_save.py --path ./image_files/test005.jpg --v v3 Using cache found in /Users/ocean/.cache/torch/hub/saahiluppal_catr_master A man holding a bottle of beer and a man in a yellow shirt. 2 electron@diynoMacBook-Pro catr % （ 入力ファイル ） （ 出力ファイル ） （ 実行画面 ） Terminal electron@diynoMacBook-Pro catr % python3 predict_output_file_save.py --path ./image_files/special.jpg --v v3 Using cache found in /Users/ocean/.cache/torch/hub/saahiluppal_catr_master A picture of a bunch of books and a teddy bear. 2 electron@diynoMacBook-Pro catr % 改修後のスクリプトファイル predict_output_file_save.py import torch from transformers import BertTokenizer from PIL import Image import argparse import cv2 import matplotlib.pyplot as plt from models import caption from datasets import coco, utils from configuration import Config import os import textwrap parser = argparse.ArgumentParser(description='Image Captioning') parser.add_argument('--path', type=str, help='path to image', required=True) parser.add_argument('--v', type=str, help='version', default='v3') parser.add_argument('--checkpoint', type=str, help='checkpoint path', default=None) args = parser.parse_args() image_path = args.path version = args.v checkpoint_path = args.checkpoint config = Config() if version == 'v1': model = torch.hub.load('saahiluppal/catr', 'v1', pretrained=True) elif version == 'v2': model = torch.hub.load('saahiluppal/catr', 'v2', pretrained=True) elif version == 'v3': model = torch.hub.load('saahiluppal/catr', 'v3', pretrained=True) else: print("Checking for checkpoint.") if checkpoint_path is None: raise NotImplementedError('No model to chose from!') else: if not os.path.exists(checkpoint_path): raise NotImplementedError('Give valid checkpoint path') print("Found checkpoint! Loading!") model,_ = caption.build_model(config) print("Loading Checkpoint...") checkpoint = torch.load(checkpoint_path, map_location='cpu') model.load_state_dict(checkpoint['model']) tokenizer = BertTokenizer.from_pretrained('bert-base-uncased') start_token = tokenizer.convert_tokens_to_ids(tokenizer._cls_token) end_token = tokenizer.convert_tokens_to_ids(tokenizer._sep_token) image = Image.open(image_path) image = coco.val_transform(image) image = image.unsqueeze(0) def create_caption_and_mask(start_token, max_length): caption_template = torch.zeros((1, max_length), dtype=torch.long) mask_template = torch.ones((1, max_length), dtype=torch.bool) caption_template[:, 0] = start_token mask_template[:, 0] = False return caption_template, mask_template caption, cap_mask = create_caption_and_mask( start_token, config.max_position_embeddings) @torch.no_grad() def evaluate(): model.eval() for i in range(config.max_position_embeddings - 1): predictions = model(image, caption, cap_mask) predictions = predictions[:, i, :] predicted_id = torch.argmax(predictions, axis=-1) if predicted_id[0] == 102: return caption caption[:, i+1] = predicted_id[0] cap_mask[:, i+1] = False return caption output = evaluate() result = tokenizer.decode(output[0].tolist(), skip_special_tokens=True) #result = tokenizer.decode(output[0], skip_special_tokens=True) final_caption_text = str(result.capitalize()) #print(result.capitalize()) print(final_caption_text) #print(result) 上と同じ文字列が出力される # Pythonで文字列を折り返し・切り詰めして整形する # https://note.nkmk.me/python-textwrap-wrap-fill-shorten/ # 生成されたキャプション文を40文字で改行する。-> List型が返却される caption_text_wrap_list = textwrap.wrap(str(final_caption_text), 40) print(len(caption_text_wrap_list)) # 入力した画像ファイルに、生成したキャプション文の文字列を埋込む img = cv2.imread(image_path) #https://qiita.com/Daiki_P/items/91acd5cc208fedd16ee9 img_RGB = cv2.cvtColor(img, cv2.COLOR_BGR2RGB) #BGRをRGBに変換. if len(caption_text_wrap_list) == 1: cv2.putText(img, str(final_caption_text), (0, 50), cv2.FONT_HERSHEY_TRIPLEX, 1, (0, 0, 255), 1, cv2.LINE_AA) elif len(caption_text_wrap_list) ==2 : cv2.putText(img, str(caption_text_wrap_list[0]), (0, 50), cv2.FONT_HERSHEY_TRIPLEX, 1, (0, 0, 255), 1, cv2.LINE_AA) cv2.putText(img, str(caption_text_wrap_list[1]), (0, 90), cv2.FONT_HERSHEY_TRIPLEX, 1, (0, 0, 255), 1, cv2.LINE_AA) elif len(caption_text_wrap_list) >= 3: cv2.putText(img, str(caption_text_wrap_list[0]), (0, 50), cv2.FONT_HERSHEY_TRIPLEX, 1, (0, 0, 255), 1, cv2.LINE_AA) cv2.putText(img, str(caption_text_wrap_list[1]), (0, 90), cv2.FONT_HERSHEY_TRIPLEX, 1, (0, 0, 255), 1, cv2.LINE_AA) cv2.putText(img, str(caption_text_wrap_list[2]), (0, 130), cv2.FONT_HERSHEY_TRIPLEX, 1, (0, 0, 255), 1, cv2.LINE_AA) #output_file_name = "result.jpeg" # ファイルパスの最期の.jpg, .png, .jpegを削除 output_file_name = "{0}_解析結果.jpeg".format(image_path[0:-1]) cv2.imwrite(output_file_name,img) #cv2.putText(img_RGB, str(final_caption_text), (0, 100), cv2.FONT_HERSHEY_TRIPLEX, 1, (0, 0, 255), 1, cv2.LINE_AA) #plt.imshow(img_RGB) #plt.show() ```

Blobに配置した画像をNotebooksで触る必要があったので備忘録として。 import os import matplotlib.pyplot as plt from azure.storage.blob import BlobServiceClient, BlobClient, ContainerClient, __version__ from PIL import Image from io import BytesIO container_name = "container" target_image_folder = "hoge/fuga" STORAGE_ACCOUNT_CONNECTION_STRING = "" blob_service_client = BlobServiceClient.from_connection_string(STORAGE_ACCOUNT_CONNECTION_STRING) with blob_service_client: container_client = blob_service_client.get_container_client(container_name) blobs_list = [] img_list = [] for blob in container_client.list_blobs(): blobs_list.append(blob) for blob in blobs_list: if os.path.dirname(blob.name) != target_image_folder: continue print("Downloading blob: ",blob.name) blob_client = blob_service_client.get_blob_client(container=container_name, blob=blob.name) stream = blob_client.download_blob() data = stream.readall() img = Image.open(BytesIO(data)) img_list.append(img) print(len(img_list)) for img in img_list: plt.imshow(img) plt.show

概要 競プロ典型90問の解説記事です。 解説の画像を見ても分からない（理解力が足りない）ことが多々あったので、後で解き直した際に確認できるようまとめました。 ※順次、全ての問題の解説記事を挙げていく予定です。 ※★５以上の問題は難易度的に後回しにしているため、投稿時期が遅くなる可能性があります。（代わりに丁寧に解説してくれる方いたらぜひお願いします） 問題002-Encyclopedia of Parentheses 問題概要 長さNの「正しい」カッコ列を辞書順で出力する。 正しいとは、 - ()のペアが揃っている。 - 右から見て行った時に、"("より")"の方が多くなってはいけない - ()以外の文字列が含まれてはいけない 制約 1 <= N <= 20 \\ Nは整数 解き方 Nが小さく、文字列のシンプルさから、工夫して全探索で求められそう。 一応、同じものを含んだ並び替えで全パターン列挙することも考えたが、パターンが最大で $20!$ となってしまうため、 しんどそう。（厳密には、$20!/(10!*10!)$） そこで今回は、 1. "(" を1, ")" を0とみなして、ビット全探索 2. １で出てきた文字列において "(" より先に ")" が多く出現しないか確認 といった流れで実装していく。 １に関して、0か1かの２通り × N文字 なので、$O(2^N)$の計算量。 （ビット全探索に関してはこちら） 2に関して、N文字を左から検索するだけなので、$O(N)$。 2^Nのパターンにおいて、N回の検索を行うので、 合計で$O(N*2^N)$。これは、10^7のオーダーになるので、十分に間に合う。 引用元：競プロ典型90問 Github 実装 answer.py # 入力の受け取り N = int(input()) # カッコ列を格納する配列 pare = [] # ビット全探索 for i in range(1 << N): # 1の時 "(" を、0の時 ")" を文字列に追加していく l = "" for j in range(N): if i >> j & 1: l += "(" else: l += ")" # 文字列の長さがNの時のみ、左から検索する。その時cntで"("と")"の出現回数を計算する。cntが0より小さくならなかった＋cntが0(出現回数が等しかった)文字列のみをpareに格納する if len(l) == N: cnt = 0 flag = True for k in range(N): if l[k] == "(": cnt += 1 else: cnt -= 1 if cnt < 0: flag = False if flag and cnt == 0: pare.append(l) # pareを辞書順にソート、順番に出力する pare.sort() for i in range(len(pare)): print(pare[i]) ※ちなみに、上記のコードは、Pythonで提出するとTLEになり、pypyだとACになります。 やたらと長い配列や再起関数を用いる場合を除いては、pypyの方が基本的に早いので、提出はpypyがおすすめです。 参考記事 最後に 問題の解説を読んでも他の人のコードを見てもさっぱり分からないという方の 力に少しでもなれれば幸いです。 コードの改善点やその他、ご意見などあれば、気軽にコメントください。 ここまでお読みいただきありがとうございました。

与えられた10分間の観測データを使い、10分後のRealized Volatilityを予測するコンペ Optiver Realized Volatility Predictにて知っておくべき周辺知識をまとめた。 Realized Volatilityとは ボラティリティの推定値の一種 ボラティリティは資産価格の変動度合いを表し、この指標が高い商品は高リスクと判定できる。 Realized Volatilityは、高頻度データによる株価リターンを使用し、日時換算したボラティリティを指す。 ごちゃごちゃ書いているが、要は日中に観測されたリターンの2乗を観測回数分足して、 日時換算したものがRealized Volatility 定義式はt日の観測回数i回のRVが下記式 $$ RV_t = \sum_{i=1}^{n}r^2_{-1+i/n}\quad $$ $\lim_{n \to \infty}RV_t$が一日単位で観測された累積ボラティリティと 一致していることが示されているので、 観測回数を増やせば増やすほど真のボラティリティに近づく。 なお、一般的なボラティリティ自体は$dS_{t}=\mu S_{t}\, dt+\sigma S_{t}\, dB_{t} $ におけるシグマを指すようです。（Wikipedia定義） 名前がややこしいが、あくまでRVは統計量なので平均値の仲間のようなもの（と思っている） 今回のコンペにおけるRealized Volatility 上記は金融工学におけるやや厳格なRealized Volatilityを記載した。 今回のコンペでは金融知識を持たないKagglerに配慮し、もっとシンプルなRealized Volatilityの定義が用いられている。 本題に入る前に、下準備とし市場の流動性と株式評価の上で重要なOrder book statisticsとLog Returnについて確認する。 bid/ask spread Weighted averaged price Log Return bid/ask spread 原文 As different stocks trade on different level on the market we take the ratio of best offer price and best bid price to calculate the bid-ask spread. 異なる水準の株価のトレードをする際の最高買値と最安売値の比であり、以下で定義される。 $$ BidAskSpread=BestOffer/BestBid-1 $$ 株の流動性の1つの指標として用いられようで、これが1に近いと取引費用は小さく 流動性が高くなります。（つまりボラは低い） Weighted averaged price 読んで字のごとく重みづけされた平均価格。定義見たほうが早い $$ \frac{WAP=BidPrice*AskSize+AskPrice*BidSize}{BidSize+AskSize} $$ Realized Volatilityを計算する上での価格はこれを用いるので、 関数化必須 Log returns 異なる時系列の株価を比較するときに、logを使うと計算が楽になります（あたりまえ体操） $$ r_{t_1,t_2} = \log(\frac{S_{t_2}}{S_{t_1}}) $$ Realized Volatility ようやくRealized Volatilityの定義に入る。 年単位で正規化された株のログリターンの、年換算標準偏差がVolatility WAPによって定義された株価$S_t$をログリターンにぶち込み、 各tにおける二乗和の平方根を取った値σこそ、今回求めるRealized Volatility $$ σ = \sqrt{\sum_{t}r^2_{t-1,t}} $$ 頑張ってこれを予測していきましょう 終わりに RMSPEについて 指標としてはRMSPE（Root Mean Squared Percentage Error）が使われる。 誤差を正解値で割った値の二乗の平均の平方根（長い） RMSPEを0に近づけていくことがコンペの目標だが、RMSPEだからと言って留意すべきことは特段なさそう 特徴量の検討方法について考察できたら追記します。 参考文献 Kaggleのチュートリアル https://www.kaggle.com/jiashenliu/introduction-to-financial-concepts-and-data 高頻度データによるボラティリティの推定：Realized Volatilityのサーベイと日本の株価指数および株価指数先物の実証分析 https://www1.econ.hit-u.ac.jp/finmodel/pdf-2/2014-7Wakaki.pdf Realized volatility のショックの影響力についての検証 https://www.imes.boj.or.jp/research/papers/japanese/07-J-14.pdf

動かしたコード CAption TRansformerのPyTorch実装がありましたので、git cloneして動かしてみた。 CPTR: FULL TRANSFORMER NETWORK FOR IMAGE CAPTIONING GitHub saahiluppal / catr 動かしてみた git cloneすると、学習済みのpre-trained modelが同梱されています・ 今回は、学習済みモデルに、画像ファイルを食わせて、その画像ファイルに対応する説明文（キャップション文）を出力させてみます。 学習済みモデルを使うモデル推論工程は、predict.py*を使います。 食わせる画像ファイルは、コマンドライン引数で、画像のパスを渡します。 Terminal electron@diynoMacBook-Pro catr % cat predict.py | grep image_path image_path = args.path image = Image.open(image_path) electron@diynoMacBook-Pro catr % Google画像検索で、画像ファイルを取ってくる Terminal electron@diynoMacBook-Pro catr % mkdir image_files electron@diynoMacBook-Pro catr % cd image_files electron@diynoMacBook-Pro image_files % ls test001.png electron@diynoMacBook-Pro image_files % 取得した画像ファイルを渡して、キャプション文を出してみる 学習済みモデルは、v1, v2, v3を選べますが、まずはv2を選びます。 実行すると、v2モデルの学習済みの重みパラメータファイルのダウンロードが始まります。 Terminal electron@diynoMacBook-Pro catr % python3 predict.py --path ./image_files/test001.png --v v2 Downloading: "https://github.com/saahiluppal/catr/archive/master.zip" to /Users/electron/.cache/torch/hub/master.zip Downloading: "https://download.pytorch.org/models/resnet101-5d3b4d8f.pth" to /Users/electron/.cache/torch/hub/checkpoints/resnet101-5d3b4d8f.pth 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 170M/170M [00:13<00:00, 13.6MB/s] Downloading: "https://github.com/saahiluppal/catr/releases/download/0.2/weight389123791.pth" to /Users/electron/.cache/torch/hub/checkpoints/weight389123791.pth 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 322M/322M [00:26<00:00, 13.0MB/s] Downloading: 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████| 232k/232k [00:00<00:00, 792kB/s] Downloading: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████| 28.0/28.0 [00:00<00:00, 7.09kB/s] Downloading: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████| 466k/466k [00:00<00:00, 1.06MB/s] Downloading: 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 570/570 [00:00<00:00, 144kB/s] Traceback (most recent call last): Using cache found in /Users/electron/.cache/torch/hub/saahiluppal_catr_master A man riding on the back of a brown horse. electron@diynoMacBook-Pro catr % [結果] A man riding on the back of a brown horse. A man riding on the back of a brown horse.という文が、Terminalに標準出力されました。 なお、v2を最初に選択したときに、v2モデルの学習済み重みファイルは、torch.hubからダウンロードされます。 predict.py if version == 'v1': model = torch.hub.load('saahiluppal/catr', 'v1', pretrained=True) elif version == 'v2': model = torch.hub.load('saahiluppal/catr', 'v2', pretrained=True) elif version == 'v3': model = torch.hub.load('saahiluppal/catr', 'v3', pretrained=True) 他の画像で試してみます。 Terminal electron@diynoMacBook-Pro catr % ls ./image_files test001.png test002.jpg electron@diynoMacBook-Pro catr % Terminal electron@diynoMacBook-Pro catr % python3 predict.py --path ./image_files/test002.jpg --v v2 Using cache found in /Users/electron/.cache/torch/hub/saahiluppal_catr_master A group of men laying on top of a bed. electron@diynoMacBook-Pro catr % [結果] A group of men laying on top of a bed. 今度は、A group of men laying on top of a bed.という文が返されました。 学習済みモデルのv3を使ってみる Terminal electron@diynoMacBook-Pro catr % python3 predict.py --path ./image_files/test002.jpg --v v3 Using cache found in /Users/electron/.cache/torch/hub/saahiluppal_catr_master Downloading: "https://github.com/saahiluppal/catr/releases/download/0.2/weight493084032.pth" to /Users/electron/.cache/torch/hub/checkpoints/weight493084032.pth 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 322M/322M [00:24<00:00, 13.9MB/s] A group of men laying on top of a beach. electron@diynoMacBook-Pro catr % 同じ文が返されました。 v1ではどうでしょうか。 Terminal electron@diynoMacBook-Pro catr % python3 predict.py --path ./image_files/test002.jpg --v v1 Using cache found in /Users/electron/.cache/torch/hub/saahiluppal_catr_master Downloading: "https://github.com/saahiluppal/catr/releases/download/0.1/weights_9348032.pth" to /Users/electron/.cache/torch/hub/checkpoints/weights_9348032.pth 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 322M/322M [00:23<00:00, 14.4MB/s] A man in a wetsuit laying on a beach. electron@diynoMacBook-Pro catr % [結果] A man in a wetsuit laying on a beach. 今度は、違う文が返されました。 返された文は、A man in a wetsuit laying on a beach.です。 v4はないみたいです。 Terminal electron@diynoMacBook-Pro catr % python3 predict.py --path ./image_files/test002.jpg --v v4 Checking for checkpoint. Traceback (most recent call last): File "/Users/electron/Desktop/catr/predict.py", line 32, in <module> raise NotImplementedError('No model to chose from!') NotImplementedError: No model to chose from! electron@diynoMacBook-Pro catr % 最初の画像を、v3モデルに与えてみます。 Terminal electron@diynoMacBook-Pro catr % python3 predict.py --path ./image_files/test001.png --v v3 Using cache found in /Users/electron/.cache/torch/hub/saahiluppal_catr_master A man is riding a horse in a field electron@diynoMacBook-Pro catr % [結果] A man is riding a horse in a field 最初のv2モデルとは違う文が提案されました。 なお、今度は最期にピリオドが出力されていないのが面白いです。 A man is riding a horse in a field 第３の画像を与えてみる Terminal electron@diynoMacBook-Pro catr % ls ./image_files test001.png test002.jpg test003.jpg electron@diynoMacBook-Pro catr % v3モデルを選びます。 Terminal electron@diynoMacBook-Pro catr % python3 predict.py --path ./image_files/test003.jpg --v v3 Using cache found in /Users/electron/.cache/torch/hub/saahiluppal_catr_master A group of young people posing for a picture. electron@diynoMacBook-Pro catr % [結果] A group of young people posing for a picture. 以下の文が返されました。 A group of young people posing for a picture. 環境構築 Terminal electron@diynoMacBook-Pro Desktop % git clone https://github.com/saahiluppal/catr.git catrディレクトリの中身 Terminal electron@diynoMacBook-Pro Desktop % ls catr LICENSE catr_demo.ipynb datasets finetune.py main.py predict.py README.md configuration.py engine.py hubconf.py models requirements.txt electron@diynoMacBook-Pro Desktop % requirement.txtに列挙されているImportするモジュール一覧 Terminal electron@diynoMacBook-Pro Desktop % cd catr electron@diynoMacBook-Pro catr % cat requirements.txt torch torchvision numpy transformers tqdm% electron@diynoMacBook-Pro catr % importする Terminal electron@diynoMacBook-Pro catr % pip3 install -r requirements.txt （ 省略 ） Successfully installed click-8.0.1 huggingface-hub-0.0.12 joblib-1.0.1 packaging-21.0 regex-2021.8.3 sacremoses-0.0.45 tokenizers-0.10.3 transformers-4.9.1 WARNING: You are using pip version 21.2.1; however, version 21.2.2 is available. You should consider upgrading via the '/usr/local/opt/python@3.9/bin/python3.9 -m pip install --upgrade pip' command. electron@diynoMacBook-Pro catr % Transformerを用いた画像キャプション文生成モデルということで、 Transformerの実装コードとして、Huggingfaceのものが入った。 Terminal electron@diynoMacBook-Pro catr % ls LICENSE catr_demo.ipynb datasets finetune.py main.py predict.py README.md configuration.py engine.py hubconf.py models requirements.txt electron@diynoMacBook-Pro catr % 今回は、自前の学習データ（画像とキャプション文のペア集合）を学ばせるモデル学習は行わない。 すでに学習済みのpre-trainedモデルに対象画像を与えて、入力した画像を説明するキャプション文を出力させる「推論」工程だけを行います。 「推論」工程は、predict.pyを使います。 Pythonのどのモジュールをimportしているか、確認します。 Terminal electron@diynoMacBook-Pro catr % head -50 predict.py import torch from transformers import BertTokenizer from PIL import Image import argparse from models import caption from datasets import coco, utils from configuration import Config import os parser = argparse.ArgumentParser(description='Image Captioning') parser.add_argument('--path', type=str, help='path to image', required=True) parser.add_argument('--v', type=str, help='version', default='v3') parser.add_argument('--checkpoint', type=str, help='checkpoint path', default=None) args = parser.parse_args() image_path = args.path version = args.v checkpoint_path = args.checkpoint config = Config() if version == 'v1': model = torch.hub.load('saahiluppal/catr', 'v1', pretrained=True) elif version == 'v2': model = torch.hub.load('saahiluppal/catr', 'v2', pretrained=True) elif version == 'v3': model = torch.hub.load('saahiluppal/catr', 'v3', pretrained=True) else: print("Checking for checkpoint.") if checkpoint_path is None: raise NotImplementedError('No model to chose from!') else: if not os.path.exists(checkpoint_path): raise NotImplementedError('Give valid checkpoint path') print("Found checkpoint! Loading!") model,_ = caption.build_model(config) print("Loading Checkpoint...") checkpoint = torch.load(checkpoint_path, map_location='cpu') model.load_state_dict(checkpoint['model']) tokenizer = BertTokenizer.from_pretrained('bert-base-uncased') start_token = tokenizer.convert_tokens_to_ids(tokenizer._cls_token) end_token = tokenizer.convert_tokens_to_ids(tokenizer._sep_token) image = Image.open(image_path) image = coco.val_transform(image) image = image.unsqueeze(0) electron@diynoMacBook-Pro catr % 以下のmodels, datasets, configuratoon*は、git cloneした同梱されている資材と思われる。 predict.pyを別のディレクトリに移動したり、別のスクリプトファイルにコードを移植する場合は、これらの資材も、そのスクリプトファイルの格納先と同じ場所に置く必要がありそうです。 from models import caption from datasets import coco, utils from configuration import Config 次に、対象となる画像ファイルをどこで読み込んでいるのか、を見ていきます。 Terminal electron@diynoMacBook-Pro catr % cat predict.py | grep image parser.add_argument('--path', type=str, help='path to image', required=True) image_path = args.path image = Image.open(image_path) image = coco.val_transform(image) image = image.unsqueeze(0) predictions = model(image, caption, cap_mask) electron@diynoMacBook-Pro catr % argparseを使って、コマンドライン引数で、画像のパスを受け取っていることが分かりました。 Terminal electron@diynoMacBook-Pro catr % cat predict.py | grep image_path image_path = args.path image = Image.open(image_path) electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % mkdir image_files electron@diynoMacBook-Pro catr % cd image_files electron@diynoMacBook-Pro image_files % ls test001.png electron@diynoMacBook-Pro image_files % electron@diynoMacBook-Pro image_files % cd .. electron@diynoMacBook-Pro catr % cat predict.py | grep args args = parser.parse_args() image_path = args.path version = args.v checkpoint_path = args.checkpoint electron@diynoMacBook-Pro catr % cat predict.py | grep argument parser.add_argument('--path', type=str, help='path to image', required=True) parser.add_argument('--v', type=str, help='version', default='v3') parser.add_argument('--checkpoint', type=str, help='checkpoint path', default=None) electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % python3 predict.py --path /image_files/test001.png --v v2 Downloading: "https://github.com/saahiluppal/catr/archive/master.zip" to /Users/electron/.cache/torch/hub/master.zip Downloading: "https://download.pytorch.org/models/resnet101-5d3b4d8f.pth" to /Users/electron/.cache/torch/hub/checkpoints/resnet101-5d3b4d8f.pth 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 170M/170M [00:13<00:00, 13.6MB/s] Downloading: "https://github.com/saahiluppal/catr/releases/download/0.2/weight389123791.pth" to /Users/electron/.cache/torch/hub/checkpoints/weight389123791.pth 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 322M/322M [00:26<00:00, 13.0MB/s] Downloading: 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████| 232k/232k [00:00<00:00, 792kB/s] Downloading: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████| 28.0/28.0 [00:00<00:00, 7.09kB/s] Downloading: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████| 466k/466k [00:00<00:00, 1.06MB/s] Downloading: 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 570/570 [00:00<00:00, 144kB/s] Traceback (most recent call last): File "/Users/electron/Desktop/catr/predict.py", line 46, in <module> image = Image.open(image_path) File "/usr/local/lib/python3.9/site-packages/PIL/Image.py", line 2968, in open fp = builtins.open(filename, "rb") FileNotFoundError: [Errno 2] No such file or directory: '/image_files/test001.png' electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % python3 predict.py --path ./image_files/test001.png --v v2 Using cache found in /Users/electron/.cache/torch/hub/saahiluppal_catr_master A man riding on the back of a brown horse. electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % ls ./image_files/ test001.png electron@diynoMacBook-Pro catr % ls /image_files/ ls: /image_files/: No such file or directory electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % ls /image_files/test001.png ls: /image_files/test001.png: No such file or directory electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % ls ./image_files/test001.png ./image_files/test001.png electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % ls LICENSE catr_demo.ipynb engine.py image_files predict.py README.md configuration.py finetune.py main.py requirements.txt __pycache__ datasets hubconf.py models electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % ls ./image_files test001.png test002.jpg electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % python3 predict.py --path /image_files/test002.jpg --v v2 Using cache found in /Users/electron/.cache/torch/hub/saahiluppal_catr_master Traceback (most recent call last): File "/Users/electron/Desktop/catr/predict.py", line 46, in <module> image = Image.open(image_path) File "/usr/local/lib/python3.9/site-packages/PIL/Image.py", line 2968, in open fp = builtins.open(filename, "rb") FileNotFoundError: [Errno 2] No such file or directory: '/image_files/test002.jpg' electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % python3 predict.py --path ./image_files/test002.jpg --v v2 Using cache found in /Users/electron/.cache/torch/hub/saahiluppal_catr_master A group of men laying on top of a bed. electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % python3 predict.py --path ./image_files/test002.jpg --v v3 Using cache found in /Users/electron/.cache/torch/hub/saahiluppal_catr_master Downloading: "https://github.com/saahiluppal/catr/releases/download/0.2/weight493084032.pth" to /Users/electron/.cache/torch/hub/checkpoints/weight493084032.pth 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 322M/322M [00:24<00:00, 13.9MB/s] A group of men laying on top of a beach. electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % python3 predict.py --path ./image_files/test002.jpg --v v1 Using cache found in /Users/electron/.cache/torch/hub/saahiluppal_catr_master Downloading: "https://github.com/saahiluppal/catr/releases/download/0.1/weights_9348032.pth" to /Users/electron/.cache/torch/hub/checkpoints/weights_9348032.pth 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 322M/322M [00:23<00:00, 14.4MB/s] A man in a wetsuit laying on a beach. electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % python3 predict.py --path ./image_files/test002.jpg --v v4 Checking for checkpoint. Traceback (most recent call last): File "/Users/electron/Desktop/catr/predict.py", line 32, in <module> raise NotImplementedError('No model to chose from!') NotImplementedError: No model to chose from! electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % python3 predict.py --path ./image_files/test001.png --v v3 Using cache found in /Users/electron/.cache/torch/hub/saahiluppal_catr_master A man is riding a horse in a field electron@diynoMacBook-Pro catr % ls ./image_files test001.png test002.jpg test003.jpg electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % python3 predict.py --path ./image_files/test003.jpg --v v3 Using cache found in /Users/electron/.cache/torch/hub/saahiluppal_catr_master A group of young people posing for a picture. electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % cat predict.py | import zsh: command not found: import electron@diynoMacBook-Pro catr % cat predict.py | grep import import torch from transformers import BertTokenizer from PIL import Image import argparse from models import caption from datasets import coco, utils from configuration import Config import os electron@diynoMacBook-Pro catr % ls LICENSE catr_demo.ipynb engine.py image_files predict.py README.md configuration.py finetune.py main.py requirements.txt __pycache__ datasets hubconf.py models electron@diynoMacBook-Pro catr %

動かしたコード CAption TRansformerのPyTorch実装がありましたので、git cloneして動かしてみた。 CPTR: FULL TRANSFORMER NETWORK FOR IMAGE CAPTIONING GitHub saahiluppal / catr 動かしてみた git cloneすると、学習済みのpre-trained modelが同梱されています・ 今回は、学習済みモデルに、画像ファイルを食わせて、その画像ファイルに対応する説明文（キャップション文）を出力させてみます。 学習済みモデルを使うモデル推論工程は、predict.py*を使います。 食わせる画像ファイルは、コマンドライン引数で、画像のパスを渡します。 Terminal electron@diynoMacBook-Pro catr % cat predict.py | grep image_path image_path = args.path image = Image.open(image_path) electron@diynoMacBook-Pro catr % Google画像検索で、画像ファイルを取ってくる Terminal electron@diynoMacBook-Pro catr % mkdir image_files electron@diynoMacBook-Pro catr % cd image_files electron@diynoMacBook-Pro image_files % ls test001.png electron@diynoMacBook-Pro image_files % 取得した画像ファイルを渡して、キャプション文を出してみる 学習済みモデルは、v1, v2, v3を選べますが、まずはv2を選びます。 実行すると、v2モデルの学習済みの重みパラメータファイルのダウンロードが始まります。 Terminal electron@diynoMacBook-Pro catr % python3 predict.py --path ./image_files/test001.png --v v2 Downloading: "https://github.com/saahiluppal/catr/archive/master.zip" to /Users/electron/.cache/torch/hub/master.zip Downloading: "https://download.pytorch.org/models/resnet101-5d3b4d8f.pth" to /Users/electron/.cache/torch/hub/checkpoints/resnet101-5d3b4d8f.pth 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 170M/170M [00:13<00:00, 13.6MB/s] Downloading: "https://github.com/saahiluppal/catr/releases/download/0.2/weight389123791.pth" to /Users/electron/.cache/torch/hub/checkpoints/weight389123791.pth 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 322M/322M [00:26<00:00, 13.0MB/s] Downloading: 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████| 232k/232k [00:00<00:00, 792kB/s] Downloading: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████| 28.0/28.0 [00:00<00:00, 7.09kB/s] Downloading: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████| 466k/466k [00:00<00:00, 1.06MB/s] Downloading: 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 570/570 [00:00<00:00, 144kB/s] Traceback (most recent call last): Using cache found in /Users/electron/.cache/torch/hub/saahiluppal_catr_master A man riding on the back of a brown horse. electron@diynoMacBook-Pro catr % [結果] A man riding on the back of a brown horse. A man riding on the back of a brown horse.という文が、Terminalに標準出力されました。 なお、v2を最初に選択したときに、v2モデルの学習済み重みファイルは、torch.hubからダウンロードされます。 predict.py if version == 'v1': model = torch.hub.load('saahiluppal/catr', 'v1', pretrained=True) elif version == 'v2': model = torch.hub.load('saahiluppal/catr', 'v2', pretrained=True) elif version == 'v3': model = torch.hub.load('saahiluppal/catr', 'v3', pretrained=True) 他の画像で試してみます。 Terminal electron@diynoMacBook-Pro catr % ls ./image_files test001.png test002.jpg electron@diynoMacBook-Pro catr % Terminal electron@diynoMacBook-Pro catr % python3 predict.py --path ./image_files/test002.jpg --v v2 Using cache found in /Users/electron/.cache/torch/hub/saahiluppal_catr_master A group of men laying on top of a bed. electron@diynoMacBook-Pro catr % [結果] A group of men laying on top of a bed. 今度は、A group of men laying on top of a bed.という文が返されました。 学習済みモデルのv3を使ってみる Terminal electron@diynoMacBook-Pro catr % python3 predict.py --path ./image_files/test002.jpg --v v3 Using cache found in /Users/electron/.cache/torch/hub/saahiluppal_catr_master Downloading: "https://github.com/saahiluppal/catr/releases/download/0.2/weight493084032.pth" to /Users/electron/.cache/torch/hub/checkpoints/weight493084032.pth 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 322M/322M [00:24<00:00, 13.9MB/s] A group of men laying on top of a beach. electron@diynoMacBook-Pro catr % 同じ文が返されました。 v1ではどうでしょうか。 Terminal electron@diynoMacBook-Pro catr % python3 predict.py --path ./image_files/test002.jpg --v v1 Using cache found in /Users/electron/.cache/torch/hub/saahiluppal_catr_master Downloading: "https://github.com/saahiluppal/catr/releases/download/0.1/weights_9348032.pth" to /Users/electron/.cache/torch/hub/checkpoints/weights_9348032.pth 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 322M/322M [00:23<00:00, 14.4MB/s] A man in a wetsuit laying on a beach. electron@diynoMacBook-Pro catr % [結果] A man in a wetsuit laying on a beach. 今度は、違う文が返されました。 返された文は、A man in a wetsuit laying on a beach.です。 v4はないみたいです。 Terminal electron@diynoMacBook-Pro catr % python3 predict.py --path ./image_files/test002.jpg --v v4 Checking for checkpoint. Traceback (most recent call last): File "/Users/electron/Desktop/catr/predict.py", line 32, in <module> raise NotImplementedError('No model to chose from!') NotImplementedError: No model to chose from! electron@diynoMacBook-Pro catr % 最初の画像を、v3モデルに与えてみます。 Terminal electron@diynoMacBook-Pro catr % python3 predict.py --path ./image_files/test001.png --v v3 Using cache found in /Users/electron/.cache/torch/hub/saahiluppal_catr_master A man is riding a horse in a field electron@diynoMacBook-Pro catr % [結果] A man is riding a horse in a field 最初のv2モデルとは違う文が提案されました。 なお、今度は最期にピリオドが出力されていないのが面白いです。 A man is riding a horse in a field 第３の画像を与えてみる Terminal electron@diynoMacBook-Pro catr % ls ./image_files test001.png test002.jpg test003.jpg electron@diynoMacBook-Pro catr % v3モデルを選びます。 Terminal electron@diynoMacBook-Pro catr % python3 predict.py --path ./image_files/test003.jpg --v v3 Using cache found in /Users/electron/.cache/torch/hub/saahiluppal_catr_master A group of young people posing for a picture. electron@diynoMacBook-Pro catr % [結果] A group of young people posing for a picture. 以下の文が返されました。 A group of young people posing for a picture. 環境構築 Terminal electron@diynoMacBook-Pro Desktop % git clone https://github.com/saahiluppal/catr.git catrディレクトリの中身 Terminal electron@diynoMacBook-Pro Desktop % ls catr LICENSE catr_demo.ipynb datasets finetune.py main.py predict.py README.md configuration.py engine.py hubconf.py models requirements.txt electron@diynoMacBook-Pro Desktop % requirement.txtに列挙されているImportするモジュール一覧 Terminal electron@diynoMacBook-Pro Desktop % cd catr electron@diynoMacBook-Pro catr % cat requirements.txt torch torchvision numpy transformers tqdm% electron@diynoMacBook-Pro catr % importする Terminal electron@diynoMacBook-Pro catr % pip3 install -r requirements.txt （ 省略 ） Successfully installed click-8.0.1 huggingface-hub-0.0.12 joblib-1.0.1 packaging-21.0 regex-2021.8.3 sacremoses-0.0.45 tokenizers-0.10.3 transformers-4.9.1 WARNING: You are using pip version 21.2.1; however, version 21.2.2 is available. You should consider upgrading via the '/usr/local/opt/python@3.9/bin/python3.9 -m pip install --upgrade pip' command. electron@diynoMacBook-Pro catr % Transformerを用いた画像キャプション文生成モデルということで、 Transformerの実装コードとして、Huggingfaceのものが入った。 Terminal electron@diynoMacBook-Pro catr % ls LICENSE catr_demo.ipynb datasets finetune.py main.py predict.py README.md configuration.py engine.py hubconf.py models requirements.txt electron@diynoMacBook-Pro catr % 今回は、自前の学習データ（画像とキャプション文のペア集合）を学ばせるモデル学習は行わない。 すでに学習済みのpre-trainedモデルに対象画像を与えて、入力した画像を説明するキャプション文を出力させる「推論」工程だけを行います。 「推論」工程は、predict.pyを使います。 Pythonのどのモジュールをimportしているか、確認します。 Terminal electron@diynoMacBook-Pro catr % head -50 predict.py import torch from transformers import BertTokenizer from PIL import Image import argparse from models import caption from datasets import coco, utils from configuration import Config import os parser = argparse.ArgumentParser(description='Image Captioning') parser.add_argument('--path', type=str, help='path to image', required=True) parser.add_argument('--v', type=str, help='version', default='v3') parser.add_argument('--checkpoint', type=str, help='checkpoint path', default=None) args = parser.parse_args() image_path = args.path version = args.v checkpoint_path = args.checkpoint config = Config() if version == 'v1': model = torch.hub.load('saahiluppal/catr', 'v1', pretrained=True) elif version == 'v2': model = torch.hub.load('saahiluppal/catr', 'v2', pretrained=True) elif version == 'v3': model = torch.hub.load('saahiluppal/catr', 'v3', pretrained=True) else: print("Checking for checkpoint.") if checkpoint_path is None: raise NotImplementedError('No model to chose from!') else: if not os.path.exists(checkpoint_path): raise NotImplementedError('Give valid checkpoint path') print("Found checkpoint! Loading!") model,_ = caption.build_model(config) print("Loading Checkpoint...") checkpoint = torch.load(checkpoint_path, map_location='cpu') model.load_state_dict(checkpoint['model']) tokenizer = BertTokenizer.from_pretrained('bert-base-uncased') start_token = tokenizer.convert_tokens_to_ids(tokenizer._cls_token) end_token = tokenizer.convert_tokens_to_ids(tokenizer._sep_token) image = Image.open(image_path) image = coco.val_transform(image) image = image.unsqueeze(0) electron@diynoMacBook-Pro catr % 以下のmodels, datasets, configuratoon*は、git cloneした同梱されている資材と思われる。 predict.pyを別のディレクトリに移動したり、別のスクリプトファイルにコードを移植する場合は、これらの資材も、そのスクリプトファイルの格納先と同じ場所に置く必要がありそうです。 from models import caption from datasets import coco, utils from configuration import Config 次に、対象となる画像ファイルをどこで読み込んでいるのか、を見ていきます。 Terminal electron@diynoMacBook-Pro catr % cat predict.py | grep image parser.add_argument('--path', type=str, help='path to image', required=True) image_path = args.path image = Image.open(image_path) image = coco.val_transform(image) image = image.unsqueeze(0) predictions = model(image, caption, cap_mask) electron@diynoMacBook-Pro catr % argparseを使って、コマンドライン引数で、画像のパスを受け取っていることが分かりました。 Terminal electron@diynoMacBook-Pro catr % cat predict.py | grep image_path image_path = args.path image = Image.open(image_path) electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % mkdir image_files electron@diynoMacBook-Pro catr % cd image_files electron@diynoMacBook-Pro image_files % ls test001.png electron@diynoMacBook-Pro image_files % electron@diynoMacBook-Pro image_files % cd .. electron@diynoMacBook-Pro catr % cat predict.py | grep args args = parser.parse_args() image_path = args.path version = args.v checkpoint_path = args.checkpoint electron@diynoMacBook-Pro catr % cat predict.py | grep argument parser.add_argument('--path', type=str, help='path to image', required=True) parser.add_argument('--v', type=str, help='version', default='v3') parser.add_argument('--checkpoint', type=str, help='checkpoint path', default=None) electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % python3 predict.py --path /image_files/test001.png --v v2 Downloading: "https://github.com/saahiluppal/catr/archive/master.zip" to /Users/electron/.cache/torch/hub/master.zip Downloading: "https://download.pytorch.org/models/resnet101-5d3b4d8f.pth" to /Users/electron/.cache/torch/hub/checkpoints/resnet101-5d3b4d8f.pth 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 170M/170M [00:13<00:00, 13.6MB/s] Downloading: "https://github.com/saahiluppal/catr/releases/download/0.2/weight389123791.pth" to /Users/electron/.cache/torch/hub/checkpoints/weight389123791.pth 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 322M/322M [00:26<00:00, 13.0MB/s] Downloading: 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████| 232k/232k [00:00<00:00, 792kB/s] Downloading: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████| 28.0/28.0 [00:00<00:00, 7.09kB/s] Downloading: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████| 466k/466k [00:00<00:00, 1.06MB/s] Downloading: 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 570/570 [00:00<00:00, 144kB/s] Traceback (most recent call last): File "/Users/electron/Desktop/catr/predict.py", line 46, in <module> image = Image.open(image_path) File "/usr/local/lib/python3.9/site-packages/PIL/Image.py", line 2968, in open fp = builtins.open(filename, "rb") FileNotFoundError: [Errno 2] No such file or directory: '/image_files/test001.png' electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % python3 predict.py --path ./image_files/test001.png --v v2 Using cache found in /Users/electron/.cache/torch/hub/saahiluppal_catr_master A man riding on the back of a brown horse. electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % ls ./image_files/ test001.png electron@diynoMacBook-Pro catr % ls /image_files/ ls: /image_files/: No such file or directory electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % ls /image_files/test001.png ls: /image_files/test001.png: No such file or directory electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % ls ./image_files/test001.png ./image_files/test001.png electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % ls LICENSE catr_demo.ipynb engine.py image_files predict.py README.md configuration.py finetune.py main.py requirements.txt __pycache__ datasets hubconf.py models electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % ls ./image_files test001.png test002.jpg electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % python3 predict.py --path /image_files/test002.jpg --v v2 Using cache found in /Users/electron/.cache/torch/hub/saahiluppal_catr_master Traceback (most recent call last): File "/Users/electron/Desktop/catr/predict.py", line 46, in <module> image = Image.open(image_path) File "/usr/local/lib/python3.9/site-packages/PIL/Image.py", line 2968, in open fp = builtins.open(filename, "rb") FileNotFoundError: [Errno 2] No such file or directory: '/image_files/test002.jpg' electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % python3 predict.py --path ./image_files/test002.jpg --v v2 Using cache found in /Users/electron/.cache/torch/hub/saahiluppal_catr_master A group of men laying on top of a bed. electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % python3 predict.py --path ./image_files/test002.jpg --v v3 Using cache found in /Users/electron/.cache/torch/hub/saahiluppal_catr_master Downloading: "https://github.com/saahiluppal/catr/releases/download/0.2/weight493084032.pth" to /Users/electron/.cache/torch/hub/checkpoints/weight493084032.pth 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 322M/322M [00:24<00:00, 13.9MB/s] A group of men laying on top of a beach. electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % python3 predict.py --path ./image_files/test002.jpg --v v1 Using cache found in /Users/electron/.cache/torch/hub/saahiluppal_catr_master Downloading: "https://github.com/saahiluppal/catr/releases/download/0.1/weights_9348032.pth" to /Users/electron/.cache/torch/hub/checkpoints/weights_9348032.pth 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 322M/322M [00:23<00:00, 14.4MB/s] A man in a wetsuit laying on a beach. electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % python3 predict.py --path ./image_files/test002.jpg --v v4 Checking for checkpoint. Traceback (most recent call last): File "/Users/electron/Desktop/catr/predict.py", line 32, in <module> raise NotImplementedError('No model to chose from!') NotImplementedError: No model to chose from! electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % python3 predict.py --path ./image_files/test001.png --v v3 Using cache found in /Users/electron/.cache/torch/hub/saahiluppal_catr_master A man is riding a horse in a field electron@diynoMacBook-Pro catr % ls ./image_files test001.png test002.jpg test003.jpg electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % python3 predict.py --path ./image_files/test003.jpg --v v3 Using cache found in /Users/electron/.cache/torch/hub/saahiluppal_catr_master A group of young people posing for a picture. electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % electron@diynoMacBook-Pro catr % cat predict.py | import zsh: command not found: import electron@diynoMacBook-Pro catr % cat predict.py | grep import import torch from transformers import BertTokenizer from PIL import Image import argparse from models import caption from datasets import coco, utils from configuration import Config import os electron@diynoMacBook-Pro catr % ls LICENSE catr_demo.ipynb engine.py image_files predict.py README.md configuration.py finetune.py main.py requirements.txt __pycache__ datasets hubconf.py models electron@diynoMacBook-Pro catr %

はじめに 高校生という生き物はいつだって欲望に忠実である。 「星が見たい！」と言い出したのはいつだったか。 よーしわかった。プラネタリウムを作ろう。 目標 高校生が可能な範囲内でそれなりに良いものを作る。 作るものはプラネタリウムのドームと，恒星球と呼ばれる星の投影機。 準備するもの ドーム 段ボール ロール紙 ガムテープ 恒星球 PC（winでもmacでも，pythonとblenderが使えればよい） 厚紙 針 のり 豆電球 制作手順 ドーム 今回はJAXAのダンボールドームプロジェクトを参考にした。 やり方はすべて書いてあるので，ここでは実際に作ってみて感じたことを書いていこうと思う。 まず，組み立てる時にはそれなりの人数が必要になってくる点。 寸法通り作るとかなりでかいものとなる。 教室に入れるために土台となるA8は半分の450mmの高さとしたが，それでも5人くらいはいないと難しい。 また，JAXAの手順では紙を貼るときにスプレーのりを使用しているが，今回は予算的に厳しかったため，紙を折り返してガムテープで貼り付けた。 そうするとかなり安く抑えることができたのだが，ドームの見た目が少し残念なことになってしまった。 最後に，ダンボールの状態。 今回は家電量販店のごみをリユースしたのだが，そうすると折り目がついた部品が出てきてしまい，ドームの耐久性に少し不安が残る形となった。 余談だが，2020年に制作したときは，竹を使ってドームを制作した（方法はこちら）。 そのドームに白のカーテンを掛けてスクリーンとしたのだが，どうしても内部の竹が邪魔で投影できない星が出てきてしまった。 ただ，制作時間的には圧倒的に短く終えられるので，時間がない場合はこちらでも良いだろう。 恒星球 qiitaの記事なので，ここからが本題。 ネット上にはいろいろな人がプラネタリウムを制作している様子をアップロードしているが，恒星球は市販のものを使っているケースが多い。 だが，自分たちで作るからこそ価値があるというものだろう。 pythonの準備 最初にPCにpython3系をインストールする。 （分からない場合はこちら） ちゃんとpythonがPCで使えるようになっているか確認。 「win + r」でファイル名を指定して実行を立ち上げ，そこにcmdと打ち込むと黒い画面が出てくる。 そこに，pythonと入力 cmd C:\Users\user>python enter。 python3系が立ち上がり cmd >>> のようなものがでればOK。 出ない場合は恐らくPATHが通っていないので，環境変数を確認。 次にmatplotlibをインストールする。 さきほどの黒い画面に cmd >>>exit() と入力してpythonの対話モードを終了。 そのまま， cmd C:\Users\user>python -m pip install matplotlib と入力すれば，matplotlibがインストールできる。 次に，実際にpythonのプログラムを組むための環境を作っていく。 私のおすすめのエディタはatomなので，これを使うと想定しながら進める。 まず，こちらからatomをダウンロードする。 インストールが完了したら，atomを立ち上げ，左上のファイル，設定，インストールと進んでいく。 そうすると検索ボックスが出てくるので，そこに，atom-runnerと入力し，インストール。 これで，atom上でpythonを実行できるようになる。 ためしに新しいファイルtest.pyを作成し，以下のように入力してみる。 test.py print("Hello, World!!!!!") そのまま，「Alt + r」で実行してみると，下のほうに Hello, World!!!!! と表示されるだろう。 二次元座標系への星のプロット プログラムを組めるようにはなったが，星のデータがなければなにもできない。 今回はヒッパルコス星表に含まれる3215個分の星のデータを使用した。 ページに飛んでダウンロードを押すと，なんかよくわからない数字がたくさん出てくるとおもう。 これをすべてコピーして，atomに戻り，新しいファイルstar_plot.txtを作成してそこにペースト，保存する。 ここからは実際にプログラムを組んでいく。 まずは，先ほど保存した星の位置データを二次元座標系にプロットした画像データを作る。 実際に組んだものがこちら。 star_plot.py # -*- coding: utf-8 -*- import matplotlib.pyplot as plt import math import sys #入力は # 0 HIP番号 # 1 赤経時 # 2 赤経分 # 3 赤経秒 # 4 赤緯符号 # 5 赤緯度 # 6 赤緯分 # 7 赤緯秒 # 8 等級 #星の数 N = 3215 plot_data = [] with open('star_plot.txt', 'r', encoding='utf-8') as star: star_data = star.readlines() for i in range(N): star_data[i] = star_data[i].split(',') a = float(float(star_data[i][1])*15 + float(star_data[i][2])*0.25 + float(star_data[i][3])*0.0042) b = float(float(star_data[i][5]) + float(star_data[i][6])/60 + float(star_data[i][7])/3600) c = int(round(float(star_data[i][8]), 0)) if star_data[i][4] == "0": b = -b if c <= 5: plot_data.append([a,b,c,i]) star.close() x = [] y = [] blt = [] for i in range(len(plot_data)): x_ = round(float(plot_data[i][0]), 1) y_ = round(float(plot_data[i][1]), 1) blt_ = plot_data[i][2] x.append(x_) y.append(y_) blt.append(blt_) plt.figure(figsize=(36, 9), dpi=300) plt.xlim(0, 360) plt.ylim(-90, 90) plt.axis("on") for i in range(len(plot_data)): s_ = (6 - blt[i])**3 plt.scatter(x[i], y[i], s=s_, c=s_, cmap='spring', vmin=-50, vmax=343) plt.savefig("star_fig.png") 最初はこのプログラムをそのままコピペして使用してみて，いろいろいじくっていけばなんとなく分かってくると思う。 それから，私は別にプログラミングの専門家ではないので，汚い部分もあるがそこは許してほしい。 注意点として，星のデータを保存したstar_plot.txtとこのstar_plot.pyは同じディレクトリに入れておかないと動かない。ひとまず，デスクトップ等にplanetariumとかいうフォルダを作って，そこに２つとも入れておけばよい。 プログラムを「Alt + r」で実行すれば，以下のような画像ファイルができる。 画像の左上からsinカーブに沿って，星が密集している部分があると思うが，これが天の川である。 また，x座標70，y座標0の付近にオリオン座が見える。 気づいた人もいるかもしれないが，この画像は実際の星と比べて左右が反転している。 しかし，丸める際にこの印刷面を球の表にすれば，スクリーンに映るときは実際の星空の向きになるので問題ない。 正多面体上に星をプロット しかし，二次元座標上に星をプロットしただけでは，完成ではない。 例えば，メルカトル図法で描かれた地球を丸めても，球体の地球にはならない。 そこで今回は，正多面体の展開図にこの二次元座標で描かれた星をリプロットする方法をとった。 このリプロットに関してだが，しっかり計算すれば星が正しい位置に来る展開図を作成できると思う。 しかし，私はその方法をうまく実装することができなかった（うまく計算できた方はぜひやり方を教えていただきたい）。 私がとった方法は，3DCGソフトのBlenderを用いて正多面体をUV展開し，その展開図を作成した画像に合わせて，ペーパークラフト作成用のアドオンを用いて展開図を作るというものである（友人のKNにこのアイデアを頂き，Blenderの使い方まで教えてもらった。感謝）。 まずは，Blenderをここからインストールする。 Blenderを立ち上げ，新規ファイル，全般を選ぶと，以下のような画面になる。 この真ん中にある物体を右クリックして，削除。 「Shift + a」でメッシュ，ico球を選択すると，次はこんなものが出てくる。 これで正多面体が用意できた。 ここでおまじないとして，右上の「オプション」の下にある球のアイコンのうち，右から二番目を選択しておく（画像では右から三番目が選択された状態となっている）。 次に，画面を２分割し，どちらかの画面でUV展開を行う（今回は左の画面）。 左上にあるBlenderマークから少し左下にカーソルを動かすと，カーソルが十字型になる場所があるので，そこから右へドラッグすれば，２画面に分割することができる。 次に，Blenderマークのすぐ下にある井戸の井の字とボールのマークをクリックすると，ドロップダウンが開くため，そこからUVエディターを選択する。 そうするとこのような画面になる。 できたら，右の画面に戻り，先ほどの井戸とボールのマークの下にある「オブジェクト」をクリックし，編集に変更する。すると， この展開図ではいささか問題があるので，展開図を変更する。 右の正多面体にカーソルを合わせ，「a」を押し，その後，「u」を押すと，UV展開の方法がいくつか出てくる。 そこで円筒状投影を選択すれば， こうなる。 できたら，左の画面の「開く」というところから，pythonで作成したstar_fig.pngを開く。 左の画面の1つ1つの三角形が右の正多面体の面を表している。 このままではすべての面がうまく画像に合っていないので，点をうまいこと動かして画像の枠内にピッタリと収まるように変形していく。 Blenderを用いた方法で誤差が生じるのはココの操作である。 そのため，なるべく誤差を少なくするために，スケール変換で大まかな形は合わせてしまったほうが良い。 ショートカットキー a：全選択，s：スケール変換，s->x：横軸のみスケール変換，s->y：縦軸のみスケール変換，g：移動 ピッタリ合わせた後のものがこちら。 見るとよくわかると思うが，てっぺん付近や左右の端っこはかなりゆがんでしまっている。 簡単にやろうと思うとここら辺が限界。 次に，またBlenderマークの左下にカーソルを合わせて，十字型のカーソルを出した後，下にドラッグして３画面にする。 同様に，Blenderマークの下のアイコンをクリック，今度はシェーダーエディターを選択。 できたら，今出したウィンドウの右上あたりにある新規をクリックする。 次に，シェーダーエディターの追加ボタンを押して，「UVマップ」と「画像テクスチャ」の２つを追加する。 画面テクスチャの「開く」からstar_fig.pngを開く。 UVマップのUVの右にある青い点を，画像テクスチャのベクトルの左にある青い点にドラッグ＆ドロップ。 また，画像テクスチャの右のカラーからプリンシプルBSDFのベースカラーにドラッグ＆ドロップ。 ここまで完了すれば，一番最初のICO球を用意した画面に恒星球の完成形が映し出されているだろう。 正多面体の展開図を作る 恒星球はできたので，この展開図を作っていく。 参考にしたのはこちら この記事で解説されているので詳しいことは書かないが，アドオンが見つからない場合はgithubからインストールしなければならない場合があることと，シームの付け方を注意しないとペーパークラフトができないことがあるので注意。 完成したものがこちら。 ピンホール式プラネタリウム 完成したものを印刷して，光を通さない厚紙等にのりで貼り付ける。 組み立てる前に星の部分に針で穴を空けていく。 この時，星の等級によって穴の大きさを変えるのだが，思ったより大きく空けた方が良い。 その後，下の方の一面をはさみで切り取ってから，のりしろを組み合わせて正多面体をクラフトする。 完成したら正多面体の面がない部分から豆電球を入れて，光らせる。 うまくできていれば，こんな風に映るはず。 ポイントとして，豆電球はなるべくフィラメントの短い，明るいものが良い。 ただし，この２つの要素は相反する。 私たちの実験では，3.8Vの豆電球に定電圧電源で電流を流すと良い感じに投影できることが分かっている。 最後に ぶっちゃけ高校生が自力でここまでできるのかは疑問が残るところだが，部分的にでもプラネタリウムを作るときの参考にして頂ければ大変嬉しく思っている。 また，小中学生のお父さんお母さんが夏休みの自由研究を行う時の参考にしてもらってもいいだろう（誰の自由研究？）。 実際にドームの中で恒星球を光らせてみるとかなりの感動ものである。 こればっかりは体験してみないと分からないので，興味が湧いた方はぜひ作ってみてほしい。