はじめに 新型コロナの感染者数が増えてきて、東京、大阪、京都、兵庫で緊急事態宣言が発出されました。感染者数の情報は厚生労働省のページからダウンロードすることができます。まずは陽性者数のデータをダウンロードしてグラフ化してみます。 import pandas as pd import matplotlib.pyplot as plt import japanize_matplotlib import seaborn as sns def get_data(): url = 'https://www.mhlw.go.jp/content/pcr_positive_daily.csv' df = pd.read_csv(url, usecols=[0, 1], names=['date', 'positives'], skiprows=1, parse_dates=['date'], index_col='date' ) return df def main(): df = get_data() df.plot() plt.title('COVID19 日別感染者数') plt.show() if __name__ == '__main__': main() 素人の私が見ても、これは第四波が来てると思ってしまいます。 この先どうなるのでしょう？ 今回はProphetとARIMAを使用して、1ヶ月後までの感染者推移数を予測してみます。実はLSTMを使ったモデルでも試してみたのですが、うまく予想できなかったため別の方法で試してみることにしました。 Prophetで予測してみる ProphetはFacebookが開発した予測モデルです。 使い方に特長がありますが、モデルを作るのは簡単です。 入力データにはPandasのDataFrameで指定します。このとき時系列の列名をds、測定値の列名をyとしておく必要があります。更にcapという列を追加する必要があり、ここには取り得る最大値を指定します。今回は10000としておきます。データの取得とDataFrameの作成は下記のようになります。 url = 'https://www.mhlw.go.jp/content/pcr_positive_daily.csv' df = pd.read_csv(url, usecols=[0, 1], names=['ds', 'y'], skiprows=1, parse_dates=['ds'], ) df['cap'] = 10000 モデルの作成は下記のように行います。 model = Prophet() 何も指定しないと線形回帰になるので、今回は下記のようにlogistic回帰を指定します。 model = Prophet(growth='logistic') あとはデータを指定してフィッティングさせます。 model.fit(df) 将来の予想を行うには、将来の日付を指定したDataFrameを作成します。make_future_dataframeというメソッドを使用すると、指定した先までのデータフレームを作ってくれます。ここでもcap列を追加しておく必要があるので注意です。 df_future = model.make_future_dataframe(periods=30) df_future['cap'] = 10000 予測値を求めてみます。 predicts = model.predict(df_future) 下記のようにすると、結果をグラフで表示できます。 plt.xkcd()とすると、手書き風のグラフを作成できます。(すみません、やってみたかっただけです。) plt.xkcd() model.plot(predicts) plt.tight_layout() plt.title('Prediction COVID-19 by Prophet Model') plt.show() 増加傾向にあるのはわかります。1週間の周期性を持っているのもわかります。でも、実際の値にはフィットしていませんね... ARIMAで予測してみる 次はARIMAで予測してみます。感染者数は1週間の周期性があるのがわかっているので、今回はSARIMA(Seasonal ARIMA model)を使ってみます。 ARIMAモデルを理解しておらず説明がうまくできないので、ソースコードだけ書いておきます。 import datetime import pandas as pd from matplotlib import pylab as plt import japanize_matplotlib import statsmodels.api as sm def get_data(): url = 'https://www.mhlw.go.jp/content/pcr_positive_daily.csv' df = pd.read_csv(url, usecols=[0, 1], names=['date', 'positives'], skiprows=1, parse_dates=['date'], index_col='date' ) return df def main(): plt.xkcd() df = get_data() diff = df['positives'].diff() diff = diff.dropna() params = sm.tsa.arma_order_select_ic(diff, ic='aic', trend='nc') aic_order = params['aic_min_order'] ''' orderはarma_order_select_icで求めた値を指定。 seasonal_orderは1週間周期なので4番目に7を指定。 ''' model = sm.tsa.SARIMAX( df, order=(aic_order[0], 1, aic_order[1]), seasonal_order=(1, 1, 1, 7) ).fit() ''' 30日分の予想をしてみる ''' predict_period_from = df.index.max() predict_period_to = df.index.max() + datetime.timedelta(days=30) predict = model.predict(predict_period_from, predict_period_to) plt.plot(df, label='real') plt.plot(predict, label='predict') plt.title('Prediction COVID-19 by SARIMA Model') plt.savefig('covid19-arima2.png') plt.show() if __name__ == '__main__': main() なんか、いい感じで予測できていますね。 終わりに 今回はDeepLearningをあきらめて、別の方法でCOVID19の予測をしてみました。 ARIMAモデルのすごさもわかりました。DeepLearningと並行して勉強していきたいと思います。 ソースはGitHubに置いています。

QuoraにGoogleの ”Jin Yamanaka” 氏が投稿していた、 ラマヌジャンの公式を使った円周率の計算法を自分でも動かしてみた。 普通にmathライブラリをインポートしたpythonを使って計算すると、 20桁程度の精度しかない。しかもラマヌジャンの式は収束が速いようで、ほんの３ループほどで計算値の違いがなくなった。 1 : 3.141592 7300133056603 2 : 3.141592653589793 878 3 : 3.1415926535897932385 4 : 3.1415926535897932385 5 : 3.1415926535897932385 6 : 3.1415926535897932385 7 : 3.1415926535897932385 8 : 3.1415926535897932385 9 : 3.1415926535897932385 任意精度の結果が得られるmpmathライブラリを使って書き直してみた。 もちろん、まだ入れてない場合は追加でのインストールが必要だ。 import mpmath mpmath.mp.prec = 224 def get_one_over_pie(total_n): v = 0 for n in range(total_n): v += mpmath.factorial(4 * n) * (1103 + 26390 * n) / ((mpmath.factorial(n) ** 4) * (396 ** (4 * n) )) return mpmath.sqrt(8) / 9801 * v for n in range(1, 10): print(n, ": ", 1 / get_one_over_pie(n)) 1 : 3.141592 73001330566031399618902521551859958160711003355965653629013 2 : 3.141592653589793 87799890582630601309421664502932284887917396379151 3 : 3.14159265358979323846264 906570275889815667748046233478116839959564 4 : 3.1415926535897932384626433832795 5527315997421042037991121670389601 5 : 3.141592653589793238462643383279502884197 663818133030623976165591 6 : 3.14159265358979323846264338327950288419716939937 984683274351250728 7 : 3.14159265358979323846264338327950288419716939937510582 102093324247 8 : 3.141592653589793238462643383279502884197169399375105820974944592 76 9 : 3.14159265358979323846264338327950288419716939937510582097494459231 おそろしく収束が速い式だし、なんでこんな式が浮かんだのか想像もつかない。 天才って本当に居るんだな。

python のインストール手順 　MacにPythonの３系をインストールするのに、ハマったのでメモ phythonのインストール こちらのHomebrew + pyenvを使用した手順でインストールを行う Homebrewとpyenvのインストールは略 ターミナルで下記コマンドを実行して、 最初にインストール可能なヴァージョンの確認 $ pyenv install --list 実行結果 Available versions: 2.1.3 2.2.3 2.3.7 2.4.0 2.4.1 2.4.2 2.4.3 2.4.4 2.4.5 　・ 　・ 　・ 今回は3.9.4をインストールしたいので、下記を実行 $ pyenv install 3.9.４ 成功↓ python-build: use openssl@1.1 from homebrew python-build: use readline from homebrew Downloading Python-3.9.4.tar.xz... -> https://www.python.org/ftp/python/3.9.4/Python-3.9.4.tar.xz Installing Python-3.9.4... python-build: use readline from homebrew python-build: use zlib from xcode sdk Installed Python-3.9.4 to {path} 念の為、versionsの確認 $ pyenv versions ↓のようにインストールしたversionsが出ていればよい system * 3.9.4 (set by {path}) インストールしただけでは、適用されないので、 下記コマンドを実行 $ pyenv global 3.9.4 実行したら、versionの確認 $ python --version ↓mac標準のversionsのまま。なんでさ Python 2.7.16 https://github.com/pyenv/pyenv#homebrew-on-mac-os-x によると以下手順が必要とのこと $ pyenv init ↓実行結果の指示通り、.zshrcに"eval "$(pyenv init -)"を追記または、 記載していることを確認 # Load pyenv automatically by appending # the following to ~/.zshrc: eval "$(pyenv init -)" 下記、実行してversion確認 $ source ~/.bash_profile $ python -V ↓インストールしたものが適用されました。 Python 3.9.4 これで、phytonのインストールは完了

はじめに tensorflow2.xでKerasのModel.predictを繰り返し実行したら、メモリ使用量がもりもり増えていったので、メモ。 Inputサイズが大きいというのもありますが、128Gメモリ、バッチサイズ4000、4回のpreditほどで124G張り付き。。。 ソースまで深追いはしていませんが、現状メモリーリーク(memory leak）があるようです。 ubuntu==18.04 python==3.8.5 tensorflow==2.4.1 nvidia-driver==465.19.01 cuda==11.0 対策：うまく行かなかったパターン マニュアルにメモリをクリアしてみる。 疑似コードとしては以下のように修正 for x in dataset: y_pred = mode.predct(x) # 明示的に削除とGC del y_pred y_pred = None gc.collect() 結果としては、多少（数ループ）改善したものの8000回のループ回したい要件としてはNG 対策：うまく行ったパターン 参考にしたのは、↓にリンクされていたgithubのissueをよく読んでみました。 するとpredict_on_batchでは問題なく動作するとのこと 疑似コードとしては以下のように修正 for x in dataset: y_pred = mode.predct_on_batch(x) 20〜30Gぐらいで均衡を保つようになりました

はじめに 　本投稿では趣味で遊んでいる King's Raid のPT構成を最適化しようと組んだ実行コードを紹介する。趣味にプログラムを応用するのはモチベーションも高くできて、いい勉強になるしおススメ.+.(ﾉ･ω･)ﾉ*. 1. 準備 1.1 Import はじめのおまじないを唱えます(∩｡･ｏ･｡)っ.ﾟ☆｡'` import itertools import os import numpy as np from matplotlib import pyplot as plt import pandas as pd 1.2 関数定義 今回使う関数を定義しておく def list2dic(x): """ 数値のリストを辞書に変換（キャラとか攻撃上昇などのラベルを付ける） :param x: リスト :return: 辞書 """ dic = { 'キャラ': x[0], '攻撃上昇': x[1], '攻撃速度': x[2], '攻撃値': x[3], 'クリダメ': x[4], '防御貫通': x[5], '与ダメ': x[6], '被ダメ': x[7], '防御低下': x[8] } return dic def damage_rate(df, a=1.0): """ DataFrame（PT構成）からダメージ上昇率をざっくり計算する ※ 攻撃値加算がダメージに与える影響が不明だったので係数aで調整できるようにしておいた(´っ･ω･)っ :param df: DataFrame（PT構成） :param a: 攻撃値加算がダメージに与える影響 :return: ダメージ上昇率 """ f = lambda x: sum(x) # 上昇率を計算する return (1 + f(df['攻撃上昇']) + f(df['攻撃値']) * a) * (2 + f(df['クリダメ'])) * (1 + f(df['与ダメ'])) * (1 + f(df['被ダメ'])) def damage_rate_check(df, a=1.0): """ DataFrame（PT構成）から各上昇率をざっくり計算する ※ 攻撃値加算がダメージに与える影響が不明だったので係数aで調整できるようにしておいた(´っ･ω･)っ :param df: DataFrame（PT構成） :param a: 攻撃値加算がダメージに与える影響 :return: 各上昇率 """ f = lambda x: sum(x) # 上昇率を計算する data = np.array([(1 + f(df['攻撃上昇']) + f(df['攻撃値']) * a), (2 + f(df['クリダメ'])), (1 + f(df['与ダメ'])), (1 + f(df['被ダメ']))]) a = np.array([1/2, 1/4, 1]) plt.figure(figsize=(10, 4)) plt.bar(['ATK', 'CD', 'GSD', 'TD'], data, width=0.3) plt.bar(['ATK', 'CD', 'GSD'], data[:3]*a, width=0.3) plt.show() return dict(zip(['ATK', 'CD', 'GSD', 'TD'], data)) def pick(chara, df): """ 指定したキャラをピックしてDataFrame（PT構成）を作る :param chara: ピックするキャラ :param df: キャラ全体 :return: PT構成 """ return df[[c in chara for c in df['キャラ']]] def chara_combinations(chara, num): """ キャラ全体からnum人ピックしてできるPTの全パターンを作る :param chara: キャラ全体 :param num: ピック人数 :return: PTの全パターン """ return list(itertools.combinations(chara, num)) def concat_skill(df, skill): """ スキルを疑似的なPTメンバとして追加する :param df: PT構成 :param skill: スキル :return: PT構成＋スキル """ dic = list2dic(['サポスキル']+skill.tolist()) return pd.concat([df, pd.DataFrame([dic])], ignore_index=True) def analysis(df, base, relics, skill, num=6): """ 最適PTを分析する :param df: PT構成 :param base: 基礎バフ :param relics: DD遺物 :param skill: DDスキル :return: 全パターンのDataFrame, 全パターンの倍率, 倍率高い順のインデックス """ pt_df = [pick(pt, df) for pt in chara_combinations(df['キャラ'], num=num)] all_df = [pd.concat([pt, base, r, s]) for pt in pt_df for r in relics for s in skill] all_rate = np.array([damage_rate(df) for df in all_df]) sorted_index = all_rate.argsort()[::-1] return all_df, all_rate, sorted_index # スキル(atk:攻撃上昇, cd:クリダメ, gsd:与ダメ) atk = np.array([0.05, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00]) cd = np.array([0.00, 0.00, 0.00, 0.10, 0.00, 0.00, 0.00, 0.00]) gsd = np.array([0.00, 0.00, 0.00, 0.00, 0.00, 0.025, 0.00, 0.00]) def add_skill(df, skill_num): """ 不足分を考慮し最適なスキルを追加する :param df: スキルを付与する構成 :param skill_num: 追加するスキル数 :return: 最適なスキルを付与した構成 """ skill = np.array([0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00]) for i in range(skill_num): out = [damage_rate(concat_skill(df, skill+op)) for op in [atk, cd, gsd]] skill += [atk, cd, gsd][np.argmax(out)] return concat_skill(df, skill) # リストから辞書を生成 list2dict = lambda x: dict(zip([data[0] for data in x], range(len(x)))) 2. DDの定義 2.1 基礎バフ 常時のっているバフはほぼほぼこんな感じなのかな… # 0.キャラ（名前）, 1.攻撃上昇, 2.攻撃速度, 3.攻撃値, 4.クリダメ, 5.防御貫通, 6.与ダメ, 7.被ダメ, 8.防御低下 buff = [ ['オディ', 0.35, 0.00, 0.00, 1.10, 0.35, 0.00, 0.25, 0.00], # オディ君は固定バフ扱い∑(･`ω･ﾉ)ﾉ ['魔導セット', 0.25, 0.00, 0.00, 0.00, 0.00, 0.25, 0.00, 0.00], ['クラスバフ', 0.25, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00], ['DD超越', 0.65, 0.00, 0.00, 0.50, 0.20, 0.00, 0.00, 0.00] ] pt_relics = [ ['羽', 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.25, 0.00], ['鞭', 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.50, 0.00], ['王冠', 0.75, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00], ['優秀賞', 0.20, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00], ['かぼちゃ', 0.50, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00], ['ユノデザ', 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.25, 0.00], ['吐息', 0.00, 0.00, 0.00, 0.00, 0.00, 0.125, 0.00, 0.00], ] base = pd.DataFrame(map(list2dic, buff + pt_relics)) 2.2 魔導スキル DDの魔導スキルは攻撃特化(防御25攻撃50×4), 与ダメ特化(9525×4), バランスを準備(防御25攻撃50×2, 9525×2)。 # 0.キャラ（名前）, 1.攻撃上昇, 2.攻撃速度, 3.攻撃値, 4.クリダメ, 5.防御貫通, 6.与ダメ, 7.被ダメ, 8.防御低下 dd_skill_list = [ ['攻撃特化：スキル', 2.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00], ['与ダメ特化：スキル', 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 1.00, 0.00], ['バランス：スキル', 1.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.50, 0.00], ] dic_dd_skill = list2dict(dd_skill_list) dd_skill = [pd.DataFrame([list2dic(skill)]) for skill in dd_skill_list] 2.3 遺物 DDの遺物候補をいくつか用意 # 0.キャラ（名前）, 1.攻撃上昇, 2.攻撃速度, 3.攻撃値, 4.クリダメ, 5.防御貫通, 6.与ダメ, 7.被ダメ, 8.防御低下 dd_relics_list = [ ['大福', 1.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00], ['爪飾り', 0.00, 0.00, 0.00, 0.00, 0.00, 0.50, 0.00, 0.00], ['火の扉', 0.00, 0.00, 0.00, 0.00, 0.00, 0.25, 0.00, 0.00], ['獣人本', 0.00, 0.00, 0.00, 0.00, 0.00, 0.30, 0.00, 0.00], ['コイン', 0.50, 0.00, 0.00, 0.50, 0.00, 0.00, 0.00, 0.00], # ['シャク剣', 0.00, 0.00, 0.00, 0.50, 0.00, 0.50, 0.00, 0.00], # 使いにくいのでコメントアウト ] dic_dd_relics = list2dict(dd_relics_list) dd_relics = [pd.DataFrame([list2dic(relics)]) for relics in dd_relics_list] 3. PT検討 3.1. キャラリストの作成 キャラリストはCSVファイルから読み込めるようにした。（一応元のリストも残しておく） # chara = [ # ['ジェーン', 0.05, 0.00, 0.00, 0.00, 0.00, 0.15, 1.95, 0.50], # ['プリシラ', 0.40, 0.10, 0.33, 1.20, 0.00, 0.00, 0.00, 0.00], # [ '堕フレ', 0.40, 0.70, 0.07, 0.00, 0.50, 0.15, 1.10, 0.00], # ['ラブリル', 0.90, 0.00, 0.00, 2.75, 0.30, 0.00, 1.31, 0.25], # [ 'シア', 0.60, 0.40, 0.05, 1.00, 0.00, 0.20, 0.15, 0.20], # ['エステル', 0.25, 0.45, 0.00, 0.90, 0.15, 0.45, 0.00, 0.00], # ['ベロニカ', 0.15, 0.00, 0.00, 1.50, 0.00, 0.95, 0.75, 0.00], # ['アネット', 0.30, 0.20, 0.05, 0.60, 0.00, 0.00, 1.05, 0.00], # # ['オディ', 0.35, 0.00, 0.00, 1.10, 0.35, 0.00, 0.25, 0.00], # [ 'メイ', 0.90, 1.20, 0.00, 0.00, 0.96, 0.00, 0.25, 0.00], # ['メディ', 0.60, 0.00, 0.60, 0.50, 0.00, 0.00, 0.00, 0.00], # ['エスカー', 0.99, 0.00, 0.00, 0.00, 0.35, 0.04, 0.75, 0.00] # ] # chara = [ # ['シャクメ', 0.05, 0.00, 0.00, 0.00, 0.00, 0.20, 2.11, 0.00], # [ 'ロマン', 0.05, 0.00, 0.00, 0.00, 0.00, 0.00, 1.65, 0.00], # ['プリシラ', 0.40, 0.10, 0.33, 1.20, 0.00, 0.00, 0.00, 0.00], # [ '堕フレ', 0.40, 0.70, 0.07, 0.00, 0.50, 0.15, 1.10, 0.00], # ['ラブリル', 0.90, 0.00, 0.00, 2.75, 0.30, 0.00, 0.00, 0.00], # [ 'シア', 0.60, 0.40, 0.05, 1.00, 0.00, 0.20, 0.15, 0.20], # ['エステル', 0.25, 0.45, 0.00, 0.90, 0.15, 0.45, 0.70, 0.00], # ['ベロニカ', 0.15, 0.00, 0.00, 1.50, 0.00, 0.95, 0.75, 0.00], # # ['オディ', 0.35, 0.00, 0.00, 1.10, 0.35, 0.00, 0.00, 0.00], # [ 'メイ', 0.90, 1.20, 0.00, 0.00, 0.96, 0.00, 0.00, 0.00], # ['メディ', 0.60, 0.00, 0.70, 0.50, 0.00, 0.00, 0.45, 0.00], # ['エスカー', 0.99, 0.00, 0.00, 0.00, 0.35, 0.04, 0.00, 0.00], # ['ルシアス', 0.65, 0.00, 0.00, 0.00, 0.00, 0.20, 0.25, 0.00] # ] # df = pd.DataFrame(map(list2dic, chara)) df = pd.read_csv('MagicPT.csv') # df = pd.read_csv('PhysicalPT.csv') df 3.2. 最適PTの確認 最適なPTを取得する。best, second, third = [all_df[i] for i in sorted_index[:3]]　こんな感じで、上位3位まで取ってくるとかもできる。print(damage_rate(best))でどれだけ強いかスコアを表示する。この構成は2327.4らしい。また、棒グラフは 攻撃力（ATK）, クリダメ（CD）, 与ダメ（GSD）, 被ダメ（TD）を可視化したもの。 補足 　与ダメや攻撃、クリダメなどは盛りやすさで補正してどれだけ盛るのが最適かが変わります。例えば与ダメを+150%すると基礎値の100%+150%=250%で250%になります。サポの魔導スキルで与ダメは1箇所2.5%、クリダメは1箇所10.0%上がることを考慮すると、与ダメ+150%した場合はクリダメ+800%し、基礎値の200%+800%=1000%にする必要があります。わかりやすく表にして、最適（2.5倍×10倍）から与ダメスキルをクリダメスキルに変更すると下がる様子を表にしておきます。 　このことを考慮して、単純な攻撃やクリダメ、与ダメの倍率では何が足りないか判断が難しい…（棒グラフの青いバー）。そこで！この盛りやすさで補正をかけたのがオレンジのバーになります.+.(ﾉ･ω･)ﾉ.この高さが等しくなるとつよい！つまり、現状クリダメが足りていないことがわかる！ num = 6 # ピック人数 all_df, all_rate, sorted_index = analysis(df, base, dd_relics, dd_skill, num=num) best, second, third = [all_df[i] for i in sorted_index[:3]] print(damage_rate(best)) print(damage_rate_check(best)) best >>> 2327.4137737499996 3.3. スキルを割り振る 最高倍率目指して魔導のスキルを割り振るよ。とりあえずピックした全キャラ4枠使う想定… (add_skill(best, (num+1)*4)の(num+1)*4を書き換えれば追加するサポのスキル数を変更可能)。この構成は2995.8らしい。（さっきの棒グラフで確認した通り、クリダメが不足しているのでサポのスキルはクリダメ全盛になった） best_add_skill = add_skill(best, (num+1)*4) print(damage_rate(best_add_skill)) print(damage_rate_check(best_add_skill)) best_add_skill >>> 2995.79926775 3.4. オリジナル構成 このキャラは必須、このDD遺物やスキルは使えないなどのためにある程度操作できるようにした。試しにDDを与ダメ特化のスキルにして、大福持たせてみる。2595.3…よわい(´・ω・｀) my_pt = list(best['キャラ'][:num+1]) my_pt = pick(my_pt, df) my_dd_skill = dd_skill[dic_dd_skill['与ダメ特化：スキル']] my_dd_relics = dd_relics[dic_dd_relics['大福']] my_pt = pd.concat([my_pt, base, my_dd_skill, my_dd_relics]) my_pt_add_skill = add_skill(my_pt, (num+1)*4) print(damage_rate(my_pt_add_skill)) print(damage_rate_check(my_pt_add_skill)) my_pt_add_skill >>> 2595.30555025 4. おわりに ある程度綺麗に整理できたのでなんとなーく投稿してみた。たまにはこういう簡単なのもいいね(´っ･ω･)っ

はじめに こんにちは。 アメフトチーム、ノジマ相模原ライズの姫野です。 コーチの傍ら、アメフト界のデータ利用を推進するべく活動しています。 ↓前回記事です。 スポーツコーチがチームで使える検温アプリをDjangoで作ってみた 今回は、選手の全身写真から3Dモデリングを生成し、体脂肪率を算出しました。 結論から言うと実用レベルには至らない精度でしたがご紹介します。 課題 スポーツチームとなると選手の身体が資本なので、選手のコンディションチェックは重要なわけです。 体重が一番重要な要素の一つではありますが、その次に大事になるのは体脂肪率です。 ただ、体脂肪率を管理しようとしても、市販の数千~数万円の体組成計は精度が悪く、同じ日に測定しても5%とかは普通にズレます。 なにかいい方法がないかな～と探していたところに見つけたのが、コチラ！ 3Dボディスキャニング・システム『BODYGEE』の日本語版がリリース 『BODYGEE』最大の強みは、iPadと特殊な3Dカメラのみで、極めて精密なリアル3Dアバターを作成することができる、その”手軽さ”です。 「体組成」の重要な一要素である体脂肪率は、スキャンした3Dモデルとターンテーブルに内蔵された体重計から取得されるデータから算出されます。この3Dデータを活用した体脂肪率測定アプローチは、現在一般的なインピータンス方式を採用する体組成計に対して精度で優っていることが、大学の研究で明らかになっています。 月額数万円とお高くて手が出せませんが、手法を真似て3Dモデリングを自作していきます。(自分で作ってない) 3Dモデルの生成(PIFu) PIFu PIFuとは、FacebookとUSC(南カリフォルニア大)との共同研究で開発された、全身写真から3Dモデルを生成する技術です。(元論文) 最近PIFuHDという高解像度化されたものが新しく開発されたので、使っていきます。 ありがたいことに、写真を用意してコチラのGoogle Colabをポチポチするだけで3Dモデルができてしまいます。 やってみた 協力してくれたのは背番号4番の星野選手！ LBという守備の要のポジションで、声がやたらでかいチームのムードメーカーです。 できたモデルがこちら！ (元写真は半裸できたないエロいので載せてません) 足元がちょっと変ですが簡単に3Dモデルができました。 モデリングに使ったのは前からの写真ですが背中側もキレイにできています。 目がキュートですね。 ここで出来たobjファイルをBlenderで体積に換算していきます。 Blenderで体積に換算 Blenderとはフリーの3Dモデルソフトです。 体積を計算する機能があるので使っていきます。 こちらのサイトを参考にしました。 Blender上では91706cm3と算出されました。 ここでBlender上では縮尺が実世界とは違うので補正する必要があります。 blenderのメジャーで測ると星野選手は184.4cmだそうです。 ホントは星野選手は171cmなので、Blenderの縮尺は184.4/171となります。 縮尺の逆数の3乗をBlenderでの体積にかけて補正します。 体積計算.py blender_v = 91706 blender_h = 184.42 real_h = 171 calc_v = blender_v * (real_h/blender_h)**3 print('星野選手の体積: {} cm3'.format(round(calc_v))) 実行結果.py 星野選手の体積: 73108 cm3 体積に換算できたので、次にいよいよ体脂肪率の計算を行います。 体脂肪率の計算 体脂肪率の算出には、市販の体重計に付属されている測定法であるインピーダンス法のほかに、密度法があります。 現実ではプールに沈めて水中での体重を測ることで換算可能です。 密度から体脂肪率に変換するには以下の経験式が知られています。 fat = \frac{4.570}{D} - 4.124 \hspace{50pt} Brozek.et.al.\hspace{4pt}1963 この場合のDは密度(g/cm3)です。 星野選手の体重は87kgなので、こちらの式を使って体脂肪率を計算していきます。 結果 体脂肪率計算.py blender_v = 91706 blender_h = 184.42 real_h = 171 real_w = 87 calc_v = blender_v * (real_h/blender_h)**3 D = real_w / calc_v * 1000 fat = 4.570 / D - 4.124 print('星野選手の体脂肪率: {} %'.format(round(fat*100))) 実行結果.py 星野選手の体脂肪率: -28 % ん？ 人間の限界超えました。 星野選手の反応 気に入ってくれました。 なぜうまくいかなかったか 今回やったことは、 ①写真から3Dモデリング生成 ②Blenderで実体積に換算 ③密度法で体脂肪率に換算 です。 どこかのプロセスに原因がありそうですが、①では見た目でそれほど違和感ないモデルできているし(足元は汚いけどマイナスになる程か？)、③に至っては医科学会の権威？のような経験則ですので文句のつけようがありません。 ②が怪しそうです。 髪の毛を除いて再計算 Blenderの身長をもう一度測ります。 髪の毛がもっこりしてるので目視で目減らすと、182.16cmとなりました。 もう一度計算してみます。 体脂肪率計算.py blender_v = 91706 blender_h = 182.16 real_h = 171 real_w = 87 calc_v = blender_v * (real_h/blender_h)**3 D = real_w / calc_v * 1000 fat = 4.570 / D - 4.124 print('星野選手の体脂肪率: {} %'.format(round(fat*100))) 実行結果.py 星野選手の体脂肪率: -14 % なんと-28%から-14%まで変化しました。 依然としてマイナスですが、身長が大きく効いているのは間違いなさそうです。 縮尺の3乗で掛け算してくるのでちょっとのズレが影響大きいんですね。 縮尺のずれが精度不足になっていたようです。 髪の毛の考慮でも大きく変わりますが、足元はうまく撮れてませんし、写真の画角などでも大きく変わるので原理的には難しそうです。。。 今後 選手にはキャップをかぶせて定点カメラで写真を撮り、絶対尺度となるルーラーなどを用意しようと思います。 また、ノジマ相模原ライズでは一緒に活動してくれるエンジニアを募集中です！ スタッフ申し込みはこちらのフォームからお願いします。

はじめに 最近技術調査でHuggingFaceのサンプルコードを動かしているのですが、なかなか動かないことが多く(公式のものでも)、今回は動くようにするところまでできたので、まとめたいと思います。 問題 こちらの記事を参考にハンズオンをGoogle Colabratory`で行っていたところ、 import tensorflow as tf from transformers import TFGPT2LMHeadModel, GPT2Tokenizer tokenizer = GPT2Tokenizer.from_pretrained("gpt2") # 警告を避けるために、EOS トークンを PAD トークンとして追加 model = TFGPT2LMHeadModel.from_pretrained("gpt2", pad_token_id=tokenizer.eos_token_id) この最初の箇所で以下のようなエラーが発生しました。 --------------------------------------------------------------------------- ImportError Traceback (most recent call last) <ipython-input-23-4d671ac52c7c> in <module>() 1 import tensorflow as tf ----> 2 from transformers import TFGPT2LMHeadModel, GPT2Tokenizer 3 4 5 tokenizer = GPT2Tokenizer.from_pretrained("gpt2") ImportError: cannot import name 'TFGPT2LMHeadModel' from 'transformers' (unknown location) また、いくつかのサイトをみたところtensorflow==2.0.1を入れると使えるよというissueがあったが動きませんでした。 原因 同様の報告がバグとして報告され、installationの記事が追加されたようです。 https://github.com/huggingface/transformers/issues/3396 This repository is tested on Python 3.6+, PyTorch 1.0.0+ (PyTorch 1.3.1+ for examples) and TensorFlow 2.0. You should install ? Transformers in a virtual environment. If you're unfamiliar with Python virtual environments, check out the user guide. 他のライブラリが影響しているようでColabratoryのようなライブラリがたくさんある環境では、ほかのライブラリの影響のせいでimportできないらしいです。 詳しい原因はまだわかっていないが、自分で作成した環境で実行したところ記事通りのコードで動作した。 おわりに FuggingFaceを理解するためにとりあえずColabで動かすという流れをしていますが、Colabで動かないことがあることがよくわかりました。引き続きいろいろ試してHuggingFaceを理解したいと思います。 参考サイト https://note.com/npaka/n/n5d296d8ae26d https://github.com/huggingface/transformers/issues/3396

pythonを使用してExcelファイルの操作を勉強しています。 本日の気づき(復習)は、日付の代入や入力内容の削除に関してです。 pythonでExcelを操作するため、openpyxlというパッケージを使用しています。 発注書とかの書類をエクセルで作ってある場合って多いと思います。 新しいブックを作成するたびに前回のブックをコピーして内容を消して・・・。 そこまでの手間ではありませんが、プログラムで行った方がミスも減るのでお勧めかなとは思います。 と言う訳で、上の様なブック「申請書」を この様に初期化して、日付だけ当日のモノを入力しておきたいです。 date.todayメソッド cell.value = date.today() 日付の入力にはdatetimeモジュールを使います。 時刻も記入したい場合は cell.value = datetime.now() 上記のようにdatetime.nowメソッドを使用します。 入力内容削除 cell.value = None Noneを代入することで、内容を削除します セルを削除する訳では無いのでそこは安心です。 最終的なコード from datetime import date, datetime from openpyxl import load_workbook wb = load_workbook('申請書.xlsx') ws = wb.active ws['D4'].value = '営業部' ws['D5'].value = '田中一郎' ws['D6'].value = date.today() for row in ws.iter_rows(min_row=9, max_row=ws.max_row-6, min_col=2, max_col=7): row[0].value = None row[4].value = None row[5].value = None wb.save('申請書_初期化.xlsx') セルの記述を「=TODAY()」にしておいてもいいのですが、 そちらだと、日付をまたいだ時がめんどそうですね。

MacBook(m1)のpython環境構築 WindowsPCが壊れて新しくMacBookを購入 調べて環境構築したときの自分用メモです 環境 ・PC：13インチMacBook Pro (Apple M1チップ) ・OSバージョン：11.2.3 ・ターミナル：iTerm2 ・ARMアーキテクチャ 情報だと、homebrewを使用するためにはRosettaに切り替えないといけないとのことだったので(m1買ったのに...)condaで構築しました 構築にはこちらのサイトを参考にしました condaの導入 1. 「Rosetta上なのか」「ARMアーキテクチャ上なのか」を確認する ・ターミナル上でuname -mと入力 ・arm64と出力 : ARMアーキテクチャ　(こちらで構築) ・x86_64と出力 : Rosetta or Intelアーキテクチャ 2. miniforgeインストール arm64に対応しているminiforgeをインストールしてconda環境を構築していきます(理由は参考サイト参照) github から「Miniforge3-MacOSX-arm64」をダウンロード ターミナルでbash Miniforge3-MacOSX-arm64.sh と打つ すべて「yes」で答える ターミナルのシェルがbashからzshになったらOK 3. condaの初期化処理 コマンド source ~/.zshrc で .zshrc を読み込む ・するとcondaのbase環境が有効化される ・コマンドラインの頭に(base)がつく ターミナルを立ち上げるたびにcondaのbase環境が起動しないようにする(任意) ・conda config --set auto_activate_base false conda環境操作メモ 1.起動〜停止 ・conda activate : conda環境(base)起動 ・conda create -n [name] python=3.9 : 環境を作る(pythonのバージョンは各自合ったものを) ・conda activate [name]　: [name]というconda環境を起動する ・conda deactivate : 起動しているconda環境停止 2.その他(パッケージ系) conda info -e : 起動しているconda環境の情報をみる conda list : インストール済みのパッケージの確認 conda install [pakcage name] : パッケージのインストール jupyter notebook : jupyterNotebook 起動 ※jupyter終了時には必ずブラウザで「終了」を押す(qが使えなかった...) ※あくまで自分用メモなのであまり詳しく書いてません。詳しく知りたい方は参考サイトか他でググってください [参考]M1 MacにPythonインストールして開発環境構築してみた[https://zenn.dev/osuzuki/articles/380be0f682d72d]

1.ボタン名でクリック CssSelectorを使用すると、forループ等をせずにエレメントを見つけられます。 HTML.html <input type="submit" value="登録" /> C#.cs driver.FindElement(By.CssSelector("input[value='登録']")).Click(); または driver.FindElement(By.CssSelector("input[value='登録']")).Submit(); Java.java driver.findElement(By.cssSelector("input[value='登録']")).click(); または driver.findElement(By.cssSelector("input[value='登録']")).submit(); Python.py driver.find_element_by_css_selector("input[value='登録']").click() または driver.find_element_by_css_selector("input[value='登録']").submit()

1.ボタン名(ボタンの表示文字列)でクリック CssSelectorを使用すると、forループ等をせずに表示文字列等でエレメントを見つけられます。 HTML.html <input type="submit" value="登録" /> C#.cs driver.FindElement(By.CssSelector("input[value='登録']")).Click(); または driver.FindElement(By.CssSelector("input[value='登録']")).Submit(); Java.java driver.findElement(By.cssSelector("input[value='登録']")).click(); または driver.findElement(By.cssSelector("input[value='登録']")).submit(); Python.py driver.find_element_by_css_selector("input[value='登録']").click() または driver.find_element_by_css_selector("input[value='登録']").submit() 記述方法は異なりますが、XPathで同じことができます。 C#.cs driver.FindElement(By.XPath("//input[@value='登録']")).Click(); または driver.FindElement(By.XPath("//input[@value='登録']")).Submit(); Java.java driver.findElement(By.xpath("//input[@value='登録']")).click(); または driver.findElement(By.xpath("//input[@value='登録']")).submit(); Python.py driver.find_element_by_xpath("//input[@value='登録']").click() または driver.find_element_by_xpath("//input[@value='登録']").submit() XPathの参考：Seleniumで要素を選択する方法まとめ - Qiita

はじめに Pythonで正規表現による文字データの検索方法について調べたので、備忘として記録しています。 パターンマッチングのためのPython関数 Pythonにはパターンマッチングを使用するための標準ライブラリとしてreモジュールが用意されています。 match()関数 第一引数に検索するパターンを指定し、第二引数にソース文字列を指定します。match()関数はオブジェクトを返します。パターンマッチングした場合は、group()関数で確認することができます。マッチングしなかった場合は、Noneを返します。 sample01 >>> import re >>> source = 'This table lists the character \ reference names that are supported by HTML, \ and the code points to which they refer. \ It is referenced by the previous sections.' >>> result = re.match(r'This', source) >>> result.group() 'This' >>> result = re.match(r'That', source) >>> result is None True >>> search()関数 match()関数は、ソース文字列の先頭を検索するものでしたが、ソース文字列の途中にある文字列とマッチングをすることができます。 sample02 >>> import re >>> source = 'This table lists the character \ reference names that are supported by HTML, \ and the code points to which they refer. \ It is referenced by the previous sections.' >>> result = re.search(r'lists', source) >>> result.group() 'lists' >>> findall()関数 search()関数は、先頭からマッチングして１つマッチしたらそこで処理が終わりでしたが、findall()関数では末尾までマッチングを行います。 sample03 >>> import re >>> source = 'This table lists the character \ reference names that are supported by HTML, \ and the code points to which they refer. \ It is referenced by the previous sections.' >>> result = re.findall(r'th', source) >>> result ['th', 'th', 'th', 'th', 'th'] >>> マッチングの高速化 事前にパターンをコンパイルしておくことで、マッチング処理の高速化を図ることができます。compile()関数を使います。 sample04 >>> import re >>> source = 'This table lists the character \ reference names that are supported by HTML, \ and the code points to which they refer. \ It is referenced by the previous sections.' >>> ptn = re.compile('th') >>> result = ptn.findall(source) >>> result ['th', 'th', 'th', 'th', 'th'] >>> 正規表現の例 Pythonで利用できる正規表現をリストにしてみました。 文字 解説 例 マッチングされる文字列 . 任意の１文字 x.z xyz, xaz, xxz * 0回以上の繰り返し xz* x, xz, xzz, xzzz + 1回以上の繰り返し xz+ xz, xzz, xzzz ? 0回または1回 xz? x, xz ^ 文字列の先頭 ^xyz xyz, xyzabc $ 文字列の末尾 xyz$ abcxyz, xyz [xyz] x or y or z [xyz] x, y, z [^xyz] not (x or y or z) [^xyz] a, b, c {m} m回連続 x{3} xxx {m,n} m〜n回の連続 x{3,5} xxx, xxxx, xxxxx (expr) 文字列 (abc){2} abcabc x｜y x or y x｜y x, y 特殊文字 文字 説明 別の書き方 \d 任意の数字１文字 [0-9] \D 任意の数字１文字以外 [^0-9] \s 任意の空白１文字 [\t\n\r\f\v] \S 任意の空白１文字以外 [^\t\n\r\f\v] \w 任意の英数字 [a-zA-Z0-9_] \W 任意の英数字以外 [^a-zA-Z0-9_] sample05 >>> import string >>> printable = string.printable >>> printable '0123456789abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ!"#$%&\'()*+,-./:;<=>?@[\\]^_`{|}~ \t\n\r\x0b\x0c' >>> len(printable) 100 >>> import re >>> re.findall(r'\d', printable) ['0', '1', '2', '3', '4', '5', '6', '7', '8', '9'] >>> re.findall(r'\s', printable) [' ', '\t', '\n', '\r', '\x0b', '\x0c']

はじめに Pythonでアプリを作りました（以下main.py）。このmain.pyでは、他のデータファイルの参照をしています。普通はカレントディレクトリからの相対パス、main.pyからの相対パスなどで管理すると思います。 # カレントディレクトリからの相対パス PATH = 'file.txt' # main.pyからの相対パス PATH = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'file.txt') # python3.8以前 PATH = os.path.join(os.path.dirname(__file__), 'file.txt') # python3.9以降 PATH = Path(__file__).parent // 'file.txt' # pathlibを使用した場合 さて、このmain.pyを実行ファイル化しました（以下main.exe）。しかしmain.exeを実行すると、バグが生じます。どうやらパス周りがバグっていてファイルがないといわれました。今回はこのバグの原因を調べたところ、どうやら__file__に格納してあるパスが原因のようです。 解決方法 => 検証結果まとめ 検証 実行ファイルにおける__file__, os.getcwd(), sys.executableの挙動を確認しました。 検証内容 以下の内容を検証します。 __file__: main.pyの相対パス(python < 3.9)、絶対パス(python >= 3.9) os.getcwd(): カレントディレクトリの絶対パス os.path.abspath(__file__): main.pyの絶対パス Path(__file__).parent: main.pyがあるディレクトリ sys.executable: Python インタプリタの実行ファイル(python.exe)の絶対パスを示す 検証環境 Windows10 Python 3.9.4 Poetry version 1.1.6 pyproject.toml python = "^3.9,<3.10" pyinstaller = "^4.3" cx-Freeze = "^6.6" py2exe = "^0.10.4" pythonのバージョンは、デフォルト指定で"^3.9"で、これは3.8<=pythonのバージョン<4.0を示しています。しかしpy2exeが"^3.9,<3.10"という指定のため上限も指定しています。 その他 デフォルトのカレントディレクトリC:\test ターゲットファイルC:\test\test.py 実行ファイル名: test.exe としています。 検証結果 pyinstaller, cx_Freeze, py2txtそれぞれでexe作成を行い、実行ファイルを実行します。 実行の仕方として、 C:\testからコンソール実行 test.exeのあるディレクトリに移動してコンソール実行 エクスプローラーから本体をダブルクリックで実行 デスクトップ作成したショートカットをダブルクリックで実行 の3種類試しています。 ソースコードは test.py import os import sys from pathlib import Path print(f"1. __file__ => {__file__}") print(f"2. os.getcwd() => {os.getcwd()}") print(f"3. os.path.abspath(__file__) => {os.path.abspath(__file__)}") print(f"4. Path(__file__).parent => {Path(__file__).parent}") print(f"5. sys.executable => {sys.executable}") 0. pythonでそのまま実行 コントロールとして、exe化せずにそのまま実行しました。 # poetry run pytnon test.py で実行 1. __file__ => C:\test\test.py 2. os.getcwd() => C:\test 3. os.path.abspath(__file__) => C:\test\test.py 4. Path(__file__).parent => C:\test 5. sys.executable => C:\test\.venv\Scripts\python.exe 1~5で不審な点は見当たりません。poetryの仮想環境で検証しているため5.のパスは仮想環境内のpythonのパスを示しています。 1. pyinstallerの場合 以下のコマンドでexeの作成を行いました。 poetry run pyinstaller test.py C:\test\dist\test\test.exeが実行ファイルのパスです。 結果 # dist\test\test.exe での実行 1. __file__ => C:\test\dist\test\test.py 2. os.getcwd() => C:\test 3. os.path.abspath(__file__) => C:\test\dist\test\test.py 4. Path(__file__).parent => C:\test\dist\test 5. sys.executable => C:\test\dist\test\test.exe # cd dist\test && test.exe での実行 1. __file__ => C:\test\dist\test\test.py 2. os.getcwd() => C:\test\dist\test 3. os.path.abspath(__file__) => C:\test\dist\test\test.py 4. Path(__file__).parent => C:\test\dist\test 5. sys.executable => C:\test\dist\test\test.exe # 本体をダブルクリックで実行 1. __file__ => C:\test\dist\test\test.py 2. os.getcwd() => C:\test\dist\test 3. os.path.abspath(__file__) => C:\test\dist\test\test.py 4. Path(__file__).parent => C:\test\dist\test 5. sys.executable => C:\test\dist\test\test.exe # ショートカットをダブルクリックで実行 __file__ => C:\test\dist\test\test.py os.getcwd() => C:\test\dist\test os.path.abspath(__file__) => C:\test\dist\test\test.py Path(__file__).parent => C:\test\dist\test sys.executable => C:\test\dist\test\test.exe の結果は変です。4通りともC:\test\dist\test\test.pyというパスを表示しましたが、C:\test\dist\testにあるのはtest.exeでありtest.pyは存在しません。 の結果で不審な点はありません。ダブルクリックの際は本体があるディレクトリを示しています。 の結果で不審な点はありません（python3.9で実行しているので当たり前です）。 の結果で不審な点はありません（1.のファイル名はおかしいですが親ディレクトリまでのパスはtest.exeのパスと等しいので使えるっちゃあ使える？） これはびっくり。python.exeでもNoneでもなく、作成した実行ファイルのパスが格納されています。このパスをもとに相対パスなどを指定してやるのが安全で便利そうです。 2. cx_Freezeの場合 cx_Freezeでexeを作成する場合、setup.pyという設定ファイルを記述する必要があります。今回はpy2txtのsetup.pyと差別化するため、setup_cx_freeze.pyという名前にしました。 ソースコードは以下のようになります。 setup_cx_freeze.py import sys from cx_Freeze import setup, Executable build_exe_options = {"packages": ["os", "sys", "pathlib"], "excludes": ["tkinter"]} setup( name = "test", version = "0.1", description = "My GUI application!", options = {"build_exe": build_exe_options}, executables = [Executable("test.py")] ) 見よう見まねで作成しました。 exeの作成は以下のコマンドにより行いました。 poetry run python setup_cx_freeze.py build C:\test\build\exe.win-amd64-3.9\test.exeが実行ファイルのパスです。 結果 # build\exe.win-amd64-3.9\test.exe での実行 1. __file__ => test.py 2. os.getcwd() => C:\test 3. os.path.abspath(__file__) => C:\test\test.py 4. Path(__file__).parent => . 5. sys.executable => C:\test\build\exe.win-amd64-3.9\test.exe # cd build\exe.win-amd64-3.9 && test.exe での実行 1. __file__ => test.py 2. os.getcwd() => C:\test\build\exe.win-amd64-3.9 3. os.path.abspath(__file__) => C:\test\build\exe.win-amd64-3.9\test.py 4. Path(__file__).parent => . 5. sys.executable => C:\test\build\exe.win-amd64-3.9\test.exe # 本体をダブルクリックで実行 1. __file__ => test.py 2. os.getcwd() => C:\test\build\exe.win-amd64-3.9 3. os.path.abspath(__file__) => C:\test\build\exe.win-amd64-3.9\test.py 4. Path(__file__).parent => . 5. sys.executable => C:\test\build\exe.win-amd64-3.9\test.exe # ショートカットをダブルクリックで実行 1. __file__ => test.py 2. os.getcwd() => C:\test\build\exe.win-amd64-3.9 3. os.path.abspath(__file__) => C:\test\build\exe.win-amd64-3.9\test.py 4. Path(__file__).parent => . 5. sys.executable => C:\test\build\exe.win-amd64-3.9\test.exe の結果は変です。ファイル名でなく相対パスとして考えると、カレントディレクトリ中のtest.pyを指していますが、実際にそのようなファイルはありません。 の結果で不審な点はありません。ダブルクリックの際は本体があるディレクトリを示しています。 の結果は変です。os.path.join(os.getcwd(), __file__)という雰囲気です。 1.の結果の親を表すため、カレントディレクトリの相対パス.が表示されています。python.exeからの相対パスを知りたい場合にこのように書くのは危ないかな。 pyinstallerと同様に、作成した実行ファイルのパスが格納されています。便利。 3. py2exeの場合 py2exeでexeを作成する場合、setup.pyという設定ファイルを記述する必要があります。今回はcx_Freezeのsetup.pyと差別化するため、setup_py2exe.pyという名前にしました。 ソースコードは以下のようになります。 setup_py2exe.py from distutils.core import setup import setuptools import py2exe setup( console = [{'script': 'test.py'}] # exe化するファイル ) 見よう見まねで作成しました。 exeの作成は以下のコマンドにより行いました。 poetry run python setup_py2txt.py py2txt C:\test\dist\test.exeが実行ファイルのパスです。 結果 # dist\test.exe での実行 1. __file__ => # 2. os.getcwd() => C:\test 3. os.path.abspath(__file__) => # 4. Path(__file__).parent => # 5. sys.executable => C:\test\dist\test.exe # cd dist && test.exe での実行 1. __file__ => # 2. os.getcwd() => C:\test\dist 3. os.path.abspath(__file__) => # 4. Path(__file__).parent => # 5. sys.executable => C:\test\dist\test.exe # 本体をダブルクリックで実行 1. __file__ => # 2. os.getcwd() => C:\test\dist 3. os.path.abspath(__file__) => # 4. Path(__file__).parent => # 5. sys.executable => C:\test\dist\test.exe # ショートカットをダブルクリックで実行 1. __file__ => # 2. os.getcwd() => C:\test\dist 3. os.path.abspath(__file__) => # 4. Path(__file__).parent => # 5. sys.executable => C:\test\dist\test.exe __file__は未定義でエラーとなりました。そのため__file__を用いる1. 3. 4. をコメントアウトしました。 # の結果で不審な点はありません。ダブルクリックの際は本体があるディレクトリを示しています。 # # pyinstaller, cx_Freezeと同様に、作成した実行ファイルのパスが格納されています。便利。 まとめ 今回はpyinstaller, cx_Freeze, py2exe3種類の実行ファイル生成方法における、__file__, os.getcwd(), sys.executableの挙動を確認しました。 分かったこと __file__は実行ファイルの生成方法によって指し示す先が異なったり、相対・絶対パスが変わったりするので使わない方が良い。 os.getcwd()は実行ファイルの生成方法に依らず、正しい挙動を示した。 ダブルクリックの際の挙動はショートカット、本体で全く同じ挙動だった。 sys.executableは、実行ファイルの場合は自分自身の絶対パスを返すという挙動だった。 ということで、カレントディレクトリからの相対パスでパスを表現したい場合は今まで通りで問題なさそうだが、自分自身の場所からの相対パスで表したい場合はsys.executableを使う必要がある。しかし、.pyのままだとsys.executableはpython.exeの場所を示してしまう。どうしよう。 実は cx_freezeの公式ドキュメントに書いてありました。 Applications often need data files besides the code, such as icons. Using a setup script, you can list data files or directories in the include_files option to build_exe. They’ll be copied to the build directory alongside the executable. Then to find them, use code like this: def find_data_file(filename): if getattr(sys, "frozen", False): # The application is frozen datadir = os.path.dirname(sys.executable) else: # The application is not frozen # Change this bit to match where you store your data files: datadir = os.path.dirname(__file__) return os.path.join(datadir, filename) An alternative is to embed data in code, for example by using Qt’s resource system. サンプルコードでは、任意のファイルにアクセスしたい場合に、実行ファイルからの相対パスを引数として、絶対パスを返してくれる関数です。 実行ファイル化された場合はif文のコードが実行され、Pythonのまま実行した場合はelse以降が実行されます。 cx_freezeの公式ドキュメントの内容ですが、pyinstaller, py2exeを用いた場合も同様の対処で解決可能だと考えられます。

はじめに 住んでいる地域について天気予報を得るために、無料でAPIを提供している OpenWeatherMapを使ってみる 概要 OpenWeatherMapはGUIで天気予報を取得することができるのですが、今回はAPIを使って天気や気温を取ってきます。 初期設定 いくつかプランがあるので、無料のプランに登録し、APIkeyを手に入れます。無料枠だと5日後までの天気予報を取得出来て、1分間に60回までのAPIコール制限があるみたいです。 天気予報を取得する そもそもAPIを試したいならcurlを使えば簡単に試すことはできるのですが、今後システムに組み込んでいく予定なのでPythonで実装します。 エンドポイントはapi.openweathermap.org/data/2.5/weather パラメータは下記の通り パラメーター 説明 q 都市名 appid APIkey mode 応答方式 units 測定単位 lang 出力言語 都市名を指定してAPIリクエストを送る apicall.py import json import requests #パラメーター params={"q":"Kobe","appid":apikey} url="http://api.openweathermap.org/data/2.5/forecast" res=requests.get(url,params=params) k=res.json() jsonText = json.dumps(k["list"][0], indent=4) print(jsonText) Kobe(神戸)で試してみた結果 { "dt": 1619406000, "main": { "temp": 289.15, "feels_like": 287.75, "temp_min": 288.84, "temp_max": 289.15, "pressure": 1020, "sea_level": 1020, "grnd_level": 1016, "humidity": 36, "temp_kf": 0.31 }, "weather": [ { "id": 801, "main": "Clouds", "description": "few clouds", "icon": "02d" } ], "clouds": { "all": 20 }, "wind": { "speed": 8.04, "deg": 3, "gust": 9.57 }, "visibility": 10000, "pop": 0, "sys": { "pod": "d" }, "dt_txt": "2021-04-26 03:00:00" } 天気予報は取得できたのですが14時に取得したデータが"dt_txt": "2021-04-26 03:00:00"なのでリアルタイムの情報ではなく日単位での取得になるっぽいです。まあ無料枠なので当然と言えば当然。一番近い時間帯の情報を取得することで解決。 天気予報を表示 情報自体は取得できたので気温、湿度、天気、風速、取得日時を抜き出して天気予報っぽく表示します。 #省略 jsonText=res.json() print("都市：Kobe",end="\n\n") for i in range(0,5): print("気温：",jsonText["list"][i]["main"]["temp"])#気温 print("湿度：",jsonText["list"][i]["main"]["humidity"])#湿度 print("天気：",jsonText["list"][i]["weather"][0]["main"])#天気 print("風速：",jsonText["list"][i]["wind"]["speed"])#風速 print("取得日時：",jsonText["list"][i]["dt_txt"],end="\n\n")#取得日時 結果 気温： 289.14 湿度： 27 天気： Clear 風速： 7.34 取得日時： 2021-04-26 06:00:00 気温： 286.29 湿度： 39 天気： Clouds 風速： 6.69 取得日時： 2021-04-26 09:00:00 気温： 283.53 湿度： 50 天気： Clouds 風速： 4.59 取得日時： 2021-04-26 12:00:00 気温： 282.36 湿度： 58 天気： Clouds 風速： 3.12 取得日時： 2021-04-26 15:00:00 気温： 281.99 湿度： 60 天気： Clouds 風速： 2.57 取得日時： 2021-04-26 18:00:00 温度表示がケルビンなのと最低間隔が3時間という制限はありますが、工夫次第では十分に使えそうです。表示設定を日本語にしたり、ほかにも取得できる情報はたくさんあるので、詳しい内容はドキュメントを参照してください。 さいごに OpenWeatherMapを使うことで、比較的簡単に情報を取ってくることができました。今後はこれで夕立を通知してくれるBOTを開発したいです。

はじめに 久しぶりにtensorflowをインストールしたのですが、なぜかインポートできず(tensorflowがないよと言われる)、解決方法がすぐに見つからなかったため記事にしました。普段意識していないところが原因でしたので、まとめておきます。 環境 python3.8 Docker jupyter notebook 問題 !pip install tensorflow import tensorflow as tf を実行すると、以下のよエラーが発生する この記事を確認したところ--userをつけることで権限をつけることでインストール先を~./localに変更できる。 https://boook24.com/?p=518 !pip install tensorflow --user import tensorflow as tf を実行するとなぜかno moduleエラーとなる。インストールしているのになぜとなった。 解決方法 !pip istallはPython2にインストールするが、python3のライブラリを読み込んでいるので、!pip3 installとする必要があった。 !pip3 install tensorflow import tensorflow as tf で無事読み込むことができました。 終わりに pipとcondaの共存やpipとpip3の違いなど、どこにエラーの原因が潜んでいるかなかなか難しいです。

前回の『量子化学入門　Roothaan SCF計算（続報）』の続報2です。 SCF計算の非経験的な解法(ab initio 法)でプログラムに組み込むために、解法ページのMyLibraryのガウス型への変更を、少し試行してみました。 ガウス型関数を具体的に次のように置きます。 φ_{A1}φ_{A1} = φ_{P1} = cαe^{-αr_{A1}^2} \\ φ_{B2}φ_{B2} = φ_{Q2} = c'α'e^{-α'r_{B2}^2} \iintφ_{A1}φ_{A1}\frac{1}{r_{12}}φ_{B2}φ_{B2} \, dτ_1dτ_2 = \iintφ_{P1}\frac{1}{r_{12}}φ_{Q2} \, dτ_1dτ_2 = \\ \intφ_{Q2}\left(\intφ_{P1}\frac{1}{r_{12}}dτ_1\right) \, dτ_2 = I J=\intφ_{P1}\frac{1}{r_{12}}dτ_1 と置きます。 最初に、この$J$、電子1の積分を求めます。ガウス型関数を代入します。 J=\int cαe^{-αr_{A1}^2}\frac{1}{r_{12}}dτ_1 $r_{A1}, r_{12}$を代入します。解法ページの座標系A1（$r_{A2}$は定数）をそのまま使えます。 r_{A1}=(coshμ_1-cosν_1)\frac{r_{A2}}{2} \\ r_{12}=(coshμ_1+cosν_1)\frac{r_{A2}}{2} \\ J=cα\int e^{-\frac{1}{4}αr_{A2}^2(coshμ_1-cosν_1)^2}\frac{1}{(coshμ_1+cosν_1)\frac{r_{A2}}{2}}dτ_1 体積素片を代入します。 dτ_1=(\frac{r_{A2}}{2})^3sinhμ_1sinν_1(cosh^2μ_1-cos^2ν_1) \, dμ_1dν_1dφ_1 \\ J=cα\int e^{-\frac{1}{4}αr_{A2}^2(coshμ_1-cosν_1)^2}\frac{1}{(coshμ_1+cosν_1)\frac{r_{A2}}{2}} \\ (\frac{r_{A2}}{2})^3sinhμ_1sinν_1(cosh^2μ_1-cos^2ν_1) \, dμ_1dν_1dφ_1 整理すると J= \\ \frac{1}{4}cαr_{A2}^2\iiint e^{-\frac{1}{4}αr_{A2}^2(coshμ_1-cosν_1)^2} sinhμ_1sinν_1(coshμ_1-cosν_1) \, dμ_1dν_1dφ_1 ここでβを置きます。 β= \frac{1}{4}αr_{A2}^2 \\ J=cβ\iiint e^{-β(coshμ_1-cosν_1)^2} sinhμ_1sinν_1(coshμ_1-cosν_1) \, dμ_1dν_1dφ_1 変数$ν_1$だけの積分を分けると J=cβ\int dφ_1\int sinhμ_1 \left\{\int e^{-β(coshμ_1-cosν_1)^2}sinν_1(coshμ_1-cosν_1) \, dν_1\right\} dμ_1 変数$ν_1$だけの積分を求めます。 \int e^{-β(coshμ_1-cosν_1)^2}sinν_1(coshμ_1-cosν_1) \, dν_1 この積分は指数部分を微分してみると求まります。 \int e^{-β(coshμ_1-cosν_1)^2}sinν_1(coshμ_1-cosν_1) \, dν_1 \\ =-\frac{1}{2β}e^{-β(coshμ_1-cosν_1)^2}+C 0からπの定積分は \int e^{-β(coshμ_1-cosν_1)^2}sinν_1(coshμ_1-cosν_1) \, dν_1 \\ =-\frac{1}{2β} \left[ e^{-β(coshμ_1-cosν_1)^2} \right]_0^π \\ =-\frac{1}{2β} \left( e^{-β(coshμ_1+1)^2}-e^{-β(coshμ_1-1)^2}\right) 変数$ν_1$だけの積分結果を、元の式に代入すると J=2πcβ \int sinhμ_1 \left\{ -\frac{1}{2β} \left( e^{-β(coshμ_1+1)^2}-e^{-β(coshμ_1-1)^2}\right) \right\} dμ_1 整理すると J =-πc \int sinhμ_1 \left\{ e^{-β(coshμ_1+1)^2}-e^{-β(coshμ_1-1)^2} \right\} dμ_1 \\ =-πc \int sinhμ_1 \left\{ e^{-β(cosh^2μ_1+2coshμ_1+1)}-e^{-β(cosh^2μ_1-2coshμ_1+1)} \right\} dμ_1 \\ =-πc \int sinhμ_1 \left\{ e^{-βcosh^2μ_1}e^{2βcoshμ_1}e^{-β}-e^{-βcosh^2μ_1}e^{-2βcoshμ_1}e^{-β} \right\} dμ_1 \\ =-2πce^{-β} \int e^{-βcosh^2μ_1} \left\{ \frac{e^{2βcoshμ_1}-e^{-2βcoshμ_1}}{2} \right\} sinhμ_1 \, dμ_1 \\ sinhχ=\frac{e^{χ}-e^{-χ}}{2} さらにハイパー関数を考えると、$J$は次のように書き換えれます。 J =-2πce^{-β} \int e^{-βcosh^2μ_1} sinh(2βcoshμ_1) sinhμ_1 \, dμ_1 積分部分だけを求めます。 \int e^{-βcosh^2μ} sinh(2βcoshμ) sinhμ \, dμ さらに$ξ=coshμ$と置くと、$dξ=sinhμ \, dμ$ \int e^{-βξ^2} sinh(2βξ) \, dξ 難解になってきました・・・・ 以上、ガウス型で少し試行しました。スレーター型で解いてある解法ページは、よく知られた解であり、解に間違いはありません。しかしハイパー関数で解いてあるのは初めて見ました。それゆえ解法ページのMyLibraryをガウス型に変更するのは難解です。非経験的な解法のプログラムに組み込みたかったのですが。 突然ですが、このシリーズは今回で最終回です。短い間でしたが、最後までご覧いただきありがとうございました。

pythonデータ分析試験の概要 一般社団法人Pythonエンジニア育成協会が始めた認定試験。 ぶっちゃけて言うと、pythonエンジニア基礎認定試験の続き？のような資格試験です。 データ分析の部分に特化したpython試験だと思って間違えないかと。 ・受検日：通年 ・受検方法：全国のオデッセイコミュニケーションズCBTテストセンター ・受験料：一万円（税別）、ただし学生は5000円（税別） ・合格ライン：７０％ pythonエンジニア基礎認定試験とほぼ同じです。 （もちろん試験範囲は異なりますが。） 受験しようとしたきっかけ 一番の理由は、E資格試験を、勉強する上でpythonの知識をもっと極めないと前に進まないといった状況から。もう一つは、pythonのプログラミング能力を上げるのに、必要かなぁと感じたことから。 （競プロとかやってみたいなーとほのかな野望を持ちつつ） 勉強方法 あとにあげる推薦図書を写経したのち、ウェブで公開されてる模擬問題を解きつつ、 弱い部分を実際にコードを組んで頭に叩きこんでました。 使用した模擬問題サイト 1.PRIME STUDY 2.DIVE INTO EXAM 平均して７〜８割取得できるようになるまでやりこみました。 模擬試験を行うには、無料ですがメールアドレスの登録が必要。 模擬試験以外の問題もたくさんありますが、模擬試験でいいスコアを取れるように していけば大丈夫かなと思います。 （不安な人は、下記にあげる分析ノックもやってみると理解が深まるかなと） 参考図書 おすすめはやはり、推薦図書に挙げられているこの１冊 Pythonによる新しいデータ分析の教科書 この一冊をやりこめば、大丈夫かと。 上記の一冊でもの足りなかったら下記の一冊もおすすめ。 Python実践データ分析100本ノック 実際にコードを動かすので、頭に叩きこみやすい。 あると良さげな知識 言わずかもがな,Pythonエンジニア基礎認定試験の知識と、所々に出てくるG検定の知識でした。 分析方法なんかは、G検定の試験範囲とかぶるところが出てきます。 感想 E資格試験でぼんやりとわからなかった部分が少しずつ、わかるようになってきました。 まだまだE資格試験までたくさん山がありますが、こなしていけるようになりたいです。 無事合格できてよかったです。

はじめに キューにグローバルAISACを設定した後とかに、違うAISACに差し替えたいといったときに、 AISACを外すスクリプトを作ってみた グローバルAISACを複数のキューに追加する場合は、複数キューを選択してグローバルAISACへドロップ（なんか逆な感じがするけど） することでまとめて設定できるのですが、 間違えてしまったり、　あとで変更したいと思ったときに　Undoができないところまで進行していると、削除するのに苦労する場合があります。 （ツリーで消したいキューを選択して、リストでLinkAisacを選択して消すとかはできそうですが、数が多いと選択がそもそも大変だったり） 選択に関しては 自分のやりかただと、コメントに特定の文字列を入れて、リストエディターでフィルターして表示とかしているのですが、 ここで、キューまでは表示できるのですが、 AisacLinkの場合、キューの下の項目なので、リストのフィルターでも表示ができない・・・ そこで、ここまでキューのフィルターができている状態で、キューを選択して、 AisacLinkを消すというスクリプトになります。 スクリプト # --Description:[tatmos][Remove]選択しているキューのAisacLinkをすべて削除 import sys import cri.atomcraft.project import cri.atomcraft.project as acproject import cri.atomcraft.debug as acdebug import cri.atomcraft.preview as acpreview # 選択しているキューを得る selected_Cues = acproject.get_selected_objects("Cue")["data"] if not selected_Cues : acdebug.warning("キューを選択してください.") sys.exit() for cue in selected_Cues: aisacLinks = acproject.find_objects(cue, "AisacLink")["data"] for aisacLink in aisacLinks: acproject.delete_object(aisacLink)

Grove Multi Channel Gas sensor (v1) Grove - Multichannel Gas Sensor ・Cで書かれたサンプルコードをAdruinoに流し込んだ時はうまく動くもののPythonで書かれたサンプルコードをRaspberryPi3Bで動かすと認識しない ・Cのコードを読んでPythonに書き直そうと頑張ってみるも無理 ・I2Cアドレス初期値は0x04、接続に使ったGrove Base Hat for Raspberry Pi Zeroも初期値は0x04らしいことがadc.pyに書いてあった（さっき気づいた） ・いろいろ測定できるように見えてセンサはCO(とその他諸々)、NOx、NH3の３つしかない(自動車の排ガス用センサかな？) ・結局諦めて他のセンサを買った MQシリーズとADコンバータ ・ガスセンサでぐぐると一番最初にヒットするMQシリーズを買ってみた ・ADコンバータが必要になるようなのでADS1015を買った ・ADS1015はサンプルコードがgitにあるので簡単に動いた ・MQ用のブレークアウト基盤はあった方がいい ・12bit(1024)はちょっと解像度が物足りないけどゲインを上げると壊しそうで怖い ・CO2センサーと違って校正済みではないのでガス校正が必要 ・ゼロ(N2)とスパン(最大値)でいいかな… ・この辺までは先月にやったことなので詳細を忘れた(思い出したら書く) Grove O2 sensor (MIX8410) ・センサー関係の実装を後回しにしているうちに酸素センサーも発売されていたので買ってみた ・公式のサンプルコードだとBase HatのADコンバータを介して値を取得する仕様のよう ・モジュールとしてインポートしてそのまま使えるのでは？→だめでした ・・値を取得するMix8410()が@propertyで守られていて使えないしコメントアウトしても動かない→なぜ？？？そもそもプログラミングの基本を習っていないので全然わからん。 ・サンプルコードの中でBaseHatのADコンバータの値を取得するモジュールが使用されていたのでadc.pyの中身を読んでみるとコメント内にサンプルコードが書かれていたのでそれを参考にした ・・このときBase HutのADCのアドレスの初期設定は0x04と書かれていた（んんん？？？Multi Channel Gas Sensorの時に出てきた謎の0x04はこれじゃ…） ・Base HutのADCから直接値を取得できたので公式のサンプルコードを見ながら適当に数値の換算→公式サンプルコードを直接動かしたときとほぼ同じ動作を確認 ・MQ用のADCのチャンネルまだ余ってるしこっちでいいよね… ・2個買ってしばらく動かして様子を見ると23とか28あたりで数値が安定しているのでこちらも自前でガス校正が必要っぽい ・ゼロ(N2)とスパン(大気)でいいかな ・プロトコルは以下 ・Getting Started with Raspberry Piを参考に進めた ・センサのチュートリアルにあるSetting Softwareに従ってGrove Base Hat関係のモジュールが入ったライブラリをインストール ・チュートリアルの手順に従ってサンプルコードで動作確認 ・公式のサンプルコードをモジュールとして使うことができなかったので自分で書き直したコードが以下 Python import time from grove.adc import ADC adc = ADC() while True: value = adc.read(0) voltage = value*3.3/1024.0 Mix8410Value = Voltage*0.21*100/2.0 print(Mix8410Value) time.sleep(10) とりあえずこれでしばらく動かしてセンサが公式のサンプルコードと同じように動作するか確認した。 けど値を取得するだけなら from grove.adc import ADC adc = ADC() value = adc.read(0) #今回はBaseHatのA0チャンネルを読み取っている で十分だと思う。 ２つ買って値がずれていた上に20.9%前後ではないことからセンサに個体差がありかつ公式の換算式はあてにならないと考えられるためガス工校正が必要 結論 ・ADS1015にMQ-7, MQ-8, Mix8410を接続して値を取得する(Base Hatは使わない) ・全てのセンサーでガス校正必要 ・・ゼロガス(N2)とスパンガス(1%添加or大気)で校正をすればいけるかなぁ ・アナログセンサとADCの間ははんだ付けが必要だよねやっぱり ・・ユニバーサル基盤と3Dプリンタでセンサユニット作らなきゃ… ・センサをブレークアウト基盤に刺すところは導電グリスで誤魔化せないかな…

Grove Multi Channel Gas sensor (v1) ・Cで書かれたサンプルコードをAdruinoに流し込んだ時はうまく動くもののPythonで書かれたサンプルコードをRaspberryPi3Bで動かすと認識しない ・Cのコードを読んでPythonに書き直そうと頑張ってみるも無理 ・I2Cアドレス初期値は0x04、接続に使ったGrove Base Hat for Raspberry Pi Zeroも初期値は0x04らしいことがadc.pyに書いてあった（さっき気づいた） ・いろいろ測定できるように見えてセンサーはCO(とその他諸々)、NOx、NH3の３つしかない(自動車の排ガス用センサーかな？) ・結局諦めて他のセンサーを買った MQシリーズとADコンバータ ・ガスセンサーでぐぐると一番最初にヒットするMQシリーズを買ってみた ・ADコンバータが必要になるようなのでADS1015を買った ・ADS1015はサンプルコードがgitにあるので簡単に動いた ・MQ用のブレークアウト基盤はあった方がいい ・12bit(1024)はちょっと解像度が物足りないけどゲインを上げると壊しそうで怖い ・CO2センサーと違って校正済みではないのでガス校正が必要 ・ゼロ(N2)とスパン(最大値)でいいかな… ・この辺までは先月にやったことなので詳細を忘れた(思い出したら書く) Grove O2 sensor (MIX8410) ・センサー関係の実装を後回しにしていうるうちに酸素センサーも発売されていたので買ってみた ・公のサンプルコードだとBase HatのADコンバータを介して値を取得する仕様のよう ・モジュールとしてインポートしてそのまま使えるのでは？→だめでした ・・値を取得するMix8410()が@propertyで守られていて使えないしコメントアウトしても動かない→なぜ？？？ ・サンプルコードの中でBaseHatのADコンバータの値を取得するモジュールが使用されていたのでadc.pyの中身を読んでみるとコメント内にサンプルコードが書かれていたのでそれを流用 ・・このときBase HutのADCのアドレスの初期設定は0x04と書かれていた（んんん？？？） ・Base HutのADCから直接値を取得できたので公式のサンプルコードを見ながら適当に数値の換算→公式サンプルコードを直接動かしたときとほぼおなじ動作をした ・MQ用のADCのチャンネルまだ余ってるしこっちでいいよね… ・2個買ってしばらく動かして様子を見ると23とか28あたりで数値が安定しているのでこちらも自前でガス校正が必要っぽい ・ゼロ(N2)とスパン(大気)でいいかな ・Getting Started with Raspberry Piを参考に進めた ・Setting Softwareに従ってGrove Base Hat関係のモジュールが入ったライブラリをインストール ・公式の手順に従ってサンプルコードで動作確認 ・公式のサンプルコードをモジュールとして使うことができなかったので自分で書き直したコードが以下 Python import time from grove.adc import ADC adc = ADC() while True: value = adc.read(0) voltage = value*3.3/1024.0 Mix8410Value = Voltage*0.21*100/2.0 print(Mix8410Value) time.sleep(10) とりあえずこれでしばらく動かしてセンサーが公式のサンプルコードと同じように動作するか確認した。 けど値を取得するだけなら from grove.adc import ADC adc = ADC() value = adc.read(0) #今回はBaseHatのA0チャンネルを読み取っている で十分だと思う。センサーに個体差あって校正が必要だし。 結論 ・ADS1015にMQ-7, MQ-8, Mix8410を接続して値を取得する(Base Hatは使わない) ・ガス校正必要 ・アナログセンサとADCの間ははんだ付けが必要だよねやっぱり ・センサをブレークアウト基盤に刺すところは導電グリスで誤魔化せないかな…

はじめに 画像認識を用いて解決したいタスクとしてしばしばある画像と類似する画像の検索を行いたいというタスクが表れることがあります。そのような、画像同士の類似度を学習する仕組みの中で深層学習を利用した方法に深層距離学習というものがあります。 今回はその深層距離学習のモデルの1つの「ArcFace」というモデルに関してPyTorchによる実装コードを追いながら勉強した内容を紹介したいと思います。 今回解説するコード from __future__ import print_function from __future__ import division import torch import torch.nn as nn import torch.nn.functional as F from torch.nn import Parameter import math class ArcMarginProduct(nn.Module): r"""Implement of large margin arc distance: : Args: in_features: size of each input sample out_features: size of each output sample s: norm of input feature m: margin cos(theta + m) """ def __init__(self, in_features, out_features, s=30.0, m=0.50, easy_margin=False): super(ArcMarginProduct, self).__init__() self.in_features = in_features self.out_features = out_features self.s = s self.m = m self.weight = Parameter(torch.FloatTensor(out_features, in_features)) nn.init.xavier_uniform_(self.weight) self.easy_margin = easy_margin self.cos_m = math.cos(m) self.sin_m = math.sin(m) self.th = math.cos(math.pi - m) self.mm = math.sin(math.pi - m) * m def forward(self, input, label): # --------------------------- cos(theta) & phi(theta) --------------------------- cosine = F.linear(F.normalize(input), F.normalize(self.weight)) sine = torch.sqrt((1.0 - torch.pow(cosine, 2)).clamp(0, 1)) phi = cosine * self.cos_m - sine * self.sin_m if self.easy_margin: phi = torch.where(cosine > 0, phi, cosine) else: phi = torch.where(cosine > self.th, phi, cosine - self.mm) # --------------------------- convert label to one-hot --------------------------- # one_hot = torch.zeros(cosine.size(), requires_grad=True, device='cuda') one_hot = torch.zeros(cosine.size(), device='cuda') one_hot.scatter_(1, label.view(-1, 1).long(), 1) # -------------torch.where(out_i = {x_i if condition_i else y_i) ------------- output = (one_hot * phi) + ((1.0 - one_hot) * cosine) # you can use torch.where if your torch.__version__ is 0.4 output *= self.s # print(output) return output 今回はArcFaceのモデルの概要について紹介し、その後上記のArcFaceのPyTorchによる実装(https://github.com/ronghuaiyang/arcface-pytorch/blob/master/models/metrics.py) でモデルの処理・数式が上記コードで実現できることを見ていきたいと思います。 ArcFaceの概要 ここではArcFaceの概要やそこに表れる数式に簡単に説明します。 ArcFaceの詳しい説明は、「モダンな深層距離学習 (deep metric learning) 手法: SphereFace, CosFace, ArcFace」や元論文の説明がわかりやすかったです。 ArcFaceはクラス分類問題において非常頻繁に用いられるSoftmax Cross Entropyという損失関数 $$ L_1 = -\frac{1}{N}\sum^N_{i=1}\log\frac{\exp({W^T_{y_i}x_i+b_i)}}{\sum^C_k\exp({W^T_{k}x_i+b_k})} \;\;\;\;\;\;\;\;\;\;\;\;\;\; (1) $$ に対して２つの工夫を施すことで得られる手法になります。 なお、式中に出てくるパラメーターの定義は以下であるとします。 N: バッチサイズ i: バッチ内の学習サンプル $W_k$: 最後の全結合層の重み$W$のk番目の列 $x_i$: 全結合層への入力となるD次元ベクトル $b_k$: バイアス $y_i$: i番目のサンプルの正解クラスID [１つ目の工夫] 重み・特徴ベクトルの正規化 ArcFaceにおける工夫の１つ目はSoftmax Cross Entropyの中の重みベクトル$W^T_{y_k}$と特徴ベクトル$x_i$をそれぞれ$L_2$正規化(ノルムが１にする正規化)を行うというものです。 この、正規化によって重みベクトルと特徴ベクトルはそれぞれ $W_{y_k}$->半径１の超球面上に射影されたクラス$k$の代表ベクトル $x_i$->半径１の超球面上に射影されたi番目のデータの特徴ベクトル と解釈することができます。 この時、クラス$k$の代表ベクトルとi番目のデータの特徴ベクトルがなす角度$\theta_{k}$とすると、それらの内積は $$ W^T_{k}x_i=|W^T_{k}||x_i|\cos\theta_{k}=\cos\theta_{k} $$ となります。 さらに、バイアス項$b_k=0$とすることで$L_2$正規化を施したCross entropy lossは $$ L_2 = -\frac{1}{N}\sum^N_{i=1}\log\frac{\exp(\cos\theta_{y_i})}{\sum^C_k\exp(\cos\theta_{k})} \;\;\;\;\;\;\;\;\;\;\;\;\;\; (2) $$ と表せます。 この式の意味するところは、各クラスの代表ベクトル$W_k$と特徴ベクトル$x_i$のなす角$\theta_k$が最も小さいクラスに分類する損失を表します。この工夫によって代表ベクトル$W_k$や特徴ベクトル$x_i$距離の近さを比べたいベクトル同士の成す角度の近さと解釈できるようになります。 [2つ目の工夫]角度に対するマージン 次に行うのはAngular Margin Penaltyと呼ばれる工夫です。その工夫とは「(2)式において、正解クラス$y_i$に対応する$\cos\theta_{y_i}$のみ」に対して、$\cos(θ_{y_i}+m) (m>0)$ という変更（ペナルティ）を加えるというものです。 こうようにすることで正解クラスの代表ベクトル$W_{y_i}$とのなす角$θy_i$が他のどのクラスの代表ベクトルとの角度よりも、$θ_{y_i}$がm以上のマージンを持って小さくなるように学習されます。 このペナルティにより、モデルは特徴ベクトルxのクラス内分散を小さくし、クラス間分散を大きくして、このマージンを稼ぐように学習が進みます。 最後にこのcosθの値はs (s>1)倍されてsoftmaxへ送られます。これは$\cos\theta_k$の値が小さすぎるとsoftmaxが機能しなくなるために調節を行っている処理になります。 以上のことをまとめると、ArcFaceのロスを書き下すと下記になります。 $$ L_3 = -\frac{1}{N}\sum^N_{i=1}\log\frac{\exp(s\cos(\theta_{y_i}+m)}{\exp(s\cos(\theta_{y_i}+m)+\sum^C_{k≠y_i}\exp(s\cos(\theta_k))} \;\;\;\;\;\;\;\;\;\;\;\;\; (3) $$ 上記のような処理を加えたモデルを通常のクラス分類として学習させることで、特徴ベクトルxのクラス内分散を小さくし、クラス間分散を大きくするような学習が実現されます。実際に利用する場合には、正規化された特徴ベクトル$x_i$を抽出し、それらのコサイン類似度によってサンプル間の類似度を算出することができます。 コード解説 ① コードの前半部分（コンストラクターの定義まで） from __future__ import print_function from __future__ import division import torch import torch.nn as nn import torch.nn.functional as F from torch.nn import Parameter import math class ArcMarginProduct(nn.Module): r"""Implement of large margin arc distance: : Args: in_features: size of each input sample out_features: size of each output sample s: norm of input feature m: margin cos(theta + m) """ def __init__(self, in_features, out_features, s=30.0, m=0.50, easy_margin=False): super(ArcMarginProduct, self).__init__() self.in_features = in_features self.out_features = out_features self.s = s self.m = m self.weight = Parameter(torch.FloatTensor(out_features, in_features)) nn.init.xavier_uniform_(self.weight) self.easy_margin = easy_margin self.cos_m = math.cos(m) self.sin_m = math.sin(m) self.th = math.cos(math.pi - m) self.mm = math.sin(math.pi - m) * m ここでは、コンストラクターでこれから使う変数の宣言・初期化などをしています。 また、コンストラクターの要求する引数は in_features: 特徴ベクトルの次元数 out_features: 分類クラスするクラス数 s: スケール因子 m: マージン easy_margin: 後述 を意味します。 ② cosの計算 def forward(self, input, label): # --------------------------- cos(theta) & phi(theta) --------------------------- cosine = F.linear(F.normalize(input), F.normalize(self.weight)) このdef forwardからがArcFaceのアルゴリズムが実装されている箇所になります。 まず、関数の定義部分を見るとArcFaceの前段の層からの出力(input)と正解ラベル(label)が引数であることがわかります。 そのあとF.linear(F.normalize(input), F.normalize(self.weight))という量を計算しています。 この処理は、 F.normalize()という処理が$L_2$正規化の処理で、 F.normalize(input)で入力データの$L_2$正規化 F.normalize(self.weight)で各クラスの代表ベクトルの$L_2$正規化 ​を実施します。 その後、1.で計算したそれぞれをF.linear()という$W^T x_i$を計算する関数に引数として渡しています。 この時2.の結果として[$\cos\theta_1$,・・・,$\cos\theta_{y_i}$,$\cos\theta_C$]という入力データと各クラスの代表ベクトルとのcosを要素として持つ配列が計算されます。 ③ マージンの適用 sine = torch.sqrt((1.0 - torch.pow(cosine, 2)).clamp(0, 1)) phi = cosine * self.cos_m - sine * self.sin_m 次の処理の準備のためにsinを計算します。 三角関数の加法定理$cos(A+B)=cosAcosB-sinAsinB$を使いcosineに対して各要素の角度からmラジアンを足した配列を$\phi$作ります。 $$ \phi = [\cos\theta_1 + m,・・・,\cos\theta_{y_i} + m,\cos\theta_C + m] $$ ④ cosの周期性を考慮した処理。 if self.easy_margin: phi = torch.where(cosine > 0, phi, cosine) else: phi = torch.where(cosine > self.th, phi, cosine - self.mm) $\cos\theta$は$\theta=\pi$を境に減少から増加に転じます。ここでマージンmの目的は入力データと答えのクラスの代表ベクトルを遠ざけることが目的であったことを思い返すと、$\cos(\theta + m)<\cos\theta$である必要があります。そのため、$\theta>\pi-m$の範囲でマージンをとった場合上記の問題が発生する可能性があります。 そのため、$\theta>\pi-m$の場合はマージンの取り方を変更し角度ではなくcosそのものにマージンを適用します（cosarcのマージンのとりかた）。この処理がelseの中身です。 一方上記のように複雑に考えずにcosine < 0になったらマージンの適用をやめるという方法で上記の問題に対処するという方針をとるのがself.easy_margin=Trueの場合です。 ⑤ 正解ラベルをone hotベクトル化する。 one_hot = torch.zeros(cosine.size(), device='cuda') one_hot.scatter_(1, label.view(-1, 1).long(), 1) ⑥ 配列cosineに関して正解ラベルの要素の部分を配列$\phi$のものに置き換える。 output = (one_hot * phi) + ((1.0 - one_hot) * cosine) # you can use torch.where if your torch.__version__ is 0.4 ⑦ 最後sでスケールする。 output *= self.s 利用方法 最後に既存のPyTorchで実装されたCNNのモデルがある場合、それに対してArcFaceを適用する方法を紹介したいと思います。 ① ソースコードのダウンロード 以下からソースコードをダウンロードします。 https://github.com/ronghuaiyang/arcface-pytorch/blob/master/models/metrics.py ② ソースファイルのインポート import metrics ③ ArcFaceのインスタンス生成 例： input_data=10、分類クラス数10、 s=30.0, m=0.05, easy_margin=Trueでインスタンス生成する場合 metric = metrics.ArcMarginProduct(10, 10, s=30.0, m=0.05, easy_margin=True) ④ Optimizerのインスタンス生成部分の変更 ArcFaceの重みパラメーターをOptimizerで更新するためにOptimizerのインスタンス生成部分にArcFaceの部分を記述します。 例：既存のモデルがmodelという名前で最適化関数がAdamの場合 #optimizer = optim.Adam([{'params': model.parameters()}, {'params': metric.parameters()}], lr=0.00001) ⑤ 学習部分に関して 損失関数の計算～誤差逆伝搬によるパラメーターの更新の部分は以下のように実装できます。 例：損失関数の名前がcriterionである場合 ・・・割愛(epoch数に関するfor文やdataloaderからdataやlabelを取得する処理等)・・・ # 前回の勾配情報をリセット optimizer.zero_grad() # 予測 features = model(inputs) outputs = metric(features, labels) # 予測結果と教師ラベルを比べて損失を計算 loss = criterion(outputs, labels) # 損失に基づいてネットワークのパラメーターを更新 loss.backward() optimizer.step() 注意する点としては損失関数の入力となるoutputsを予測するための処理が まず、既存の分類モデルであるmodelに入力データを代入して特徴ベクトルfeatures得る。 次に、出力ベクトルと正解ラベルをArcFaceに代入しoutputsを得る。 という2ステップになっている点です。 以上のようにして既存の分類モデルにArcFaceを組み込むことができます。 おわりに 今回はArcFaceという深層距離学習の理論・実装・使い方を確認しました。 今回記事を書くにあたり元論文を読んだりもしたのですが難解な数式や背景知識などがあまり現れずかなり読みやすく理解しやすいという印象をうけました。また、今回紹介した実装コードも100行に満たない簡潔で理解し易いコードと感じました。 このように、ArcFaceは導入にあたって必要なラーニングコストが少ないアルゴリズムであると感じました。是非みなさんも機会があれば活用してみてはいかがでしょうか？

scikit-learnで線形回帰をしていたら、切片がとんでもない大きさ（2e+11）になった。 どうも、特徴量に大量にあったダミー変数（業種区分のようなもの）と、数値型の標準化が悪さをしていたらしい。 それの再現と対策を書く。 線形回帰はモデルの解釈がしやすい、という点で有用な一方で、（予測単体ならともかく）まともに解釈ができる状態にもっていくまでに色々な落とし穴がある。 ダミー変数は１つで、OneHot encodingすると出てくるものとする。 方針：何をするのか ボストン住宅価格データ + 適当なダミー変数（１変数）で、上記の問題を再現する。 コード：トラブルの再現 import pandas as pd import numpy as np from sklearn.linear_model import LinearRegression from sklearn.preprocessing import StandardScaler from sklearn.datasets import load_boston どれか1つでフラグが立つようなダミー変数が、フィッティングパラメータに与える影響をみたいので CHAS列を削除 dummy変数の列を用意 それをOneHot enncoding data = load_boston() df = pd.DataFrame(data['data'], columns=data['feature_names']) df['target'] = data['target'] df.drop('CHAS', axis=1, inplace=True) str_list = [f'A{i}' for i in range(20)] df['dummy'] = np.random.choice(str_list, len(df)) df = pd.get_dummies(df) あとはダミー変数以外の入力と、ターゲットをそれぞれで標準化して、LinearRegressionでフィッティングする。 num_fe_cols = df.columns.to_list() ex_fe = [f'dummy_{i}' for i in str_list] ex_fe += ['target'] num_fe_cols = [fe for fe in num_fe_cols if fe not in ex_fe] # 標準化 x_std_scaler = StandardScaler() y_std_scaler = StandardScaler() inputs = x_std_scaler.fit_transform(df[num_fe_cols]) targets = y_std_scaler.fit_transform(df[['target']]) # 念のため、DataFrameを分けてから入力作成（ダミー変数付き） df_std = df.copy() df_std.loc[:, num_fe_cols] = inputs df_std.loc[:, 'target'] = targets feature_names = df_std.columns.to_list() feature_names.remove('target') # フィッティング lm = LinearRegression() lm.fit(df_std.loc[:, feature_names], targets) すると、係数はこのような結果……。ダミー変数に対応する部分の係数が"-2.87597484e+11"というとんでもない値になっている。 lm.coef_ array([[-1.23998168e-01, 1.11775353e-01, 3.20912206e-02, -2.20312936e-01, 2.86927225e-01, -4.96117784e-03, -3.37975459e-01, 3.19861003e-01, -2.37376673e-01, -2.38160875e-01, 9.27067463e-02, -4.29321436e-01, -2.87597484e+11, -2.87597484e+11, -2.87597484e+11, -2.87597484e+11, -2.87597484e+11, -2.87597484e+11, -2.87597484e+11, -2.87597484e+11, -2.87597484e+11, -2.87597484e+11, -2.87597484e+11, -2.87597484e+11, -2.87597484e+11, -2.87597484e+11, -2.87597484e+11, -2.87597484e+11, -2.87597484e+11, -2.87597484e+11, -2.87597484e+11, -2.87597484e+11]]) 切片を見ると lm.intercept_ array([2.87597484e+11]) 係数とちょうど逆符号になっている。何が起こったのだろう。 原因：どうしてこうなったのか この現象が起きる原因は 目的変数は標準化されて平均が0になっている ダミー変数が全て0になることはない そのため、両辺で平均を取ったときに、$0 = 切片 + w_n$($w_n$をn番目のダミー変数に対応する係数)という関係になっている 切片とダミー変数用のパラメータは、逆符号かつ値が決まらないので、非常に大きくなることがある おそらく、このような理由だと考えている。 実際、何回か計算をさせると、係数が変わったり、一見まともな値になることもあった。 切片を一意に決められないことから来る計算の不安定さが、こういった現象の原因だろう。 対策：どうすればいいのか 値が決まらないなら、切片を決めてしまえばよい、ということで切片=0にしてしまう。 lm = LinearRegression(fit_intercept=False) # デフォルトはTrue これで解決する。 lm.coef_ array([[-0.12400473, 0.1117352 , 0.03217197, -0.22048052, 0.2869507 , -0.00498836, -0.33796915, 0.319808 , -0.23725577, -0.23817404, 0.09274132, -0.42924623, 0.08038865, 0.0622342 , 0.09501951, -0.08434313, 0.11698025, 0.02396338, 0.08164036, 0.0452252 , -0.19803361, -0.10092905, -0.07053491, 0.03973364, -0.04118562, -0.20115119, -0.0753999 , 0.16343603, -0.08174667, 0.22527764, -0.10547045, -0.18457662]]) lm.intercept_ 0.0 他には pd.get_dummies(df, drop_first=True)にする（オールゼロのダミー変数が作られるので、上記の問題が生じない） 目的変数が平均がゼロにならないような変換にする（MinMaxScalerなど） などの対策がある（が、本質的ではないと思う）。 今回のデータセットだと、ダミー変数が10以上でこの現象が起きやすくなり、それ未満だとあまり起こらないことが多かった。 エラーが出ないからといって、案外見過ごしているかもしれないので注意が必要。 多重共線性とは似て非なるというか、数値計算特有のトラブルであるといった印象を持った。

概要 pythonでgnuradioやってみた。 AM変調を録音。 信号は、440ｈｚのサイン波、搬送波は、5KHZ。 サンプルコード from gnuradio import gr, blocks, analog class my_top_block(gr.top_block): def __init__(self): gr.top_block.__init__(self) samp_rate = 32000 self.blocks_wavfile_sink_0 = blocks.wavfile_sink('a2.wav', 1, samp_rate, 8) self.blocks_multiply_xx_0 = blocks.multiply_vcc(1) self.blocks_multiply_const_vxx_0 = blocks.multiply_const_vff((0.3, )) self.blocks_head_0 = blocks.head(gr.sizeof_gr_complex*1, 64000) self.blocks_complex_to_real_0 = blocks.complex_to_real(1) self.blocks_add_const_vxx_0 = blocks.add_const_vcc((1, )) self.analog_sig_source_x_1 = analog.sig_source_c(samp_rate, analog.GR_COS_WAVE, 5000, 1, 0) self.analog_sig_source_x_0 = analog.sig_source_c(samp_rate, analog.GR_COS_WAVE, 440, 1, 0) self.connect((self.analog_sig_source_x_0, 0), (self.blocks_add_const_vxx_0, 0)) self.connect((self.analog_sig_source_x_1, 0), (self.blocks_multiply_xx_0, 1)) self.connect((self.blocks_add_const_vxx_0, 0), (self.blocks_multiply_xx_0, 0)) self.connect((self.blocks_complex_to_real_0, 0), (self.blocks_multiply_const_vxx_0, 0)) self.connect((self.blocks_head_0, 0), (self.blocks_complex_to_real_0, 0)) self.connect((self.blocks_multiply_const_vxx_0, 0), (self.blocks_wavfile_sink_0, 0)) self.connect((self.blocks_multiply_xx_0, 0), (self.blocks_head_0, 0)) if __name__ == '__main__': try: my_top_block().run() except KeyboardInterrupt: pass 以上。

【AWS】【超初心者ハンズオン】 シリーズ シリーズ概要 前の記事: 2. AWS API Gatewayによって、URLへのアクセスをトリガーとしてLambdaを実行する この記事の目的 パラメータを含んだURLにアクセスした際のHttpGetRequestをトリガーとして、AWS Lambdaを実行する API GatewayとLambdaの間での変数の受け渡しのイメージを理解する 手順 URLにパラメータを付加して、GETメソッドを行う GETを行った時のログを見る json形式の応答からパラメータを取り出し、Lambdaに入力する AWS API Gatewayってどういうサービスなの？（おさらい） 前回、API Gatewayを触って、"URLアクセスをトリガーにして、AWS Lambdaを実行する"ということを行いました。 API Gatewayが"HTTP APIやWebsocket APIによって、AWSに展開された他のサービスをトリガーすること"ができるサービスだということの片鱗を理解できたのではないでしょうか。 Amazon API Gateway は、あらゆる規模の REST、HTTP、WebSocket API を作成、公開、管理、モニタリング、保護するための AWS のサービスです。 Amazon API Gateway とは何ですか? URLにパラメータを付加して、GETメソッドを行う 参考: ページリンクのURLにパラメーターをつける ここではURLにて、min、maxのパラメータを指定した上で、URLにアクセスします https://hogehoge.execute-api.ap-northeast-1.amazonaws.com/random_integer_generator_001?min=100&max=1000 min = 100 max = 1000　としてみました.> 結果: 0~9までの整数しか返ってきません 原因: URLに付加されたパラメータを、Lambda側で受け取っていないため GETを行った時のログを見る URLを開いたとき、API GatewayにはGetRequestが送られ、Lambdaがトリガーされます。 Lambda内には、どのようなRequestが送られたのかをログとして表示する部分を仕込んでいるため、ログが見れるはずです。 ログはどこにprintされているかというと... 0~9の間の整数をランダムに1つ返す import json def lambda_handler(event, context) -> int: print("Received event: " + json.dumps(event, indent=2)) ログを見てみよう 下図、'モニタリング' -> 'CloudWatchのログを表示'と開いてください。 'ログストリーム'にて、これまでLambdaを実行したときのログが表示されます。 さきほどURLに $min=100&max=1000 と付加した際のログを見ると、下記のような部分がjsonオブジェクトに含まれていることがわかります。 queryStringParameters"というkeyの中に、"max", "min"が"というkeyがある "max", "min"に対応するvalueは、int型ではなくstr型になる { "queryStringParameters" : { "max": "1000", "min": "100" } } json形式の応答からパラメータを取り出し、Lambdaに入力する Lambdaに対して、URLに付加されたパラメータを入力できるよう、Lambdaを以下のように変更しました。 -lambda_handlerはeventとしてjsonオブジェクトを受け取る -logから明らかなように、event内に"queryStringParameters"が、その中に"min", "max"がある -忘れずにint型に変換する URLに受け取ったmin,maxの範囲の整数を返す import json import random print('Loading function') def lambda_handler(event, context) -> int: print("Received event: " + json.dumps(event, indent=2)) int_min = int(event["queryStringParameters"]["min"]) int_max = int(event["queryStringParameters"]["max"]) generated_int = random.randrange(int_min, int_max, 1) return generated_int もう一度URLにパラメータを付加して、GETメソッドを行う URL呼び出し結果 190 狙い通り、URLに付加したmin, maxの範囲の値が返ってきましたね。 この先やること URLに手入力でパラメータを与えるのは流石にアレなので、ユーザーフォームを作成してmin, maxの値を入力できるようにする

はじめに Pandasを使っていて、pd.DataFrameをDBにINSERTする際にto_sql()メソッドでやりたいことがうまく出来なかったのが解決できたのでメモとして残す。 ※この方法で数万レコードをINSERTしようとすると時間がかなりかかったので要注意 やりたいこと pd.DataFrameをDBにINSERTする際にNotNull制約違反や主キー制約違反が発生した場合にその行をスキップしてINSERTしたい。 具体的には、以下のようなデータフレームを主キー制約がcompanyにNotNull制約がestablished_inにかけられたテーブルにto_sql()メソッドでINSERTしようとするとエラーになり全くINSERTされないが、INSERTできるレコードはINSERTして失敗したレコードはINSERTしないようにしたい。 company established_in founder 0 A 1900.0 Alice 1 B 1925.0 Bob 2 C NaN Charlie 3 A 1995.0 Ahn 4 D 2010.0 David ちなみに、普通に、 with engine.begin() as conn: df_fail.to_sql('company', if_exists='append', index=False, con=conn) してあげると、 (psycopg2.errors.NotNullViolation) null value in column "established_in" violates not-null constraint DETAIL: Failing row contains (C, null, Charlie). [SQL: INSERT INTO company (company, established_in, founder) VALUES (%(company)s, %(established_in)s, %(founder)s)] [parameters: ({'company': 'A', 'established_in': 1900.0, 'founder': 'Alice'}, {'company': 'B', 'established_in': 1925.0, 'founder': 'Bob'}, {'company': 'C', 'established_in': None, 'founder': 'Charlie'}, {'company': 'A', 'established_in': 1995.0, 'founder': 'Ahn'}, {'company': 'D', 'established_in': 2010.0, 'founder': 'David'})] (Background on this error at: http://sqlalche.me/e/14/gkpj) となり、全くINSERTされない。 ソースコード JupyterLabで以下を実行した。なお、DB側でcompanyには主キー制約がestablished_inにはNotNull制約がかけられている。 import numpy as np import pandas as pd import psycopg2 import sqlalchemy from IPython.display import display engine = sqlalchemy.create_engine(url='postgresql://postgres:postgres@192.168.56.101/mydb1') df_fail = pd.DataFrame({'company': ['A', 'B', 'C', 'A', 'D'], 'established_in': [1900, 1925, np.nan, 1995, 2010], 'founder': ['Alice', 'Bob', 'Charlie', 'Ahn', 'David' ]}) df_fail # 主キー制約違反・NotNull制約違反 company established_in founder 0 A 1900.0 Alice 1 B 1925.0 Bob 2 C NaN Charlie 3 A 1995.0 Ahn 4 D 2010.0 David for row in range(len(df_fail)): record = df_fail.iloc[[row]] try: with engine.begin() as conn: record.to_sql('company', if_exists='append', index=False, con=conn) print('Pass') except sqlalchemy.exc.IntegrityError as e: print('IntegrityError') print(e) except Exception as e: print('Error') print(e) finally: display(record) print('\n\n') Pass company established_in founder 0 A 1900.0 Alice Pass company established_in founder 1 B 1925.0 Bob IntegrityError (psycopg2.errors.NotNullViolation) null value in column "established_in" violates not-null constraint DETAIL: Failing row contains (C, null, Charlie). [SQL: INSERT INTO company (company, established_in, founder) VALUES (%(company)s, %(established_in)s, %(founder)s)] [parameters: {'company': 'C', 'established_in': None, 'founder': 'Charlie'}] (Background on this error at: http://sqlalche.me/e/14/gkpj) company established_in founder 2 C NaN Charlie IntegrityError (psycopg2.errors.UniqueViolation) duplicate key value violates unique constraint "constraint_pkey" DETAIL: Key (company)=(A) already exists. [SQL: INSERT INTO company (company, established_in, founder) VALUES (%(company)s, %(established_in)s, %(founder)s)] [parameters: {'company': 'A', 'established_in': 1995.0, 'founder': 'Ahn'}] (Background on this error at: http://sqlalche.me/e/14/gkpj) company established_in founder 3 A 1995.0 Ahn Pass company established_in founder 4 D 2010.0 David with engine.begin() as conn: data_on_db = pd.read_sql('SELECT * FROM company;', con=conn) data_on_db company established_in founder 0 A 1900 Alice 1 B 1925 Bob 2 D 2010 David 確かに例外処理を使ってINSERTできるものだけINSERTできた（ただし、主キー制約違反の場合上の行からINSERTされてしまう点に注意）。 おわりに 失敗した行だけ集めたpd.DataFrameを作っておけば、どこがINSERTされなかったのかわかるので便利なはず。

　Jupyter Notebookを一発で開きたい。すぐ使いたい。 やること Jupyter Notebookを開く時、一度ターミナルを開く必要があります。 またディレクトリを指定して開くこともできますが、その場合はディレクトリも入力する必要があります。 これをクリック一発で行えるアプリをデスクトップに作成します。 環境 Mac automator オートメーターを使って作ります ① 検索ウィンドウに「automator」と入力し、automatorを開く。 ② メニューの「新規」から「アプリケーション」を選択。 ③ 検索窓に「a」と打ち込んで「アップルスクリプトをを実行」を選択。 ④ スクリプトを下記のように上書きする。 applescript on run {input, parameters} tell application "Terminal" delay 0.1 do script "jupyter notebook /<デフォルトで開きたいパス名>" end tell end run ⑤ 最後に名前をつけてデスクトップ等に保存する。 実行 アプリアイコンを押して実行します。 初回起動時のみ、ダイアログボックスが出ますが、２回目以降は１発で起動できるようになります。 コンフィグでパスを設定するとディレクトリ移動に制限が出るらしいので、その点でこの方法が便利かもしれません。 参考 こちらの記事を参考にしました。 MacでワンクリックでJupyter Notebookを起動させる方法 Automatorで簡単にコマンド実行アプリを作る

単回帰モデル 最小二乗法（Ordinary Least Squares）を用いて単回帰モデルとそれに伴った色々なことを出していきます。 求めたいものを知るプログラムを末尾に載っけていてそれがメインで、序盤の理論的なところはかなり怪しいです。 今度テストがあるので、その対策でプログラムを作ったので、理論的なところ、用語の使い方は怪しさ満天です。 今回のデータ import pandas as pd import matplotlib.pyplot as plt import seaborn as sns df = pd.DataFrame({ "x" : [2,3,4,5], "y" : [6.4,7.9,7.1,8.8], }) print(df) sns.set() sns.scatterplot(data=df, x="x", y="y") plt.xlim(0,10) plt.ylim(0,10) plt.show() 説明変数をx,目的変数をyとして最小二乗法で単回帰モデルを作っていきます。 散布図に対して説明性が高いいい感じの直線を引きて〜！って感じです。 ただ、直線という制約があったら、全ての点を通るのは無理そうですよね。 ですので、いい感じということになります。 つまり、一生懸命いい感じの直線を引いたとしても、多少はその直線との誤差が出ることは避けられないということです。 でも、いい感じの直線を考えたいわけです。 できるだけ誤差が小さくなるように。 単回帰モデルの式 単回帰モデルは目的変数に対して説明変数は１つだけ。 説明変数が２つ以上になると重回帰モデルになります。 今回の単回帰モデルは $$y_i = \beta_0 + \beta_1x_i + \epsilon_i\tag{1}$$ $\beta_0,\beta_1$：回帰係数 一次関数の式を思い出せば、$\beta_0$が切片、$\beta_1$が傾きに相当するところです。 切片と傾きが決まれば直線は引けたはずですから、 $\beta_0$、$\beta_1$を求めるというのが直線を引く上での大きな目的になるわけですね。 $\epsilon_i$：y軸方向の誤差（モデルで説明できなかった差） 直線という制約上全ての点を通る直線は引けないということから発生する、避けられない誤差の部分が$\epsilon_i$です。 $\epsilon_i$には仮定があって、 $\epsilon_i$〜$N(0,\sigma^2)$ としています。 平均：$0$　分散：$\sigma^2$ の正規分布に従うということです。 この仮定のおかげで後々の式変形がうまく進みます。 (1)の式を$\epsilon_i$について変形すると $$\epsilon_i =y_i - \beta_0 - \beta_1x_i $$ $\epsilon_i$がモデルで説明できなかった差であるということを考えると、 モデルを作る上では、モデルで説明できない差が小さければ小さいほど良いと考えるのが普通ですね。 ですので、各$\epsilon_i$をかき集めてきてきた和が最小になることを考えます。 ただ、そのまま$\epsilon_i$をかき集めると、符号の影響を受けてしまうので、誤差の最小の評価は二乗で考えることにします。 なんとなく最小二乗法っぽくなりました。 ここまでの話を整理して、目的は何か改めてみると $$\sum\epsilon_i^2=\sum(y_i-\beta_0-\beta_1x_i)^2$$ これを最小にしてくれる $\beta_0,\beta_1$を推定して、いい感じの直線を引こうよ！ということです。 $\beta_0,\beta_1$は偏微分で解くのですが、省略します。 『最小二乗法 偏微分』とかググればいくらでも出てきます。 最小二乗推定量 $\beta_0,\beta_1$は推定値で出てくるので、 $\hat{\beta_0},\hat{\beta_1}$と表記します。 偏微分の結論だけを書くと、推定値は $$\hat{\beta_1} = \frac{\sum(x_i-\bar{x})(y_i-\bar{y})}{\sum(x_i-\bar{x})^2}\tag{2}$$ $$\hat{\beta_0} = \bar{y}-\hat{\beta_1}\bar{x}$$ 一応、文字の確認をしておくと、 $x_i$は$x$の１つ目の数ということです。 $y_i$は$y$の１つ目の数ということです。 $\bar{x}$はxとしている数の平均です。 $\bar{y}$はyとしている数の平均です。 ここまでpythonにしてみます。 import numpy as np from scipy import stats x = np.array([2,3,4,5]) y = np.array([6.4,7.9,7.1,8.8]) # 平均 x_mu = np.mean(x) print("xの平均:{}".format(x_mu)) y_mu = np.mean(y) print("yの平均:{}".format(y_mu)) # 回帰係数の推定に必要なもの x_x_mu = x - x_mu # データと平均の差を格納 # print("x-x_mu:{}".format(x_x_mu)) numpyの配列で出てます。 y_y_mu = y - y_mu # データと平均の差を格納 # print("y-y_mu:{}".format(y_y_mu)) numpyの配列で出てます。 T_xy = np.sum(x_x_mu * y_y_mu) # 分子 T_xx = np.sum(x_x_mu**2) # 分母 T_yy = np.sum(y_y_mu**2) # 相関係数に使えるので出したおいた。 # 回帰係数の推定 beta1 = T_xy/T_xx beta0 = y_mu - beta1*x_mu print("beta1:{}".format(beta1)) print("beta0:{}".format(beta0)) print("回帰直線：y={}+{}x".format(beta0,beta1)) こんな感じで $\hat{\beta_0}=5.31$ $\hat{\beta_1}=0.64$ と推定されました。 回帰係数が推定されるとどうなるかというと、 誤差は避けられないけどいい感じの直線というのが決まるわけです。 $$\hat{y_i} = \hat{\beta_0} + \hat{\beta_1}x_i = 5.31+0.64x_i$$ $\hat{y_i}$は誤差が避けられない、いい感じ直線の式による予測値なので、誤差$\epsilon_i$は考慮されてません。 ただ、$\epsilon_i$は予測値$\hat{y_i}$と本当のデータ$y_i$の差として推定することはできそうです。 これを予測した結果、避けられなかった本当のデータとの誤差の推定値が残差$\hat{\epsilon_i}$です。 ここまでで分かっているのは、 予測値（いい感じの直線による値）: $$\hat{y_i} = \hat{\beta_0} + \hat{\beta_1}x_i$$ 残差： $$\hat{\epsilon_i}= y_i - \hat{y_i}$$ ここまで揃ったら、次は推定したいい感じの直線の精度を考えていきます。 余談：最小二乗推定量と相関係数 $$\hat{\beta_1} = \frac{\sum(x_i-\bar{x})(y_i-\bar{y})}{\sum(x_i-\bar{x})^2}\tag{2}$$ この式がちょっと相関係数を求める式 $$r = \frac{Cov[X,Y]} {\sqrt{V[X]}\sqrt{V[Y]}} =\frac{\frac{1}{n}\sum(x_i-\bar{x})(y_i-\bar{y})} {\sqrt{\frac{1}{n}\sum(x_i-\bar{x})^2}\sqrt{\frac{1}{n}\sum(y_i-\bar{y})^2}} =\frac{\sum(x_i-\bar{x})(y_i-\bar{y})}{\sqrt{\sum(x_i-\bar{x})^2}\sqrt{\sum(y_i-\bar{y})^2}}\tag{3} $$ に似てます。上が共分散っぽくなってて。 $r = U × \hat{\beta_1}$ を解くと $U＝\frac{\sqrt{\sum(x_i-\bar{x})^2}}{\sqrt{\sum(y_i-\bar{y})^2}}$ となるので、 $$r=\frac{\sqrt{\sum(x_i-\bar{x})^2}}{\sqrt{\sum(y_i-\bar{y})^2}}\hat{\beta_1}\tag{4}$$ という関係があることがわかります。 データが与えられているのであれば、相関係数はただデータを食わせれば出るので、どんな際に有用なのか、知っているとお得なことがあるのか、私にはわかりません。 # 相関係数 r_xy_1 = T_xy / np.sqrt(T_xx*T_yy) # (3)の式 print("相関係数：{}".format(r_xy_1)) r_xy_2 = np.sqrt(T_xx/T_yy)*beta1 # (4)の式 print("相関係数：{}".format(r_xy_2)) r_xy_3 = np.cov(x,y,ddof = 1)[0,1] / np.sqrt((np.var(x, ddof = 1)*np.var(y, ddof =1))) # 共分散と分散から print("相関係数：{}".format(r_xy_3)) # 相関係数行列 print("相関係数行列") print(np.corrcoef(x,y)) # データを食わせるだけ 当たり前ですが、全部同じ値が出ます。 モデルの精度を考える 準備1 単回帰モデル： $$y_i = \beta_0 + \beta_1x_i + \epsilon_i$$ 予測値（いい感じの直線）： $$\hat{y_i} = \hat{\beta_0} + \hat{\beta_1}x_i$$ 残差： $$\hat{\epsilon_i}= y_i - \hat{y_i}$$ 重相関係数： $$Corr(y,\hat{y})$$ これらを使います。 # 予測値 y_hat = beta0 + beta1*x print("予測値:{}".format(y_hat)) # 残差 residual = y - y_hat residual print("残差:{}".format(residual)) # 重相関係数 r_yy_hat = np.corrcoef(y,y_hat)[0,1] print("重相関係数:{}".format(r_yy_hat)) 準備2 回帰変動・残差変動・全変動を考えます。 変動なので、符号の影響を受けたくないので、全部二乗で考えます。 回帰変動（予測値）： $$\sum(\hat{y_i}-\bar{y})^2$$ 自由度は説明変数の個数$k$ 残差変動（残差）： $$\sum(\hat{\epsilon_i}-\bar{\epsilon})^2 = \sum\hat{\epsilon_i}^2$$ 右の形を知っておけば良いと思いますが、 ここで$\epsilon_i$〜$N(0,\sigma^2)$の仮定が生きて、$\bar{\epsilon}=0$、になって左から右に進んでる感じですかね。 ここで残差変動だけ$\bar{\epsilon}$が記載されてないのは仮定のお陰ですかね。 自由度は $n-k-1$ $n$はデータの個数です。 （回帰変動の自由度）＋（残差変動の自由度）=（全変動の自由度） となるので、これで求めても良いです。 プログラムではこっちで求めてます。 全変動（応答）： $$\sum({y_i}-\bar{y})^2 $$ 自由度はn-1 # 回帰変動 reg_v = np.sum((y_hat-y_mu)**2) print("回帰変動：{}".format(reg_v)) # 残差変動 res_v = np.sum(residual**2) print("残差変動：{}".format(res_v)) # 全変動 total_v = np.sum((y-y_mu)**2) print("全変動：{}".format(total_v)) # 回帰変動の自由度 reg_free = x.reshape(-1,1).shape[1] # 説明変数xの配列を行列にして列の数をカウントしてます。 print("回帰変動の自由度：{}".format(reg_free)) # 全変動の自由度 total_free = len(x)-1 print("全変動の自由度：{}".format(total_free)) # 残差変動の自由度 res_free = total_free - reg_free print("残差変動の自由度：{}".format(res_free)) 実は自由度の個数だけでなく、 （回帰変動）＋（残差変動）=（全変動） になっていることも分かります。 式を分解するとそうなります。 決定係数・自由度調整済み決定係数 準備ができたところで、モデルの評価をしていきます。 単回帰の場合は説明変数の選択がxのみ１通りなので、普通に決定係数で良いと思います。 自由度調整済み決定係数にすると、重回帰モデルで説明変数の個数が違うもの同士の精度を比較することができます。統計検定2級の問題でよく出てました。 決定係数 $$R^2=\frac{回帰変動}{全変動}=1-\frac{残差変動}{全変動}$$ $0≦R^2≦1$で評価されて、1であればあるほどモデルが良いということです。 説明変数を増やすほど$R^2$の値を大きくすることができてしまうという欠点があります。 それに対抗するのが自由度調整済み決定係数です。 重相関係数の2乗が決定係数に一致します。 $$R^2={Corr(y,\hat{y})}^2$$ # 決定係数 rsquared_1 = reg_v/total_v print("決定係数：{}".format(rsquared_1)) rsquared_2 = 1 - (res_v/total_v) print("決定係数：{}".format(rsquared_2)) rsquared_3 = r_yy_hat**2 # 重相関係数 print("決定係数：{}".format(rsquared_3)) 考察 いつもこういう式を見ると、なんで１だといい評価ってなるの？って気になっちゃうので、考えてみることにしてます。 ここは勝手な考察なので、信憑性は怪しいですが、自分はこれで納得してます。 回帰変動:$\sum(\hat{y_i}-\bar{y})^2$ 残差変動:$\sum\hat{\epsilon_i}^2$ 全変動:$\sum({y_i}-\bar{y})^2$ この３つの式を見てみると、全変動にはハットの部分がないので、全変動は推定モデルに関係しない式です。 yの配列が分かれば推定も何もしなくても全変動は出るということです。 対して、回帰変動と残差変動はハットの部分が含まれるので、推定しなきゃ始まらない部分です。 $\frac{回帰変動}{全変動}$について モデル関係なしに出てくる平均との差の２乗を合計した全変動 モデルを考慮して平均との差の２乗を合計した回帰変動が同じになる $（回帰変動）＝（全変動）$みたいな状況ができたらそれが理想な訳ですね。 この式の両辺を全変動で割れば、決定係数の評価の式になりますし、 左辺が1になるので理想状態で１なんだなぁ〜となんとなくわかるかなと。 $1-\frac{残差変動}{全変動}$について 残差変動はモデルが絡む部分ですが、$０$だと理想な訳です。 ということは （モデルに関係ない部分）$-$（残差変動の理想値:$0$） $=$ （モデルに関係ない部分） （全変動）$-$（残差変動) $=$ (全変動) が理想です。 これの両辺を全変動で割れば、やっぱり決定式の評価式の理想状態になって、それが１なんだなぁ〜となんとなく腑に落ちます。 本題に戻ります。 自由度調整済み決定係数 重回帰モデルで説明変数の個数が違うもの同士の精度を比較することができます。 $$\bar{R}^2=1-\frac{\frac{残差変動}{残差変動自由度}}{\frac{全変動}{全変動自由度}}$$ # 自由度調整済み決定係数 rsquared_adj = 1 - ((res_v/res_free) / (total_v/total_free)) print("自由度調整済み決定係数：{}".format(rsquared_adj)) 回帰係数の区間推定 ここまでで推定してきたわけですが、 $$\hat{\beta_1} = \frac{\sum(x_i-\bar{x})(y_i-\bar{y})}{\sum(x_i-\bar{x})^2}$$ これは点推定です。 ここから区間推定をしていきます。 最小二乗推定量の分布 $$\hat{\beta_1} 〜 N(\beta_1,\frac{\sigma^2}{\sum(x_i-\bar{x_i})^2})$$ この分布がどこからやってきたのか、あまり詳しいことは分かりません。分散の分母はなんでこんな形になってんでしょうか。 ここに書いてあるよ〜とか詳しい方いたら教えてください。 母平均の推定を思い出してみる。 平均$\mu$、分散$\sigma^2$の小標本母集団からの抽出($t$分布を使うやつ)で母平均の信頼区間を推定したときのことを思い出すと $\bar{X} 〜 N(\mu,\frac{\sigma^2}{n})$で標本平均$\bar{X}$、不偏分散を$\hat{\sigma}^2$とすると 不偏分散： $$\hat{\sigma}^2=\frac{1}{n-1}\sum(x_i-\bar{x})^2$$ 標準誤差 $$SE = \sqrt\frac{\hat{\sigma}^2}{n}$$ を用いて $$ t = \frac{\bar{X}-\mu}{SE}$$ このtが自由度$n-1$の$t$分布に従うというやつです。 $$ -t_{0.025}≦\frac{\bar{X}-\mu}{SE}≦t_{0.025}$$ $$ \bar{X}-t_{0.025}SE ≦ \mu ≦\bar{X}+t_{0.025}SE$$ この辺りから論理的に理解していないのでこの場合と同じノリで回帰係数も推定します。 回帰係数の信頼区間を推定する 回帰係数の場合は $$\hat{\beta_1} 〜 N(\beta_1,\frac{\sigma^2}{\sum(x_i-\bar{x_i})^2})$$ に従い、標本平均のところを$\hat{\beta_1}$、誤差の不偏分散を$\hat{\sigma}^2$とすると 誤差の不偏分散： $$\hat{\sigma}^2 =\frac{1}{n-k-1}\sum\hat{\epsilon_i}^2 =\frac{1}{n-k-1}\sum({y_i-\hat{y_i}})^2$$ 標準誤差： $$SE(\hat{\beta_1}) = \sqrt\frac{\hat{\sigma}^2}{\sum(x_i-\bar{x})^2}$$ を用いて $$ t = \frac{\hat{\beta_1}-\beta_1}{SE(\hat{\beta_1})}$$ このtが自由度$n-k-1$の$t$分布に従います。 $$ -t_{0.025}≦\frac{\hat{\beta_1}-\beta_1}{SE(\hat{\beta_1})}≦t_{0.025}$$ $$ \hat{\beta_1}-t_{0.025}SE(\hat{\beta_1}) ≦ \beta_1 ≦\hat{\beta_1}+t_{0.025}SE(\hat{\beta_1}) $$ これで回帰係数の推定、$\beta_1$の下限と上限が求まりました。 回帰係数の検定 帰無仮説：$H_0$：$\hat{\beta_1}=0$ 対立仮説：$H_1$：$\hat{\beta_1}≠0$ xがyに関係ないんじゃないの？っていう検定です。 $t$の式で$\hat{\beta_1}=0$を考えればいいので、 帰無仮説を想定した上での$t$値は　回帰係数の推定値 $÷$ 標準誤差 で出ます。 これも出力結果の読み取りで統計検定2級によく出てました。 # 誤差分散（不偏分散で置き換え） sigma_var = res_v/res_free print("誤差分散：{}".format(sigma_var)) # 回帰係数の標準誤差 beta1_se = np.sqrt(sigma_var/sum((x - x_mu)**2)) print("回帰係数の標準誤差:{}".format(beta1_se)) # 自由度に基づいたt分布のパーセント点 t_975 = stats.t.ppf(0.975,len(x)-x.reshape(-1,1).shape[1]-1) t_975 # 回帰係数の上限・下限 beta1_up = beta1 + t_975*beta1_se print("回帰係数の上限値:{}".format(beta1_up)) beta1_low = beta1 - t_975*beta1_se print("回帰係数の上限値:{}".format(beta1_low)) # 帰無仮説：beta1 = 0 とした場合のt値 t_value1 = beta1/beta1_se print("t_value:{}".format(t_value1)) プログラムを最後にまとめておきます。 プログラム import numpy as np from scipy import stats # ここだけ変えてね x = np.array([2,3,4,5]) y = np.array([6.4,7.9,7.1,8.8]) # 平均 x_mu = np.mean(x) print("xの平均:{}".format(x_mu)) y_mu = np.mean(y) print("yの平均:{}".format(y_mu)) # 回帰係数の推定に必要なもの x_x_mu = x - x_mu # print("x-x_mu:{}".format(x_x_mu)) y_y_mu = y - y_mu # print("y-y_mu:{}".format(y_y_mu)) T_xy = np.sum(x_x_mu * y_y_mu) T_xx = np.sum(x_x_mu**2) T_yy = np.sum(y_y_mu**2) # 回帰係数の推定 beta1 = T_xy/T_xx beta0 = y_mu - beta1*x_mu print("beta1:{}".format(beta1)) print("beta0:{}".format(beta0)) print("回帰直線：y={}+{}x".format(beta0,beta1)) # 相関係数 r_xy_1 = T_xy / np.sqrt(T_xx*T_yy) print("相関係数：{}".format(r_xy_1)) r_xy_2 = np.sqrt(T_xx/T_yy)*beta1 print("相関係数：{}".format(r_xy_2)) r_xy_3 = np.cov(x,y,ddof = 1)[0,1] / np.sqrt((np.var(x, ddof = 1)*np.var(y, ddof =1))) print("相関係数：{}".format(r_xy_3)) # 相関係数行列 print("相関係数行列") print(np.corrcoef(x,y)) # 予測値 y_hat = beta0 + beta1*x print("予測値:{}".format(y_hat)) # 残差 residual = y - y_hat residual print("残差:{}".format(residual)) # 重相関係数 r_yy_hat = np.corrcoef(y,y_hat)[0,1] print("重相関係数:{}".format(r_yy_hat)) # 回帰変動 reg_v = np.sum((y_hat-y_mu)**2) print("回帰変動：{}".format(reg_v)) # 残差変動 res_v = np.sum(residual**2) print("残差変動：{}".format(res_v)) # 全変動 total_v = np.sum((y-y_mu)**2) print("全変動：{}".format(total_v)) # 回帰変動の自由度 reg_free = x.reshape(-1,1).shape[1] print("回帰変動の自由度：{}".format(reg_free)) # 全変動の自由度 total_free = len(x)-1 print("全変動の自由度：{}".format(total_free)) # 残差変動の自由度 res_free = total_free - reg_free print("残差変動の自由度：{}".format(res_free)) # 決定係数 rsquared_1 = reg_v/total_v print("決定係数：{}".format(rsquared_1)) rsquared_2 = 1 - (res_v/total_v) print("決定係数：{}".format(rsquared_2)) rsquared_3 = r_yy_hat**2 print("決定係数：{}".format(rsquared_3)) # 自由度調整済み決定係数 rsquared_adj = 1 - ((res_v/res_free) / (total_v/total_free)) print("自由度調整済み決定係数：{}".format(rsquared_adj)) # 誤差分散（不偏分散で置き換え） sigma_var = res_v/res_free print("誤差分散：{}".format(sigma_var)) # 回帰係数の標準誤差 beta1_se = np.sqrt(sigma_var/sum((x - x_mu)**2)) print("回帰係数の標準誤差:{}".format(beta1_se)) # 自由度に基づいたt分布のパーセント点 t_975 = stats.t.ppf(0.975,len(x)-x.reshape(-1,1).shape[1]-1) t_975 # 回帰係数の上限・下限 beta1_up = beta1 + t_975*beta1_se print("回帰係数の上限値:{}".format(beta1_up)) beta1_low = beta1 - t_975*beta1_se print("回帰係数の上限値:{}".format(beta1_low)) # 帰無仮説：beta1 = 0 とした場合のt値 t_value1 = beta1/beta1_se print("t_value:{}".format(t_value1)) stats.model 一生懸命プログラムを作りましたが、実は stats.model というライブラリで全て数行で解決します。 import pandas as pd import numpy as np import scipy as sp from scipy import stats import statsmodels.formula.api as smf import statsmodels.api as sm df = pd.DataFrame({ "x" : [2,3,4,5], "y" : [6.4,7.9,7.1,8.8], }) # statsmodel # 最小二乗法(OLS)による回帰係数の推定 lm_model = smf.ols(formula = "y ~ x", data = df).fit() lm_model.summary() R-squaredが決定係数 Adj.R-squaredが自由度調整済み決定係数 Interceptの行に切片$\beta_0$の推定値・標準誤差・検定の際のt値・p値・下限・上限が並びます。 その下が$\beta_1$のそれぞれです。 重回帰の場合 ちなみにstatsmodelを使うと重回帰の際の回帰係数の推定もすぐ終わります。 $$y_i = \beta_0 + \beta_1x_{i1} + \beta_2x_{i2} + \epsilon_i$$ 先程の単回帰モデルよりも一つ説明変数を増やして、説明変数はx1,x2とします。 プログラムではx2の配列ができているところと、formulaの後の書き方が 目的変数　〜　説明変数　となっている部分に説明変数が増えたのでx1+x2となってます。 import pandas as pd import numpy as np import scipy as sp from scipy import stats import statsmodels.formula.api as smf import statsmodels.api as sm df = pd.DataFrame({ "x1" : [2,3,4,5], "x2" : [1,1,0,1], "y" : [6.4,7.9,7.1,8.8], }) # statsmodel # 最小二乗法(OLS)による回帰係数の推定 lm_model = smf.ols(formula = "y ~ x1+x2", data = df).fit() lm_model.summary() satsmodelはRライクな感じで動かせて、Rでも同じ操作感で使えます。

こんにちは，株式会社Nospareリサーチャー・千葉大学の小林です．本記事ではAn Introduction to sequential Monte Carloという逐次モンテカルロ法（sequential Monte Carlo, SMC）に関する本 に付随しているparticlesというPythonモジュールを紹介したいと思います．SMCはたくさんの変数（パーティクル）を発生させ，それらのリサンプリング，値の更新，重み付けを繰り返すことで逐次的にフィルタリング密度の近似などを行うモンテカルロ法の一種で，状態空間モデルなどに対して広く適用することができます．線形ガウス状態空間モデルの場合にはカルマンフィルターが適用できますが，非ガウス非線形モデルの場合にはシミュレーションに基づくSMCが有用となります． 状態空間モデル まず状態空間モデルについて簡単に説明します．状態空間モデルは観測されるデータ$Y_t, (t=0,\dots,T)$と観測されない状態変数$X_t, (t=0,\dots,T)$からなる時系列モデルです．$X_t$はマルコフ過程に従い，$Y_t$は$X_t$を所与として条件付き独立です．状態変数の初期分布の密度関数を$p_0^\theta(x_0)$，$X_t$の条件付き密度を$p_t^\theta(x_t|x_{t-1})$，$Y_t$の条件付き密度を$f_t^\theta(y_t|x_t)$とすると，同時密度は p^\theta_T(y_{0:T},x_{0:T})=p_0^\theta(x_0)f_t^\theta(y_0|x_0)\prod_{t=1}^Tf_t^\theta(y_t|x_t)p^\theta(x_t|x_{t-1}) で与えられます．$\theta$はモデルのパラメータです．状態変数の初期値$X_0$が$p_0^\theta$から生成され，各$X_t$の値は１期前の値$X_{t-1}=x_{t-1}$に依存して$p_t^\theta$によって生成され，$t$期の観測値$Y_t$はその期の$X_t=x_t$の値を所与として$f_t^\theta$から生成される，という構造です．以下の模式図を見ると，背後で観測されない$X_t$の過程が走っていて各$X_t=x_t$の値のもとで$Y_t$の値が生成される，ということがわかりやすいかと思います． 状態空間モデルに基づいた分析においては，（$\theta$の値を固定して）次に挙げる点について関心があります． フィルタリング（filtering）：$Y_{0:t}=y_{0:t}$を所与としたときの$X_t$の条件付き分布 状態の予測：$Y_{0:t}=y_{0:t}$を所与としたときの$X_{t+1}$の条件付き分布 スムージング（smoothing）：$Y_{0:T}$を所与としたときの$X_{t}$の条件付き分布 データの予測：$Y_{0:t}=y_{0:t}$を所与としたときの$Y_{t+1}$の条件付き分布 尤度の計算：$L_T(\theta)=\int p^\theta_T(y_{0:T},x_{0:T})dx_{0:T}$ フィルタリング密度，状態の予測密度，スムージング密度は積分の漸化式で表されます．状態空間モデルモデルのパラメータ$\theta$を推定したい場合には，最尤法の場合では上記の$L_T$を最大化する必要がありますが，尤度を計算するためにはこの$T+1$次元の積分を解く必要があります．ガウス線形モデルの場合にはカルマンフィルターによって解析的にこれらの分析を行うことができますが，一般の状態空間モデルについては解析的に解くことはできません．SMCではこれらの数量をパーティクルに基づいた近似で計算します． 例：SVモデル 本記事では状態空間モデルの例として，ファイナンス分野などでよく利用される確率的ボラティリティ変動モデル（stochastic volatility model, SVモデル）を取り上げます．まずSVモデルは次のように表記されます． Y_t|X_t=x_t \sim N(0,\exp(x_t)), \quad X_0\sim N(\mu, \sigma^2/(1-\rho^2)),\quad X_t|X_{t-1}=x_{t-1}\sim N(\mu+\rho(x_{t-1}-\mu),\sigma^2). ファイナンスでは$Y_t$は資産の（例えば日次の）対数収益などを表し，平均がゼロ，分散が$\exp(X_t)$の正規分布に従います．$X_t$はいわゆる対数ボラティリティ（分散の対数）で，平均$\mu$のAR(1)過程に従っており，初期分布はこのAR過程に基づいたものになります．$Y_t$と$X_t$の条件付き分布は両方とも正規分布ですが，$Y_t$の条件付き分散に指数関数が入っているので非線形モデルになっています． パーティクルフィルター パーティクルフィルターは状態変数のシミュレーションを行うことで関心のある数量を逐次的に近似します．冒頭でも言及したように，パーティクルフィルターはたくさんの状態変数をシミュレートし，それらに対してリサンプリング，値の更新，重みの計算を繰り返し，フィルタリング分布を代表するようなパーティクルが残るようにします．ここではパーティクルフィルター（bootstrapフィルター）の擬似コードだけを紹介しておきます．パーティクルを$X_t^n$，対応するパーティクルの重みを$w_t^n$，正規化された重みを$W_t^n$と表記します（$ n=1,\dots,N, t=1,\dots,T$）． 初期化（$n=1,\dots,N$）：パーティクル$X_0^n$を初期分布から生成，ウェイトを作成$w_0^n=f_0(y_0|X_0^n)$，正規化する$W_0^n=w_0^n/\sum_{m=1}^Nw_0^m$． $t=1,\dots,T$について次のステップを繰り返す： もし$ESS(W_{t-1}^{1:N})<ESS_{\min}$であれば，パーティクルのインデックス（$A_t^{1:n}$）を重み$W_{t-1}^{1:N}$に基づいてリサンプリングし，$\hat{w}^n_t=1$とする．そうでなければパーティクルのセットを引き続き使う（$A_t^n=n$，$\hat{w}^n_t=w^n_{t-1}$） パーティクル$X_t^n$の値を$p^\theta_t(x_t|x_{t-1}^{A_t^n})$の分布から生成して更新する． 重みを計算する：$w_t^n=\hat{w}_t^n f_t^\theta(y_t|X_t^n)$ 重みを正規化する：$W_t^n=w_t^n/\sum_{m=1}^Nw_t^m$ $ESS$は有効サンプルサイズで，これが小さくなるとごく少数のパーティクルが大きな重みを持っている状態になるので，あらかじめ設定した$ESS_{\min}$を下回った場合にパーティクルのリサンプリングを行います． パーティクルフィルターのアウトプットを使って，$t$期のフィルタリング分布は経験分布$\sum_{n=1}^N W_t^n\delta_{X_t^n}$で近似されます．また，フィルタリング分布についての$\varphi(X_t)$の期待値は$\sum_{n=1}^NW_t^n\varphi(X_t^n)$，$t$期のデータの周辺条件付き分布$L_t=f_t(y_t|y_{0:t-1})$は$\sum_{n=1}^Nw_t^n/\sum_{n=1}^N w_{t-1}^n$で近似できます（リサンプリングした場合には$N^{-1}\sum_{n=1}^Nw_t^n$）．またスムージングもパーティクルフィルターのアウトプットを使って近似することができます． Bootstrapフィルターではパーティクルの値を更新するときに状態方程式の分布から生成していますが，より効率的にするために提案分布を使うこともできます（guidedパーティクルフィルター）． particlesパッケージ 前置きが長くなりましたが，ここからパッケージの使用について紹介します． 導入 パッケージのインストールはpipなどでできます．またパッケージのドキュメンテーションはこのページにあります．このパッケージについて詳しく知りたい方はぜひ読んでみてください．githubのレポジトリには，本の中で紹介されている数値例に関するコードやデータセットが含まれています． モデルの定義と人工データ生成 パッケージの準備ができたら，標準的なライブラリと一緒にパッケージ内のモジュールをインポートします．ここではstate_space_models，distributions，Momentsというそれぞれ状態空間モデルと確率分布を扱うモジュールとフィルタリング分布のモーメントを計算するモジュールをインポートしています． import numpy as np from matplotlib import pyplot as plt import particles from particles import state_space_models as ssm # 状態空間モデルを扱うモジュール from particles import distributions as dists # 確率分布を扱うモジュール from particles.collectors import Moments # フィルタリング分布のモーメントを計算するモジュール 次にモデルのクラスを定義します．ここでは上記のSVモデルを定義しています． class StochVol(ssm.StateSpaceModel): default_params = {'mu': -1.0, 'rho': 0.9, 'sigma': 1.0 } # X_0の初期分布 def PX0(self): sig0 = self.sigma/np.sqrt(1.0-self.rho**2) return dists.Normal(loc=self.mu, scale=sig0) # X_t|X_{t-1}の分布 def PX(self, t, xp): return dists.Normal(loc=self.mu + self.rho * (xp-self.mu), scale=self.sigma) # Y_t|X_tの分布 def PY(self, t, xp, x): return dists.Normal(scale=np.exp(0.5*x)) このクラスの属性はモデルのパラメータの値（mu, rho, sigma）で，このクラスは３つのメソッド（PX0, PX, PY）が含まれています．３つのメソッドはそれぞれ$p_0^\theta(x_0)$，$p_t^\theta(x_t|x_{t-1})$，$f_t^\theta(y_t|x_t)$に対応する確率分布の記述となります．PX内のxpは$x_{t-1}$に対応し，PY内のxは$x_t$に対応します．ここで利用可能な確率分布はdistributionsモジュールに関するドキュメンテーションに記載されています．このモジュールでは 対数密度の評価（logpdfメソッド） 乱数の発生（rvsメソッド） 分位点の計算（ppfメソッド） ができます． クラスを定義したら，パラメータの値を自分で設定してインスタンス化して，データをモデルから生成します．データの生成にはインスタンスに対してsimulateメソッドを適用します． my_sv_model = StochVol(mu=-11.0, rho=0.9, sigma=0.5) #インスタンス化 np.random.seed(seed=31) T = 200 # データ数 x, y = my_sv_model.simulate(T) # モデルからデータを生成 plt.plot(y, label='data') plt.legend() 以下が生成したデータ$Y_t$のプロットになります． フィルタリング ではパーティクルフィルターを適用します． N = 200 # パーティクルの数 fk_my_sv = ssm.Bootstrap(ssm=my_sv_model, data=y) # 状態空間モデルをFeynmanbootstrapフィルターを適用する alg = particles.SMC(fk=fk_my_sv, N=N, collect=[Moments()]) # SMCオブジェクトを作成 alg.run() # アルゴリズムをT期まで実行 まず定義した状態空間モデルからFeynman–Kacモデルに基づいたFeynmanKacオブジェクト（ssm.Bootstrapはこれのサブクラス）を作成します．この本ではより一般的なFeynman–Kacモデルに対するSMCを紹介しており，particlesパッケージでもFeynman–Kacモデルに対する適用がなされます．引数ssmには先程のSVモデルのインスタンスを，dataには観測データyを指定します．particles.SMCの引数fkには作成したFeynmanKacオブジェクト，Nにはパーティクル数を指定しています．一般的には alg = particles.SMC(fk=fk_guided, N=100, qmc=False, resampling='systematic', ESSrmin=0.5, store_history=False, verbose=False, collect=None) となっており，疑似モンテカルロの使用（qmc），リサンプリングの条件（ESSrmin)，経過の保存（store_history，スムージングで使用），モーメントの計算（collect）など設定できます． algオブジェクトはに対してrunメソッドを使うと，$T$期まで一気にフィルタリングを行います（ stepメソッドでは１期づつフィルタリングを進めていきます，step自体イテレータ）．このとき，algオブジェクトの属性は alg.t：次のイタレーションのインデックス alg.X：パーティクル$X_t^{1:N}$の値 alg.W： パーティクルの重み$W_t^{1:N}$ alg.Xp： １期前のパーティクルの値$X_{t-1}^{1:N}$ alg.A： パーティクルの祖先，A[3]=12は$X_t^3$の親を意味， alg.summaries：ESS（ESSs），リサンプリングのフラグ（rs_flags），データの対数周辺条件付き分布（logLts）の履歴 などがあります．この属性にあるとおり，runでデータを全部使ってフィルタリングを行うと，store_historyがFalseの場合にはパーティクルの情報は最終時点$T$の状態しか残りません．よってcollectにモーメントを計算する関数を指定してモーメントの履歴を保存します（また，いちいち履歴をさかのぼってモーメントを計算するより便利です）．上の適用例ではMomentsの引数に何も指定していないのでフィルタリング分布の平均と分散を計算します．モーメントの値はalg.summaries.momentsに辞書型で保存されます． フィルタリングの結果（フィルタリング分布の平均）を以下に図示します． plt.plot(x, label='true state', color='k', ls=':') plt.plot([m['mean'] for m in alg.summaries.moments], label='filtered mean') plt.legend() $\log L_t$のプロットも作成してみます． plt.plot(alg.summaries.logLts, label='log likelihood') plt.legend() 最初のフィルタリングの結果の図ですが，図としてちょびっと寂しい気がするので，フィルタリング分布の2.5%と97.5%の区間も追加してみたいと思います．自分で定義した関数をcollectに指定することも可能です．分位点を計算するnp.percentileは重みを入れることができないので，このページの関数を参考に，フィルタリング分布の平均，2.5%，97.5%点を計算する関数を作成しました． def alternative_moments(W, X): # https://stackoverflow.com/questions/21844024/weighted-percentile-using-numpy X = np.array(X) quantiles = np.array([0.025, 0.975]) if W is None: W = np.ones(len(X)) W = np.array(W) assert np.all(quantiles >= 0) and np.all(quantiles <= 1), \ 'quantiles should be in [0, 1]' sorter = np.argsort(X) X = X[sorter] W = W[sorter] weighted_quantiles = np.cumsum(W) - 0.5 * W weighted_quantiles /= np.sum(W) res = np.interp(quantiles, weighted_quantiles, X) return {'mean': np.average(X, weights=W), '2.5%': res[0], '97.5%': res[1]} この関数を以下のように指定することで各期の平均と分位点を保存することができます． # mom_funcにalternative_momentsを指定 alg = particles.SMC(fk=fk_my_sv, N=N, collect=[Moments(mom_func=alternative_moments)], store_history=True) alg.run() plt.plot(x, label='true state', color='k', ls=':') plt.plot([m['mean'] for m in alg.summaries.moments], label='filtered mean') q1 = np.array([m['2.5%'] for m in alg.summaries.moments]) q2 = np.array([m['97.5%'] for m in alg.summaries.moments]) qq = np.concatenate([q1,q2[::-1]]) s = np.arange(0,T) tt = np.concatenate([s,s[::-1]]) plt.fill(tt, qq, color='b', alpha=0.2) plt.legend() ちなみに以下でスムージングを行うのでstore_historyをTrueにしています．このコードでは以下の図をプロットします．陰の部分がフィルタリング分布の2.5%から97.5%の区間になります． スムージング 下記のコードは，パーティクルフィルターの結果（フィルタリングの際にstore_history=Trueとすること）を使って後ろ向きにパーティクルのトラジェクトリーを発生させていくforward filtering backward sampling（FFBS）の適用例です．フィルタリングの履歴alg.histにbackward_samplingを適用します．引数には M：トラジェクトリーの数（スムージングでのモンテカルロサンプルの数） linear_cost：Trueにすることで計算量$O(T(M+N))$の棄却FFBSを適用（棄却サンプリングにおけるバウンドをモデル内で指定する必要あり），Falseで通常のFFBSを実行（計算量は$O(TMN)$） 例えば，M=5の場合のトラジェクトリーは以下のようになります． #トラジェクトリーの数は5，通常のFFBSを使用 trj = np.array(alg.hist.backward_sampling(M=5, linear_cost=False)) plt.plot(trj) 以下がM=Nの場合の結果です．FFBSの各トラジェクトリーは独立に生成されるので，スムージング分布の平均は単純に各時点における標本平均を計算しています．おまけとして2.5%と97.5%の区間もプロットしています． trj = np.array(alg.hist.backward_sampling(N, linear_cost=False)) plt.plot(x, label='true state', color='k', ls=':') plt.plot(np.average(trj, axis=1), label='smoothed mean') q1 = np.percentile(trj, 2.5, axis=1) q2 = np.percentile(trj, 97.5, axis=1) qq = np.concatenate([q1,q2[::-1]]) s = np.arange(0,T) tt = np.concatenate([s,s[::-1]]) plt.fill(tt, qq, color='b', alpha=0.2) plt.legend() 最後に このパッケージは状態空間モデルを定義すればあとはフィルタリングやモデルのベイズ推定を簡単に行うことができます．本記事の実装例に対するパーティクルフィルターやスムージングは驚くほど速く動作します．状態空間モデルは時系列データ分析においてとても重要なので，このパッケージはとても有用なものであると思います．本自体もかなりしっかり書かれていてとても勉強になります．本記事ではパラメータの値を固定した状態でのフィルタリングやスムージングを行いましたが，現実の分析ではパラメータの値は推定します．次の記事では，このパッケージを使ったパラメータのベイズ推定を取り上げたいと思います． 米倉先生の記事でも書かれていたように，ビジネスにおいても時系列データはちゃんと時系列モデルに基づいた分析をするべきです．本記事で紹介した状態空間モデルの枠組みは，例えば売上高と広告費が毎日観測されていて，広告費にどれほど支出すると売上高がどの程度になるかを予測したい，といった場合に利用できます．このとき，時間とともに変化する広告費あるいはアドストックの効率性や一週間・一ヶ月周期で変化する要因などを観測されない状態変数として組み込むことが考えられます．また一週間にあるSKUに対する注文が何個分入るのかの分析・予測を行う際にも同じようにモデルを立てることができます．本記事で紹介したような既存のパッケージをベースとしたビジネスデータ活用分析支援・分析システム開発支援も行っておりますのでお問い合わせください． 一緒にお仕事しませんか！ 株式会社Nospareではベイズ統計学に限らず統計学の様々な分野を専門とする研究者が所属しており，新たな知見を日々追求しています．統計アドバイザリーやビジネスデータの分析につきましては弊社までお問い合わせください．インターンや正社員も随時募集しています！