※この記事は，線形演算的に快適なM1 native Python3.9環境が欲しい方が対象です。 はじめに 最近M1 Macを買って一応brew install pythonとpip install numpyでinstallまではしたけど，Apple Siliconなのにnumpyだけやたら遅いな？と感じていませんか…？それはきっとOpenBLASがlinkされていないからです。この機会に以下の方法でnumpyにOpenBLASをlinkして，最大300倍の高速化チャンス！を獲得しましょう。 前提 本記事では以下を仮定しています M1 Mac Arm版のHomebrew Arm版のPython3.9 インストール方法 numpy/scipyをOpenBLASにlinkしてinstallするコマンドです # numpy/scipyに必要 % brew install openblas gfortran % pip3 install cython pybind11 # おまじない % export OPENBLAS="$(brew --prefix openblas)/lib/" # build from source % pip3 install --no-binary :all: --no-use-pep517 numpy # おまけでscipyも（結構長いので注意） % pip3 install --no-binary :all: --no-use-pep517 scipy 検証 以下を実行してbenchmarkを計測しました ▼ https://gist.github.com/markus-beuckelmann/8bc25531b11158431a5b09a45abd6276 Name Python Platform BLAS / LAPACK 行列積 (4096x4096) [sec] ドット積 (1x524228) [ms] SVD (2048x1024) [sec] Cholesky分解 (2048x2048) [sec] 対角化 (2048x2048) [sec] Pure NumPy aarch64 (Homebrew) - 298.54 1.09 13.27 2.12 73.81 NumPy + OpenBLAS aarch64 (Homebrew) OpenBLAS 0.95 0.28 2.49 0.11 10.27 NumPy + Intel MKL intel (Miniconda) Intel MKL 2.53 0.08 0.96 0.22 8.16 NativeなArm版OpenBLASがエミュレートされているRosetta 2 + Intel MKLに行列積以外で負けているのは何故なんでしょう？AppleのRosetta 2がすごいのかそれもと数学ライブラリの開発をサボっているからなのかわかりません。そしてPure NumPyは遅すぎて途中で計測を辞めたくなりました。 検証ログ Pure NumPy (aarch64) Dotted two 4096x4096 matrices in 298.54 s. Dotted two vectors of length 524288 in 1.09 ms. SVD of a 2048x1024 matrix in 13.27 s. Cholesky decomposition of a 2048x2048 matrix in 2.12 s. Eigendecomposition of a 2048x2048 matrix in 73.81 s. This was obtained using the following Numpy configuration: blas_mkl_info: NOT AVAILABLE blis_info: NOT AVAILABLE openblas_info: NOT AVAILABLE atlas_3_10_blas_threads_info: NOT AVAILABLE atlas_3_10_blas_info: NOT AVAILABLE atlas_blas_threads_info: NOT AVAILABLE atlas_blas_info: NOT AVAILABLE blas_info: NOT AVAILABLE blas_src_info: NOT AVAILABLE blas_opt_info: NOT AVAILABLE lapack_mkl_info: NOT AVAILABLE openblas_lapack_info: NOT AVAILABLE openblas_clapack_info: NOT AVAILABLE flame_info: NOT AVAILABLE atlas_3_10_threads_info: NOT AVAILABLE atlas_3_10_info: NOT AVAILABLE atlas_threads_info: NOT AVAILABLE atlas_info: NOT AVAILABLE lapack_info: NOT AVAILABLE lapack_src_info: NOT AVAILABLE lapack_opt_info: NOT AVAILABLE numpy_linalg_lapack_lite: language = c define_macros = [('HAVE_BLAS_ILP64', None), ('BLAS_SYMBOL_SUFFIX', '64_')] NumPy w/ OpenBLAS (aarch64) qiita@m1 ~ % python numpy_benchmark.py Dotted two 4096x4096 matrices in 0.95 s. Dotted two vectors of length 524288 in 0.28 ms. SVD of a 2048x1024 matrix in 2.49 s. Cholesky decomposition of a 2048x2048 matrix in 0.11 s. Eigendecomposition of a 2048x2048 matrix in 10.27 s. This was obtained using the following Numpy configuration: blas_mkl_info: NOT AVAILABLE blis_info: NOT AVAILABLE openblas_info: libraries = ['openblas', 'openblas'] library_dirs = ['/opt/homebrew/opt/openblas/lib/'] language = c define_macros = [('HAVE_CBLAS', None)] runtime_library_dirs = ['/opt/homebrew/opt/openblas/lib/'] blas_opt_info: libraries = ['openblas', 'openblas'] library_dirs = ['/opt/homebrew/opt/openblas/lib/'] language = c define_macros = [('HAVE_CBLAS', None)] runtime_library_dirs = ['/opt/homebrew/opt/openblas/lib/'] lapack_mkl_info: NOT AVAILABLE openblas_lapack_info: libraries = ['openblas', 'openblas'] library_dirs = ['/opt/homebrew/opt/openblas/lib/'] language = c define_macros = [('HAVE_CBLAS', None)] runtime_library_dirs = ['/opt/homebrew/opt/openblas/lib/'] lapack_opt_info: libraries = ['openblas', 'openblas'] library_dirs = ['/opt/homebrew/opt/openblas/lib/'] language = c define_macros = [('HAVE_CBLAS', None)] runtime_library_dirs = ['/opt/homebrew/opt/openblas/lib/'] NumPy w/ Intel MKL (x86_64) Dotted two 4096x4096 matrices in 2.53 s. Dotted two vectors of length 524288 in 0.08 ms. SVD of a 2048x1024 matrix in 0.96 s. Cholesky decomposition of a 2048x2048 matrix in 0.22 s. Eigendecomposition of a 2048x2048 matrix in 8.16 s. This was obtained using the following Numpy configuration: blas_mkl_info: libraries = ['mkl_rt', 'pthread'] library_dirs = ['/Users/qiita/miniconda3/lib'] define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)] include_dirs = ['/Users/user/miniconda3/include'] blas_opt_info: libraries = ['mkl_rt', 'pthread'] library_dirs = ['/Users/qiita/miniconda3/lib'] define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)] include_dirs = ['/Users/user/miniconda3/include'] lapack_mkl_info: libraries = ['mkl_rt', 'pthread'] library_dirs = ['/Users/qiita/miniconda3/lib'] define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)] include_dirs = ['/Users/user/miniconda3/include'] lapack_opt_info: libraries = ['mkl_rt', 'pthread'] library_dirs = ['/Users/qiita/miniconda3/lib'] define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)] include_dirs = ['/Users/qiita/miniconda3/include']

このプログラムではこんなことします datetimeモジュール，calendarモジュールを使用 第N曜日を文字列で受け取り，日にちに変換 既に過ぎていた場合は，翌月第N何曜日を計算 サンプル.py 元の文字列は"EveryDay"や"FirstMonday,EverySunday"などです．これを第-7回の関数separateStrで分割し入力しています． sample.py import datetime import calendar def str2date(str_list): date_list = [] for s in str_list: frequency = s[0] weekday = s[1] delta = weekday_list.index(weekday) - today.weekday() if frequency == "Every": if weekday == "Day": date = today else: if delta < 0: delta += 7 date = today + datetime.timedelta(days=delta) else: nth = frequency_list.index(frequency) + 1 weekday_first, _ = calendar.monthrange(year=today.year, month=today.month) day = 7 * (nth-1) + (weekday_list.index(weekday) - weekday_first) % 7 + 1 if day >= today.day: delta = day - today.day date = today + datetime.timedelta(days=delta) else: weekday_first, _ = calendar.monthrange(year=today.year, month=today.month+1) day = 7 * (nth-1) + (weekday_list.index(weekday) - weekday_first) % 7 + 1 date = datetime.date(year=today.year, month=today.month+1, day=1) + datetime.timedelta(days=day-1) date_list.append(date) return date_list frequency_list = ['First', 'Second', 'Third', 'Fourth'] weekday_list = ['Monday', 'Tuesday', 'Wednesday', 'Thursday', 'Friday', 'Saturday', 'Sunday'] today = datetime.date.today() 第N曜日の計算にはnkmkさんのサイトを参考にさせていただきました．

目次 ルーブックキューブをモーターでそろえる ただただ　自分メモ代わりに 　ちょっと　でも誰かのためになったり いや、やっぱり　自分がやったこと　忘れても　なんとなく　理解できるよう 環境 python3.9.1 Opencv 4.5.1 Arduino uno デイジーチェーン Visual Studio Code Windows10 ３Dプリンター(Adventure3, 材料PLA) 使用道具 秋月ステッパーモータ (バイポーラ　ステッピングモーター　ＳＭ－４２ＢＹＧ０１１) 秋月L6470モータドライバー (Ｌ６４７０使用　ステッピングモータードライブキット) カップリング ルービックキューブをモーターをつなぐ　ジョイント 直流安定化電源 30V-10A Arduino uno Arduino uno 互換マイコンボード お世話になった　サイト L6470とモータのつなぎ方 工学部生の生プレスさん arduino、ステッパーモータ、L6470モータドライバのつなぎかたなど arduinoプログラム 北の国から電子工作（仮）さん 3Dプリンター キューブをまわすボディ キューブ回転記号 キューブの色を取得コード（忘れた、申し訳ない) Python,Opencvで画像処理 カメラはスマートフォン（galaxy s9) カメラアプリ　　iriuncam NO.1 python opencv ウェブカメラで　リアルタイムで色（BGR)を取得

こちらは　NO.1 もくじはここ カメラ起動 cv2.VideoCapture(0) ( )の中身の数字を変えて　起動するカメラを指定する。 0が　スタート　複数カメラ接続がある場合　１，２など試す。 if cv2.waitKey(1) & 0xFF == ord('q'): キーボードのQボタンで　　ブレイク　終了します。 （ウィンドウが閉じる） サンプルコード（カメラ） import cv2 capture = cv2.VideoCapture(2) while True: ret, frame = capture.read() cv2.imshow('camra',frame) if cv2.waitKey(1) & 0xFF == ord('q'): break capture.release() cv2.destroyAllWindows() サンプルコード（カメラ）を実行 ウィンドウをクリックして　'Q'ボタンを押せば終了できます。 ＊右上の　＿　□　ｘ　ボタンは使えません 四角い枠で囲む opencv cv2.rectangleを使用して　四角い枠を表示する 左上が(X0, Y0)　 左右が　X方向: 上下が　Y方向: サンプルコード（四角い枠を表示） import cv2 capture = cv2.VideoCapture(2) #capture.set(cv2.CAP_PROP_FRAME_WIDTH, 1280) # 横幅 #capture.set(cv2.CAP_PROP_FRAME_HEIGHT, 720) # 縦幅 x_s = 100 y_s = 100 x_e = 400 y_e = 400 #Y縦　X横 while True: ret, frame = capture.read() colorBGR = cv2.rectangle(frame, 　　　　　　　　　(x_s, y_s), (x_e, y_e), (255, 255, 255), 2) cv2.imshow('camra',frame) if cv2.waitKey(1) & 0xFF == ord('q'): break capture.release() cv2.destroyAllWindows() #追加したコード colorBGR = cv2.rectangle(frame, 　　　　　　　　　(x_s, y_s), (x_e, y_e), (255, 255, 255), 2) 今回は使っていませんが　自分で　 ウィンドウの大きさを指定してあげることができる #capture.set(cv2.CAP_PROP_FRAME_WIDTH, 1280) # 横幅 #capture.set(cv2.CAP_PROP_FRAME_HEIGHT, 730) # 縦幅 ret, frame = capture.read() この frameを使う -----------コーーーーード---------------------- #Y縦軸　X横軸 x_s = 100 y_s = 100 x_e = 400 y_e = 400 colorBGR = cv2.rectangle(frame, (x_s, y_s), (x_e, y_e), (255, 255, 255), 2) ------------------------------------------------- ------以下説明------- ””” cv2.rectangleは四角を表示する (frame 最初のフレームは　上の赤文字のやつ (x_s, y_s)　開始地点を指定します (x_e, y_e) 終了地点を指定します (255, 255, 255)　B,G,Rの順番で指定する　今回は白色 (0, 0, 255) なら　赤色 , 2)　最後の　２は　線の太さ　０にした場合　線ではなく　指定した範囲 全部塗りつぶします　 ０の場合 ■　このように　真っ白になる ””” ＃x_s, y_s　スタートの　ｓ ＃x_e, y_s　エンドの　e サンプルコード（四角い枠を表示）を実行 下の白い線が　したのコードで表示される colorBGR = cv2.rectangle(frame, (x_s, y_s), (x_e, y_e), (255, 255, 255), 2) このように　範囲を指定してあげればいいじゃん。 指定したところのBGRを取得 サンプルコード(BGR取得） import cv2 capture = cv2.VideoCapture(2) x_s = 280 y_s = 280 x_e = 320 y_e = 320 #Y縦　X横 while True: ret, frame = capture.read() colorBGR = cv2.rectangle(frame, (x_s, y_s), (x_e, y_e), (255, 255, 255), 2) pixel = frame[300, 300] print(pixel) cv2.imshow('camra',frame) if cv2.waitKey(1) & 0xFF == ord('q'): break capture.release() cv2.destroyAllWindows() #追加したコード pixel = frame[300, 300] print(pixel) frame[Y, X]で指定すれば　その座標のBGRを取得することができる 今回は fram[３００，３００] 四角い枠は　２８０～３２０でかこんでおり その　中心　３００の座標を取得しているので　四角い枠の中心部分の BGRを取得している。 サンプルコード(BGR取得）の実行 取得したBGRを　リアルタイムで　ターミナルに表示しています 先の話、、、この枠を　９個用意してしまえば ルービックキューブ色を９個取得できる。 ちょっと　おまけ カメラ接続確認 接続されてる　カメラが不明な場合ここで確認できる　（マイク・スピーカも可能） Mac Bookで起動すると　エラーがでる可能性ある。 パソコン側がカメラの使用拒否をしてるかもしれない。 カメラの使用許可をする必要がありる パソコン側で設定はできない（はず） abort trap6: PYthonが、、、、終了しました。 もし　このエラーがでたら こちらへどうぞ ターミナルから　VS CODEを起動するとMAC　からカメラの許可を　聞いてくるはず （自分はそうでした、ほかの方はしりません)

はじめに pythonでいろいろプログラムを実行していると、もっと高速に実行することはできないかなと思うことがあります。 特に性能のいいパソコンなどを持っている人だと、タスクマネージャーを開いてCPUの使用率を眺めていると1スレッドしか使われていないじゃん、と思うことがあると思います。 何とかしてすべてのcpuコアを使用して高速化できないかな、というときに少し調べたらこんな記事が出てきました。 私もこれを参考にいろいろ高速化に成功しましたが、いろいろ落とし穴などがあったため、並列化のやり方を1から解説しようと思います。 どんなコードなら高速化できるか なんでもかんでも高速化することができるわけではありません。 大前提として、forを使用したループがある場合、並列化によって高速化ができる可能性があります。 可能性といったのは、中には難しい場合もあるからで、のちに解説します。 実際どうやるの？ 例えば以下のようなプログラムを高速化しようとします。 for i in range(30): print(i) print("の2乗は") print(i*i) print() これを高速化する場合、二つのステップがあります。 1.ループ内の関数化 マルチスレッドにするのに伴い、ループ毎に行っている内容を関数として書き直します。 つまり今回は、このように書き直します。 def square(i): print(i) print("の2乗は") print(i*i) print() for i in range(30): square(i) 2.Poolの導入 次に、並列化を実装します。 並列化では、multiprocessingというpython標準ライブラリからPoolという関数を導入します。 from multiprocessing import Pool 次に、実行方法を変えます。 if __name__ == '__main__': main() こんなプログラムを見たことはないでしょうか。 最初見たときなんだこれ、と思いましたが普通に使う分には、if __name__ == '__main__':の中に実行したいプログラムを書くことでいつものように実行することができます。 よくある「おまじない」と解説されるような内容ですが、詳しく知りたい人はこんな記事を読んでください。 とにかく今回重要なのは、このif __name__ == '__main__':を使うということです。 コードは、以下のようになります。 わかりやすいように最初のコードも載せておきますね 最初のコード.py for i in range(30): print(i) print("の2乗は") print(i*i) print() 最終的なコード.py from multiprocessing import Pool def square(i): print(i) print("の2乗は") print(i*i) print() if __name__ == '__main__': p = Pool() p.map(square, range(30)) まず、 p = Pool() の部分ですが、Pool(6)などと指定すると、6スレッド使用するように指定することができます。何も書かないとそのパソコンのCPUの最大数を使ってくれるみたいです。 そして一番大切な部分が、これです。 p.map(square, range(30)) map関数の第一変数に実行したい関数を、第二変数に関数に代入したいものを指定します。 第二変数はリストだったりrangeだったり、forでのループ時に使用していたものそのままでいいです。 a = [1,3,5,7] p.map(square, a) などという書き方もできます。 ここで注意ですが、この並列化では先頭から順番に代入してくれません。詳しくはまた後程説明します。 注意 1.実行方法の変更 これは普段エディタ上でプログラムを動かしている人向けなのですが、基本的にエディタ上では並列化のプログラムは動きません。（たぶん） 私は普段anacondaについてきたspyderを使用しているのですが、そこでは動きません。 実行時はコマンドプロンプトから実行してください。 いちいちコマンドを入力するのが面倒臭いという人はこの記事を参考にバッチファイルを作るのがおすすめです。（Windows） 2.ループ内の実行順番 並列化では順番無視でどんどん実行していきます。 たとえば、以下のような数字を順番に表示させるプログラムがあるとすると、 from multiprocessing import Pool if __name__ == '__main__': p = Pool() p.map(print, range(3000)) こんな風に最初はいいのですが、途中から順番が狂ってきます。 0 1 2 3 : : 23 24 25 26 27 47 28 48 94 : : 2819 2866 2912 2956 2913 2957 2958 2959 2960 なので並列化させたい場合は、実行順番に意味がないプログラムを使用してください。 3.変数の数 この並列化では、並列化する関数に渡せる変数は一つのみです。 この一つというのはこういうことです。 ダメな例.py from multiprocessing import Pool def plus(a,b): print(a+b) if __name__ == '__main__': a = [1,3,5,7,9] b = [2,4,6,8,10] p = Pool() p.map(plus2, a, b) #ここでエラー mapはaとbを同時に渡せない　 関数に渡せるオブジェクトが1つということで、リストにして複数渡すということができます。 この場合の並列化のプログラムはこんな感じです。 いい例.py from multiprocessing import Pool def plus2(c): print(c[0] + c[1]) if __name__ == '__main__': c = [[2,5],[4,7],[8,9],[3,4],[5,11],[4,7],[2,5],[8,6]] p = Pool() p.map(plus2, c) 複数受け取る関数を書き直したくないという方は、これを参考にしてください。 並列化できないプログラム ここまで見ると、forループさえあればどんなプログラムでも並列化できるように思えます。 しかし、forループでなく並列化することによってできなくなることがあります。 その多くがfor特有のものだったりします。 以下は私が並列化を使用していて実際にできなかった例です。 これらのほかにも並列化できない内容があると思います。 1.breakとcontinue forを使っている中で、条件を満たしたらループ終了や、実行しないようにするためにbreakや continueを使用していることがあると思います。 しかし並列化するとループという概念ではなくなるので、以上のものが使われているとエラーが起きます。 ただし、並列化された関数の中にさらにforループがあり、その中のbreakや continueは関係ありません。 2.tryとexcept for内で、エラーが生じたときの処理としてtryとexceptを書くことがありますが、それがあっても動きません。 まとめ 並列化を使用していて詰まったことを中心に、深く理解しなくても高速化する方法を書きました。 特にopencvで画像を処理する際などに使用しているのですが、使用すると非常に高速化ができます。 ぜひ使用してみてください。 参考にしたサイト

前回 Python Win32プログラミング その1 ウインドウメッセージ ウインドウメッセージを利用してウインドウの操作ができます。例えば、メモ帳を開いているとき、そのテキストをPythonで取得することができます。 from ctypes import * FindWindowEx = windll.user32.FindWindowExW SendMessage = windll.user32.SendMessageW EnumWindows = windll.user32.EnumWindows EnumChildWindows = windll.user32.EnumChildWindows # WM_SETTEXT = 0x000c WM_GETTEXT = 0x000d WM_GETTEXTLENGTH = 0x000e WNDENUMPROC = WINFUNCTYPE(c_bool, POINTER(c_int), POINTER(c_int)) def EnumChildProc(hwnd, lParam): length = SendMessage(hwnd, WM_GETTEXTLENGTH, 0, 0) buff = create_unicode_buffer(length + 1) ret = SendMessage(hwnd, WM_GETTEXT, length + 1, buff) print(buff.value) return True def EnumWindowsProc(hwnd, lParam): global ghwnd length = GetWindowTextLength(hwnd) buff = create_unicode_buffer(length + 1) GetWindowText(hwnd, buff, length + 1) if "メモ帳" in buff.value: EnumChildWindows(hwnd, WNDENUMPROC(EnumChildProc), 0) return False return True EnumWindows(WNDENUMPROC(EnumWindowsProc), 0) WM_GETTEXT,WM_GETTEXTLENGTHがウインドウメッセージです。C言語だったらマクロが定義されていますが、Pythonでは定義値を調べてグローバル変数にするといった形になります。 WM_GETTEXT https://docs.microsoft.com/ja-jp/windows/win32/winmsg/wm-gettext SendMessageやPostMessageでメッセージを送信すれば、メモ帳内でそのメッセージの処理が呼び出されます。 (実際にテキストを保持しているのは子ウインドウなので、EnumChildWindowsで探しています) これは、C言語とPythonの間でプロセス間通信ができるということです。ウインドウによって対応するメッセージが異なるので、狙い通りの操作をするのはなかなか難しいですが、C言語のソフトに対してPythonでGUIをつくるといったこともできそうです。

この記事では電波輝線データの積分強度図の作り方について紹介します。 今回、例として国立天文台の FUGIN プロジェクトで得られた野辺山45m電波望遠鏡の CO 輝線のアーカイブデータを用います。データは http://jvo.nao.ac.jp/portal/nobeyama/ から FGN_03100+0000_2x2_12CO_v1.00_cube.fits をダウンロードしました (重いです)。 全速度積分 積分強度とは、いわばスペクトルの面積のことです。まずは前回同様 fits を読み込みます。 from astropy.io import fits import numpy as np from matplotlib import pyplot as plt import aplpy hdu = fits.open("~/your/fits/dir/FGN_03100+0000_2x2_12CO_v1.00_cube.fits")[0] data の shape を確認します。 print(hdu.data.shape) # (462, 848, 848) これは、Z (速度), Y (銀緯), X (銀経) のチャンネル数を表しています。つまり、848 x 848 のピクセルがあって、その一つ一つが462 チャンネルのスペクトルを持っているということです。 では適当な場所 (424, 424) のスペクトルを見てみましょう。python では以下のように書きます。 spe = hdu.data[: ,424, 424] plt.plot(range(len(spe)), spe) plt.show() ここで「:」は「全て」を意味します。また、python では 424 と書くと実際には 424+1=425 ピクセル目をとってくることに注意です。 縦軸は [K] ですが、横軸はチャンネルです。これを速度に変換します。色々方法はありますが、今回はシンプルに、0 チャンネル目 (一番最初のチャンネル) と 461 チャンネル目 (一番最後のチャンネル) の速度を計算して、線形補完することにします。そのために、astropy.wcs を使用します。 from astropy.wcs import WCS w = WCS("~/your/fits/dir/FGN_03100+0000_2x2_12CO_v1.00_cube.fits")[0] print(w) #WCS Keywords # #Number of WCS axes: 3 #CTYPE : 'GLON-SFL' 'GLAT-SFL' 'VRAD' #CRVAL : 31.0 0.0 -99675.0 #CRPIX : 424.5 424.5 1.0 #NAXIS : 848 848 462 チャンネル数を与えると座標が返ってくる wcs_pix2world 関数を使用します。 x_start, y_start, v_start = w.wcs_pix2world(0, 0, 0, 0) #最後の 0 は気にしなくて良い x_end, y_end, v_end = w.wcs_pix2world(hdu.data.shape[2]-1, hdu.data.shape[1]-1, hdu.data.shape[0]-1, 0) v_array = np.linspace(v_start, v_end, num=hdu.data.shape[0], endpoint=True) これを使って横軸を速度にします。v_array は m/s になっているので km/s にするために1000.0で割り算します。 spe = hdu.data[: ,424, 424] plt.plot(v_array/1000.0, spe) plt.xlabel("Velocity [km/s]") plt.ylabel("Tmb [K]") plt.show() これでスペクトルが描けました。 積分強度図とはスペクトルの面積のことでした。つまり、速度軸 (pythonでいう 0 軸目) 方向に足し合わせて速度幅をかければ完成です。 integ_map = np.nansum(hdu.data, axis=0)*hdu.header["CDELT3"]/1000.0 plt.imshow(integ_map, vmin=0, vmax=500) plt.show() 指定速度積分 指定した速度の範囲だけを積分したい場合があると思います。上の例では np.nansum に hdu.data、つまり data 全部を渡していました。これをスライスして渡してあげれば完了です。 スライスするチャンネルを計算するために、先程とは反対の wcs_world2pix を使用します。 例として、視線速度 25 km/sから 125 km/sを積分する場合を考えます。 v_start_wcs, v_end_wcs = 25.0*1000.0, 125.0*1000.0 x_tempo, y_tempo, v_tempo = w.wcs_pix2world(0, 0, 0, 0) #最後の 0 は気にしなくて良い x_start_ch, y_start_ch, v_start_ch = w.wcs_world2pix(x_tempo, y_tempo, v_start_wcs, 0) x_end_ch, y_end_ch, v_end_ch = w.wcs_world2pix(x_tempo, y_tempo, v_end_wcs, 0) v_start_ch, v_end_ch = int(round(float(v_start_ch))), int(round(float(v_end_ch))) この v_start_ch, v_end_ch を使って data をスライスします。 integ_map = np.nansum(hdu.data[v_start_ch:v_end_ch+1], axis=0)*hdu.header["CDELT3"]/1000.0 plt.imshow(integ_map, vmin=0, vmax=500) plt.show() 余計なチャンネルを除いたので、その分ノイズが減ってちょっと綺麗になりました。 積分強度図を fits として保存する 上でやった処理を効率化するために以下のような関数を定義します。 def v2ch(v, w): # v(km/s)をchに変える x_tempo, y_tempo, v_tempo = w.wcs_pix2world(0, 0, 0, 0) x_ch, y_ch, v_ch = w.wcs_world2pix(x_tempo, y_tempo, v*1000.0, 0) v_ch = int(round(float(v_ch), 0)) return v_ch def del_header_key(header, keys): # headerのkeyを消す import copy h = copy.deepcopy(header) for k in keys: try: del h[k] except: pass return h def make_new_hdu_integ(hdu, v_start_wcs, v_end_wcs, w): # 積分強度のhduを作る data = hdu.data header = hdu.header start_ch, end_ch = v2ch(v_start_wcs, w), v2ch(v_end_wcs, w) new_data = np.nansum(data[start_ch:end_ch+1], axis=0)*header["CDELT3"]/1000.0 header = del_header_key(header, ["CRVAL3", "CRPIX3", "CRVAL3", "CDELT3", "CUNIT3", "CTYPE3", "CROTA3", "NAXIS3"]) header["NAXIS"] = 2 new_hdu = fits.PrimaryHDU(new_data, header) return new_hdu 25 km/s から 125 km/s の積分強度の hdu を以下で作ります。 integ_hdu = make_new_hdu_integ(hdu, 25.0, 125.0, w) 前回と同様、aplpy で plot してみます。 fig = plt.figure(figsize=(8, 8)) f = aplpy.FITSFigure(integ_hdu, slices=[0], convention='wells', figure=fig) f.show_colorscale(vmin=1, vmax=500, stretch='log', cmap="jet", aspect="equal") plt.show() 名前をつけて fits で保存するには、以下を実行します。 integ_hdu.writeto("integ_25_125.fits", overwrite=True) 以上です。 リンク 目次

はじめに 今回はDjangoを使って友達リスト的なアプリを作成しようと思います。 プロジェクト名はTomodachi 一覧表示機能実装 friendslist/models.py # Friendモデルを作成 class Friend(models.Model): name = models.CharField(default="", max_length=20) furigana = models.CharField(default="", max_length=20) nickname = models.CharField(default="", max_length=20) sex = models.IntergerField(default=0) birthday = models.DateTimeField(auto_now=False, auto_now_add=False) birthplace = models.CharField(default="", max_length=30) hobby = models.CharField(default="", max_length=30) company = models.CharField(default="", max_length=30) created_at = models.DateTimeField(auto_now=False, auto_now_add=True) updated_at = models.DateTimeField(auto_now=True, auto_now_add=True) # マイグレーションする $ python manage.py makemigrations $ python manage.py migrate friendslist/admin.py # 管理画面からfriendを登録できるようにする from friendslist.models import User, Friend admin.site.register(Friend) friendslist/views.py # Friendモデルをインポートする from django.shortcuts import render from friendslist.models import Friend def index(request): friends = Friend.objects.all() context = { 'title': 'Tomodachi', 'friends': friends, } return render(request, 'friendslist/index.html', context) friendslist/index.html # 友達一覧表示する <div class="column col-md-8"> <h2>友達一覧</h2> <ul class="list-group"> {% for friend in friends %} <li class="list-group-item">{{ friend.name }}</li> {% endfor %} </ul> </div> friendslist/models.py # モデルの必須項目を任意項目にする class Friend(models.Model): name = models.CharField(default="", max_length=20) furigana = models.CharField(blank=True, null=True, max_length=20) nickname = models.CharField(blank=True, null=True, max_length=20) sex = models.IntegerField(default=0) birthday = models.DateField(blank=True, null=True, auto_now=False, auto_now_add=False) birthplace = models.CharField(blank=True, null=True, max_length=30) hobby = models.CharField(blank=True, null=True, max_length=30) company = models.CharField(blank=True, null=True, max_length=30) created_at = models.DateTimeField(auto_now_add=True) updated_at = models.DateTimeField(auto_now=True) 詳細表示機能実装 config/urls.py # 詳細表示のURLを記述する path('<slug:pk>/', views.friend), friendslist/views.py # 詳細表示の記述をする def friend(request, pk): friend = Friend.objects.get(pk=pk) context = { 'friend': friend } return render(request, 'friendslist/friend.html', context) templates/friendslist/friend.html # 詳細ページを作成する {% extends 'friendslist/base.html' %} {% block content %} <main class="container mb-5"> <div class="row"> <div class="col col-md-4"> <h2>カテゴリ一覧</h2> <ul class="list-group"> <li class="list-group-item active" aria-current="true">仕事</li> <li class="list-group-item">小学校</li> <li class="list-group-item">中学校</li> <li class="list-group-item">高校</li> <li class="list-group-item">大学</li> </ul> </div> <div class="column col-md-8"> <h2>友達一覧</h2> <ul class="list-group"> {% for friend in friends %} <li class="list-group-item">{{ friend.name }}</li> {% endfor %} </ul> </div> </div> </main> <div class="container mb-5"> <div class="row"> <div class="col"> <a class="btn btn-danger" href="#" role="button">友達一覧</a> </div> <div class="col"> <a class="btn btn-success" href="#" role="button">友達登録</a> </div> <div class="col"> <a class="btn btn-warning" href="#" role="button">カテゴリ一覧</a> </div> <div class="col"> <a class="btn btn-info" href="#" role="button">カテゴリ登録</a> </div> </div> </div> {% endblock %} 友達作成機能実装 config/urls.py # 作成機能のURLを記述する path('create', views.create), friendslist/forms.py # 作成フォームの設定をする from django import forms from friendslist.models import Friend class FriendForm(forms.ModelForm): class Meta: model = Friend fields = ( 'name', 'furigana', 'nickname', 'sex', 'birthday', 'birthplace', 'hobby', 'company', ) friendslist/views.py # 友達作成機能の記述をする def create(request): if request.method == 'POST': form = FriendForm(request.POST) if form.is_valid(): friend = form.save(commit=False) friend.save() return redirect('/') return render(request, 'friendslist/create.html') 友達更新機能実装 friend.html # 詳細ページを編集ページに変更する(作成ページのコピーからvalueに変更する) {% extends 'friendslist/base.html' %} {% block content %} <main class="container mb-5"> <form class="row" method="POST"> {% csrf_token %} <div class="col col-md-4"> <h2>友達作成</h2> <ul class="list-group"> <li class="list-group-item"> <label class="form-label">名前</label> <input type="text" class="form-control" id="id_name" name="name" value="{{ friend.name }}"> </li> <li class="list-group-item"> <label class="form-label">ふりがな</label> <input type="text" class="form-control" id="id_furigana" name="furigana" value="{{ friend.furigana }}"> </li> <li class="list-group-item"> <label class="form-label">ニックネーム</label> <input type="text" class="form-control" id="id_nickname" name="nickname" value="{{ friend.nickname }}"> </li> <li class="list-group-item"> <label class="form-label">性別</label> <select class="form-control" id="id_sex" name="sex"> <option value="0" {% if friend.sex == 0 %}selected{% endif %}>選択しない</option> <option value="1" {% if friend.sex == 1 %}selected{% endif %}>男性</option> <option value="2" {% if friend.sex == 2 %}selected{% endif %}>女性</option> </select> </li> <li class="list-group-item"> <label class="form-label">誕生日</label> <input type="date" class="form-control" id="id_birthday" name="birthday" value="{{ friend.birthday }}"> </li> <li class="list-group-item"> <label class="form-label">出身地</label> <input type="text" class="form-control" id="id_birthplace" name="birthplace" value="{{ friend.birthplace }}"> </li> <li class="list-group-item"> <label class="form-label">趣味</label> <input type="text" class="form-control" id="id_hobby" name="hobby" value="{{ friend.hobby }}"> </li> <li class="list-group-item"> <label class="form-label">就職先</label> <input type="text" class="form-control" id="id_company" name="company" value="{{ friend.company }}"> </li> <li class="list-group-item"> <button type="submit" class="btn btn-success float-end">更新</button> </li> </ul> </div> <div class="column col-md-8"> <h2>友達メモ</h2> <ul class="list-group"> <li class="list-group-item">メモ１</li> <li class="list-group-item">メモ２</li> <li class="list-group-item">メモ３</li> </ul> </div> </form> </main> <div class="container mb-5"> <div class="row"> <div class="col"> <a class="btn btn-danger" href="#" role="button">友達一覧</a> </div> <div class="col"> <a class="btn btn-success" href="#" role="button">友達登録</a> </div> <div class="col"> <a class="btn btn-warning" href="#" role="button">カテゴリ一覧</a> </div> <div class="col"> <a class="btn btn-info" href="#" role="button">カテゴリ登録</a> </div> </div> </div> {% endblock %} friendslist/views.py # 更新機能を追記する def friend(request, pk): if request.method == 'POST': friend = Friend.objects.get(pk=pk) form = FriendForm(request.POST, instance=friend) if form.is_valid(): form.save() friend = Friend.objects.get(pk=pk) context = { 'friend': friend } return render(request, 'friendslist/friend.html', context) 友達削除機能実装 config/urls.py # 削除機能のURLを記述する path('<slug:pk>/delete', views.delete, name="delete"), friendslist/views.py # 削除機能を実装する def delete(request, pk): try: friend = Friend.objects.get(pk=pk) except Friend.DoesNotExist: raise Http404 friend.delete() return redirect('/') index.html # 削除ボタンを実装する <li class="list-group-item">{{ friend.name }}<a href="{% url 'delete' friend.pk %}" class="btn btn-danger btn-sm float-end">削除</a></li> 各画面にリンクを貼る(リンクボタンのsnippetsを作成) snippets/linkbtns.html # リンクボタンのsnippetsを作成 <div class="container mb-5"> <div class="row"> <div class="col"> <a class="btn btn-danger" href="/" role="button">友達一覧</a> </div> <div class="col"> <a class="btn btn-success" href="create/" role="button">友達登録</a> </div> <div class="col"> <a class="btn btn-warning" href="#" role="button">カテゴリ一覧</a> </div> <div class="col"> <a class="btn btn-info" href="#" role="button">カテゴリ登録</a> </div> </div> </div> base.html # リンクボタンを反映させる {% include 'friendslist/snippets/linkbtns.html' %} snippets/header.html # ヘッダナビのリンクを貼る <li class="nav-item"> <a class="nav-link active" aria-current="page" href="/">友達一覧</a> </li> <li class="nav-item"> <a class="nav-link" href="create/">友達登録</a> </li> 詳細画面へのリンクボタンを作成する index.html # 詳細画面へのリンクボタンを追加する <li class="list-group-item"> {{ friend.name }} <div class="float-end"> <a href="{% url 'friend' friend.pk %}" class="btn btn-success btn-sm">詳細</a> <a href="{% url 'delete' friend.pk %}" class="btn btn-danger btn-sm">削除</a> </div> </li> さいごに 今回はCRUD機能を実装しました。 一覧表示 詳細表示 作成 更新 削除 次回は認証機能や カテゴリ、メモのリレーションをはろうと思います。

df.copy(deep = False) deep = True（デフォルト）：元のdfの行列・値を両方保持するためメモリを消費する。 deep = False：元のdfの行列のみ保持するため、deep = Trueよりもメモリを消費しない。 df.drop inplace = False（デフォルト）：対象のオブジェクトを返すのみ。 inplace = True：削除結果を反映させる。 train_df_shallow = train_df.copy(deep=False) train_df_shallow.drop("Sale",axis = 1,inplace= True) jupyter notebookにおいて、pandasで表示が省略されるのを防ぐ pd.set_option('display.max_columns', 160) # 160番目の列まで全て強制表示 pd.set_option('display.max_columns', None) # 全て強制表示 get_dummiesの使い方 get_dummiesは、文字列のまま使用できるので、いったんlabel列を文字列に戻す。 df['label'] = le.inverse_transform(df['label']) df size price label 1 100 cola 2 150 tea 3 200 coffee get_dummiesでlabel列をダミー変数化して、元のDataFrameに結合する。 df_dummy = pd.get_dummies(df['label']) df2 = pd.concat([df.drop(['label'],axis=1),df_dummy],axis=1) df2 size price coffee cola tea 1 100 0 1 0 2 150 0 0 1 3 200 1 0 0 これでダミー変数化は完成した。 多重共線性 回帰分析などにダミー変数を使用したい場合は、変数同士の相関が高いと多重共線性という問題が発生するため、one-hotエンコーディングの配列から列の1つを削除して使用しなければならない。例えばCoffee列を削除しても、cola列とtea列がともに0であればcoffeeであることがわかるので、情報としての欠落はない。 get_dummieのdrop_firstパラメータにTrueを渡すと最初の列を削除できる。 df_dummy = pd.get_dummies(df['label'],drop_first = True) df3 = pd.concat([df.drop(['label'],axis=1),df_dummy],axis=1) df3 size price cola tea 1 100 1 0 2 150 0 1 3 200 0 0 重回帰分析などでは、カテゴリーの数がN個ある場合、N個のダミー変数を作ると予測できないため、どれか1個の変数を除外します。scikit-learnの重回帰分析 LinearRegression() では、この問題を回避する設計にしているようで、N個のダミー変数を作成しても動作します。 enumerate関数（forループでインデックスを取得できる） Pythonのenumerate()関数を使うと、forループの中でリスト（配列）などのイテラブルオブジェクトの要素と同時にインデックス番号（カウント、順番）を取得できる。 enumerate関数を使ったforループ 引数にリストなどのイテラブルオブジェクトを指定する。 インデックス番号, 要素の順に取得できる。 0以外の数値から開始したい場合は、enumerate()関数の第二引数に任意の開始数値を指定する。 for i, name in enumerate(l, 1): print(i, name) # 1 Alice # 2 Bob # 3 Charlie カレントディレクトリの取得 import os os.getcwd() ファイル選択ダイアログの表示 root = tkinter.Tk() root.withdraw() fTyp = [("", "*.csv")] iDir = os.path.abspath(current_dir) tkinter.messagebox.showinfo('売上可視化プログラム', '処理したいCSVファイルを選択してください！') file = tkinter.filedialog.askopenfilename(filetypes = fTyp,initialdir = iDir) if file == "": print("CSV選択がキャンセルされました。終了します。") sys.exit(1) 文字判定 大文字かどうかを判定（isupper) str.isupper()を使います。 文字列中の英字 全てが大文字、かつ１文字以上の場合はTrue、そうでない場合はFalseを返します。 # 全て英字大文字 >>> 'ABCDE'.isupper() True # 全て英字小文字 >>> 'abcde'.isupper() False # 英字　大文字・小文字混在 >>> 'Abcde'.isupper() False # 英数字混在 >>> 'ABC12'.isupper() True # 英数字以外も混在 >>> 'Aa12\*'.isupper() False 小文字かどうかを判定（islower) str.islower()を使います。 文字列中の英字（上記参照）全てが小文字で、かつそれが1文字以上ある場合はTrue、そうでない場合はFalseを返します。 # 全て英字大文字 >>> 'ABCDE'.islower() False # 全て英字小文字 >>> 'abcde'.islower() True # 英字　大文字・小文字混在 >>> 'Abcde'.islower() False # 英数字混在 >>> 'abc12'.islower() True # 英数字も混在 >>> 'Aa12\*'.islower() False 英数字であることを判定（isalnum） str.isalnum()を使います。 文字列中の全ての文字が英数字で、かつ1文字以上ある場合はTrue、そうでない場合はFalseを返します。 # 全て英字大文字 >>> 'ABCDE'.isalnum() True # 全て英字小文字 >>> 'abcde'.isalnum() True # 英字　大文字・小文字混在 >>> 'Abcde'.isalnum() True # 英数字混在 >>> 'abc12'.isalnum() True # 英数字以外混在 >>> 'Aa12\*'.isalnum() False 英字かどうかを判定（isalpha） str.isalpha()を使います。 文字列中の全ての文字が英字（上記参照）で、かつ1文字以上ある場合はTrue、そうでない場合はFalseを返します。 # 全て英字大文字 >>> 'ABCDE'.isalpha() True # 全て英字小文字 >>> 'abced'.isalpha() True # 英字　大文字・小文字混在 >>> 'Abced'.isalpha() True # 英数字混在 >>> 'abc12'.isalpha() False # 英数字以外混在 >>> 'Aa12\*'.isalpha() False ASCII文字列かどうかを判定（isascii） str.isascii()を使います。 文字列が空、または全ての文字が ASCII（ U+0000 〜 U+007F） である場合はTrue、それ以外の場合はFalseを返します。 Python3.7から追加されたメソッドなのでバージョンに注意。 # 全て英字大文字 >>> 'ABCDE'.isascii() True # 全て英字小文字 >>> 'abced'.isascii() True # 英字　大文字・小文字混在 >>> 'abc12'.isascii() True # 英数字以外混在 >>> 'Aa12\*'.isascii() True # ASCIIコード以外 >>> 'あ'.isascii() False タイトルケース文字列かどうかを判定（istitle） str.istitle()を使います。 タイトルケースとは、文の最初以外の単語についても先頭文字を大文字にすることです。 このメソッドは以下の２つを満たす場合にTrueを返します。 大文字は、英字以外の文字の後のみに置かれる 小文字は、英字の後にだけ続く よって、以下の点に注意が必要です。 アポストロフィー（’）の後の文字が小文字だとタイトルケースと判定されません（例：Tom’s →　Tom’S) 翻訳業界などでは4文字未満の前置詞や接続詞（on, theなど）はこれを適用しないルールがあるようですが、このメソッドはこのルールが適用されません 例 # 文章の最初の文字のみ大文字 >>> "Tom's watch is on the table.".istitle() False # istitle()でTrueとなるタイトルケース >>> "Tom'S Watch Is On The Table.".istitle() True # アポストロフィー（'）の後の文字が小文字だとFalse >>> "Tom's Watch Is On The Table.".istitle() False # onやtheの先頭文字が小文字だとFalse >>> "Tom'S Watch Is on the Table.".istitle() False 文字コードと文字を相互変換 文字をUnicodeコードポイントに変換: ord() ord()の引数に文字（1文字の文字列）を指定すると、その文字のUnicodeコードポイントが整数intで返される。 i = ord('A') print(i) # 65 print(type(i)) # <class 'int'> Unicodeコードポイントは16進数で表記されることが多い。整数を16進数表記の文字列に変換するには組み込み関数hex()を使う。 例 s = hex(i) print(s) # 0x41 print(type(s)) # <class 'str'> 組み込み関数format()を使うと、ゼロ埋めや0xの有無など、より細かい書式を指定できる。 例 print(format(i, '04x')) # 0041 print(format(i, '#06x')) # 0x0041 特定の文字の16進数表記のUnicodeコードポイントを取得する処理を一気に書くと以下のようになる。 例 print(format(ord('X'), '#08x')) # 0x000058 print(format(ord('?'), '#08x')) # 0x01f4af Unicodeコードポイントを文字に変換: chr() chr()の引数に整数を指定すると、その値がUnicodeコードポイントである文字が文字列strで返される。 print(chr(65)) # A print(type(chr(65))) # <class 'str'> Pythonでは、0xをつけると数値を16進数で記述できるので、16進数表記のUnicodeコードポイントが分かっていればそのままchr()の引数として指定できる。ゼロ埋めされていても問題ない。 print(65 == 0x41) # True print(chr(0x41)) # A print(chr(0x000041)) # A Unicodeコードポイントを表す16進数表記の文字列からそのコードポイントの文字を取得したい場合、16進数表記の文字列を整数intに変換してからchr()に渡す。 16進数表記の文字列を整数intに変換するにはint()を使う。第一引数に文字列、第二引数に基数16を指定する。 s = '0x0041' print(int(s, 16)) # 65 print(chr(int(s, 16))) # A 16進数表記の文字列をint()で整数に変換する場合、文字列に0xがついていれば第二引数は0でもOK。16進数の数値・文字列の扱いについての詳細は以下の記事を参照。 UnicodeコードポイントはU+XXXXの形で記述されることも多い。このような文字列をそのコードポイントの文字に変換するにはスライスで数値部分のみを選択すればよい。 s = 'U+0041' print(s[2:]) # 0041 print(chr(int(s[2:], 16))) # A 文字列をUnicodeコードポイントで記述: \x, \u, \U 文字列リテラルの中では\x, \u, \Uに続けて16進数表記のUnicodeコードポイントを記述でき、その文字として扱われる。 \xXX, \uXXXX, \UXXXXXXXXのように、それぞれ2桁、4桁、8桁の16進数である必要がある。桁数が合っていないとエラー。 print('\x41') # A print('\u0041') # A print('\U00000041') # A print('\U0001f4af') # ? # print('\u041') # SyntaxError: (unicode error) 'unicodeescape' codec can't decode bytes in position 0-4: truncated \uXXXX escape # print('\U0000041') # SyntaxError: (unicode error) 'unicodeescape' codec can't decode bytes in position 0-8: truncated \UXXXXXXXX escape 各コードが1文字として扱われる。文字数を返す組み込み関数len()で確認できる。 print('\u0041\u0042\u0043') # ABC print(len('\u0041\u0042\u0043')) # 3 エスケープシーケンスが無効化されるraw文字列ではそのままの文字の並びとして認識されるので注意。 print(r'\u0041\u0042\u0043') # \u0041\u0042\u0043 print(len(r'\u0041\u0042\u0043')) # 18

df.copy(deep = False) deep = True（デフォルト）：元のdfの行列・値を両方保持するためメモリを消費する。 deep = False：元のdfの行列のみ保持するため、deep = Trueよりもメモリを消費しない。 df.drop inplace = False（デフォルト）：対象のオブジェクトを返すのみ。 inplace = True：削除結果を反映させる。 train_df_shallow = train_df.copy(deep=False) train_df_shallow.drop("Sale",axis = 1,inplace= True) jupyter notebookにおいて、pandasで表示が省略されるのを防ぐ pd.set_option('display.max_columns', 160) # 160番目の列まで全て強制表示 pd.set_option('display.max_columns', None) # 全て強制表示 get_dummiesの使い方 get_dummiesは、文字列のまま使用できるので、いったんlabel列を文字列に戻す。 df['label'] = le.inverse_transform(df['label']) df size price label 1 100 cola 2 150 tea 3 200 coffee get_dummiesでlabel列をダミー変数化して、元のDataFrameに結合する。 df_dummy = pd.get_dummies(df['label']) df2 = pd.concat([df.drop(['label'],axis=1),df_dummy],axis=1) df2 size price coffee cola tea 1 100 0 1 0 2 150 0 0 1 3 200 1 0 0 これでダミー変数化は完成した。 多重共線性 回帰分析などにダミー変数を使用したい場合は、変数同士の相関が高いと多重共線性という問題が発生するため、one-hotエンコーディングの配列から列の1つを削除して使用しなければならない。例えばCoffee列を削除しても、cola列とtea列がともに0であればcoffeeであることがわかるので、情報としての欠落はない。 get_dummieのdrop_firstパラメータにTrueを渡すと最初の列を削除できる。 df_dummy = pd.get_dummies(df['label'],drop_first = True) df3 = pd.concat([df.drop(['label'],axis=1),df_dummy],axis=1) df3 size price cola tea 1 100 1 0 2 150 0 1 3 200 0 0 重回帰分析などでは、カテゴリーの数がN個ある場合、N個のダミー変数を作ると予測できないため、どれか1個の変数を除外します。scikit-learnの重回帰分析 LinearRegression() では、この問題を回避する設計にしているようで、N個のダミー変数を作成しても動作します。 enumerate関数（forループでインデックスを取得できる） Pythonのenumerate()関数を使うと、forループの中でリスト（配列）などのイテラブルオブジェクトの要素と同時にインデックス番号（カウント、順番）を取得できる。 enumerate関数を使ったforループ 引数にリストなどのイテラブルオブジェクトを指定する。 インデックス番号, 要素の順に取得できる。 0以外の数値から開始したい場合は、enumerate()関数の第二引数に任意の開始数値を指定する。 for i, name in enumerate(l, 1): print(i, name) # 1 Alice # 2 Bob # 3 Charlie カレントディレクトリの取得 import os os.getcwd() ファイル選択ダイアログの表示 root = tkinter.Tk() root.withdraw() fTyp = [("", "*.csv")] iDir = os.path.abspath(current_dir) tkinter.messagebox.showinfo('売上可視化プログラム', '処理したいCSVファイルを選択してください！') file = tkinter.filedialog.askopenfilename(filetypes = fTyp,initialdir = iDir) if file == "": print("CSV選択がキャンセルされました。終了します。") sys.exit(1) 文字判定 大文字かどうかを判定（isupper) str.isupper()を使います。 文字列中の英字 全てが大文字、かつ１文字以上の場合はTrue、そうでない場合はFalseを返します。 # 全て英字大文字 >>> 'ABCDE'.isupper() True # 全て英字小文字 >>> 'abcde'.isupper() False # 英字　大文字・小文字混在 >>> 'Abcde'.isupper() False # 英数字混在 >>> 'ABC12'.isupper() True # 英数字以外も混在 >>> 'Aa12\*'.isupper() False 小文字かどうかを判定（islower) str.islower()を使います。 文字列中の英字（上記参照）全てが小文字で、かつそれが1文字以上ある場合はTrue、そうでない場合はFalseを返します。 # 全て英字大文字 >>> 'ABCDE'.islower() False # 全て英字小文字 >>> 'abcde'.islower() True # 英字　大文字・小文字混在 >>> 'Abcde'.islower() False # 英数字混在 >>> 'abc12'.islower() True # 英数字も混在 >>> 'Aa12\*'.islower() False 英数字であることを判定（isalnum） str.isalnum()を使います。 文字列中の全ての文字が英数字で、かつ1文字以上ある場合はTrue、そうでない場合はFalseを返します。 # 全て英字大文字 >>> 'ABCDE'.isalnum() True # 全て英字小文字 >>> 'abcde'.isalnum() True # 英字　大文字・小文字混在 >>> 'Abcde'.isalnum() True # 英数字混在 >>> 'abc12'.isalnum() True # 英数字以外混在 >>> 'Aa12\*'.isalnum() False 英字かどうかを判定（isalpha） str.isalpha()を使います。 文字列中の全ての文字が英字（上記参照）で、かつ1文字以上ある場合はTrue、そうでない場合はFalseを返します。 # 全て英字大文字 >>> 'ABCDE'.isalpha() True # 全て英字小文字 >>> 'abced'.isalpha() True # 英字　大文字・小文字混在 >>> 'Abced'.isalpha() True # 英数字混在 >>> 'abc12'.isalpha() False # 英数字以外混在 >>> 'Aa12\*'.isalpha() False ASCII文字列かどうかを判定（isascii） str.isascii()を使います。 文字列が空、または全ての文字が ASCII（ U+0000 〜 U+007F） である場合はTrue、それ以外の場合はFalseを返します。 Python3.7から追加されたメソッドなのでバージョンに注意。 # 全て英字大文字 >>> 'ABCDE'.isascii() True # 全て英字小文字 >>> 'abced'.isascii() True # 英字　大文字・小文字混在 >>> 'abc12'.isascii() True # 英数字以外混在 >>> 'Aa12\*'.isascii() True # ASCIIコード以外 >>> 'あ'.isascii() False タイトルケース文字列かどうかを判定（istitle） str.istitle()を使います。 タイトルケースとは、文の最初以外の単語についても先頭文字を大文字にすることです。 このメソッドは以下の２つを満たす場合にTrueを返します。 大文字は、英字以外の文字の後のみに置かれる 小文字は、英字の後にだけ続く よって、以下の点に注意が必要です。 アポストロフィー（’）の後の文字が小文字だとタイトルケースと判定されません（例：Tom’s →　Tom’S) 翻訳業界などでは4文字未満の前置詞や接続詞（on, theなど）はこれを適用しないルールがあるようですが、このメソッドはこのルールが適用されません 例 # 文章の最初の文字のみ大文字 >>> "Tom's watch is on the table.".istitle() False # istitle()でTrueとなるタイトルケース >>> "Tom'S Watch Is On The Table.".istitle() True # アポストロフィー（'）の後の文字が小文字だとFalse >>> "Tom's Watch Is On The Table.".istitle() False # onやtheの先頭文字が小文字だとFalse >>> "Tom'S Watch Is on the Table.".istitle() False 文字コードと文字を相互変換 文字をUnicodeコードポイントに変換: ord() ord()の引数に文字（1文字の文字列）を指定すると、その文字のUnicodeコードポイントが整数intで返される。 i = ord('A') print(i) # 65 print(type(i)) # <class 'int'> Unicodeコードポイントは16進数で表記されることが多い。整数を16進数表記の文字列に変換するには組み込み関数hex()を使う。 例 s = hex(i) print(s) # 0x41 print(type(s)) # <class 'str'> 組み込み関数format()を使うと、ゼロ埋めや0xの有無など、より細かい書式を指定できる。 例 print(format(i, '04x')) # 0041 print(format(i, '#06x')) # 0x0041 特定の文字の16進数表記のUnicodeコードポイントを取得する処理を一気に書くと以下のようになる。 例 print(format(ord('X'), '#08x')) # 0x000058 print(format(ord('?'), '#08x')) # 0x01f4af Unicodeコードポイントを文字に変換: chr() chr()の引数に整数を指定すると、その値がUnicodeコードポイントである文字が文字列strで返される。 print(chr(65)) # A print(type(chr(65))) # <class 'str'> Pythonでは、0xをつけると数値を16進数で記述できるので、16進数表記のUnicodeコードポイントが分かっていればそのままchr()の引数として指定できる。ゼロ埋めされていても問題ない。 print(65 == 0x41) # True print(chr(0x41)) # A print(chr(0x000041)) # A Unicodeコードポイントを表す16進数表記の文字列からそのコードポイントの文字を取得したい場合、16進数表記の文字列を整数intに変換してからchr()に渡す。 16進数表記の文字列を整数intに変換するにはint()を使う。第一引数に文字列、第二引数に基数16を指定する。 s = '0x0041' print(int(s, 16)) # 65 print(chr(int(s, 16))) # A 16進数表記の文字列をint()で整数に変換する場合、文字列に0xがついていれば第二引数は0でもOK。16進数の数値・文字列の扱いについての詳細は以下の記事を参照。 UnicodeコードポイントはU+XXXXの形で記述されることも多い。このような文字列をそのコードポイントの文字に変換するにはスライスで数値部分のみを選択すればよい。 s = 'U+0041' print(s[2:]) # 0041 print(chr(int(s[2:], 16))) # A 文字列をUnicodeコードポイントで記述: \x, \u, \U 文字列リテラルの中では\x, \u, \Uに続けて16進数表記のUnicodeコードポイントを記述でき、その文字として扱われる。 \xXX, \uXXXX, \UXXXXXXXXのように、それぞれ2桁、4桁、8桁の16進数である必要がある。桁数が合っていないとエラー。 print('\x41') # A print('\u0041') # A print('\U00000041') # A print('\U0001f4af') # ? # print('\u041') # SyntaxError: (unicode error) 'unicodeescape' codec can't decode bytes in position 0-4: truncated \uXXXX escape # print('\U0000041') # SyntaxError: (unicode error) 'unicodeescape' codec can't decode bytes in position 0-8: truncated \UXXXXXXXX escape 各コードが1文字として扱われる。文字数を返す組み込み関数len()で確認できる。 print('\u0041\u0042\u0043') # ABC print(len('\u0041\u0042\u0043')) # 3 エスケープシーケンスが無効化されるraw文字列ではそのままの文字の並びとして認識されるので注意。 print(r'\u0041\u0042\u0043') # \u0041\u0042\u0043 print(len(r'\u0041\u0042\u0043')) # 18

この記事からは python を使った fits 画像の取り扱いについて書いていきます。今回は、astropy.io.fits, numpy, matplotlib.pyplot, aplpy を用います。例として Spitzer の 8.0 µm データを使用します。 https://irsa.ipac.caltech.edu/data/SPITZER/GLIMPSE/images/I/1.2_mosaics_v3.5/GLON_30-53/ から GLM_03000+0000_mosaic_I4.fits をダウンロードしました。 fits の開き方 まず必要な module を import します。 from astropy.io import fits import numpy as np from matplotlib import pyplot as plt import aplpy もし ModuleNotFoundError: No module named 'aplpy' などと言われてしまった場合は、pip でインストールしましょう。 pip install aplpy astropy.io.fits を使って fits を開きます。path は 絶対path でも 相対path でもどちらでも大丈夫です。 hdu = fits.open("~/your/fits/dir/GLM_03000+0000_mosaic_I4.fits")[0] 後ろの [0] は普通は気にしなくて大丈夫です。hdu には header と data が格納されており、hdu.header や hdu.data と打つとそれぞれ返ってきます。header は辞書、data は numpy array です。 fits の plot まずは最も一般的な matplotlib.pyplot を使って data を plot してみます。 plt.imshow(hdu.data) plt.show() カラーレンジが悪くてほとんど何も見えません。log にしてみましょう。 plt.imshow(np.log10(hdu.data)) plt.show() ちょっと見えました。今は軸が pixel になっているので、これを座標にします。matplotlib でも可能ですが、aplpy という module を使えば簡単にできます。 fig = plt.figure(figsize=(8, 8)) f = aplpy.FITSFigure(hdu, slices=[0], convention='wells', figure=fig) f.show_colorscale(vmin=10, vmax=1000, stretch='log', cmap="Greens", aspect="equal") plt.show() ticks が外向きなのが気になる人は以下を実行しましょう。 plt.rcParams['xtick.direction'] = 'in' plt.rcParams['ytick.direction'] = 'in' 拡大図は recenter を使います。横軸の中心 (degree), 縦軸の中心 (degree), width (degree), height (degree) の順番で渡します。 fig = plt.figure(figsize=(8, 8)) f = aplpy.FITSFigure(hdu, slices=[0], convention='wells', figure=fig) f.show_colorscale(vmin=10, vmax=1000, stretch='log', cmap="Greens", aspect="equal") f.recenter(30.5, 0.0, width=1.0, height=1.0) plt.show() その他 カラーバー, スケールバー, マーカーなど一通りの機能ももちろん備わっています。ドキュメントを参照して描きたい図を作りましょう。 https://aplpy.github.io/ 3D fits (cube) の場合も基本的には同様で、slices=[0] の 0 部分を plot したいチャンネルに変更します。 よくあるエラー ALMA の cube データがうまく開けない！ ALMA データはデフォルトでは 4軸目に Stokes parameter が入っています。しかも3軸目がfrequencyです。 casaで一回読み込んで出力すると大抵うまくいきます。casa のターミナルを起動して、 importfits(fitsimage="yourfits.fits", imagename="yourfits.im") exportfits(imagename="yourfits.im/", fitsimage="yourfits_new.fits", velocity=True, dropstokes=True, overwrite=True) で出てきた fits を使いましょう。 リンク 目次

はじめに 今回はDjangoを使って友達リスト的なアプリを作成しようと思います。 プロジェクト名はTomodachi 環境構築(プロジェクト作成) # 仮想環境 $ cd 作業スペース $ mkdir プロジェクト名 $ cd プロジェクト名 $ python -m venv venv $ source venv/bin/activate # Djangoプロジェクト作成 $ pip install django $ pip install --upgrade pip $ django-admin startproject config . なぜか何回かdjangoをインストールすることができなかったけど $ python -m pip freezeとか $ python -m pip listとかしていたら完了した。 アプリ作成とビュー設定 # アプリケーション作成 $ python manage.py startapp friendslist config/urls.py # ルーティングを実装 from friendslist import views urlpatterns = [ path('admin/', admin.site.urls), path('', views.index) ] config/settings.py # 開発環境の設定 import os　　←☆追記 ALLOWED_HOSTS = ['*']　　←☆追記 INSTALLED_APPS = [ 'django.contrib.admin', 'django.contrib.auth', 'django.contrib.contenttypes', 'django.contrib.sessions', 'django.contrib.messages', 'django.contrib.staticfiles', 'friendslist' ←☆追記 ] friendslist/views.py # ルーティング設定が出来ているかを確認する from django.http import HttpResponse def index(request): return HttpResponse('Hello Friend') # サーバー起動 $ python manage.py runserver テンプレート作成とビュー設定 # templateディレクトリとfriendslistディレクトリを作成 $ mkdir templates $ mkdir templates/friendslist config/settings.py # TEMPLATESのDIRを設定する TEMPLATES = [ { 'BACKEND': 'django.template.backends.django.DjangoTemplates', 'DIRS': [ os.path.join(BASE_DIR, 'templates'),　　←☆追記する ], 'APP_DIRS': True, 'OPTIONS': { 'context_processors': [ 'django.template.context_processors.debug', 'django.template.context_processors.request', 'django.contrib.auth.context_processors.auth', 'django.contrib.messages.context_processors.messages', ], }, }, ] friendslist/views.py # renderの第三引数に変数を追加してテンプレートで使う def index(request): context = { 'title': 'Tomodachi' } return render(request, 'friendslist/index.html', context) templates/friendslist/base.html # 共通テンプレートを作成する <!DOCTYPE html> <html lang="ja"> <head> <meta charset="UTF-8"> <meta http-equiv="X-UA-Compatible" content="IE=edge"> <meta name="viewport" content="width=device-width, initial-scale=1.0"> <title>{{ title }}</title> </head> <body> {% block content %} {% endblock %} </body> </html> templates/friendslist/index.html # base.htmlを読み込んで表示する {% extends 'friendslist/base.html' %} {% block content %} <h1>こんにちは、{{ title }}</h1> {% endblock %} テンプレートにBootstrapを導入する friendslist/base.html # Bootstrapを導入する <!DOCTYPE html> <html lang="ja"> <head> <meta charset="UTF-8"> <meta http-equiv="X-UA-Compatible" content="IE=edge"> <meta name="viewport" content="width=device-width, initial-scale=1.0"> <link href="https://cdn.jsdelivr.net/npm/bootstrap@5.0.0-beta1/dist/css/bootstrap.min.css" rel="stylesheet" integrity="sha384-giJF6kkoqNQ00vy+HMDP7azOuL0xtbfIcaT9wjKHr8RbDVddVHyTfAAsrekwKmP1" crossorigin="anonymous">　　←☆追記する <title>{{ title }}</title> </head> <body> {% block content %} {% endblock %} <script src="https://cdn.jsdelivr.net/npm/bootstrap@5.0.0-beta1/dist/js/bootstrap.bundle.min.js" integrity="sha384-ygbV9kiqUc6oa4msXn9868pTtWMgiQaeYH7/t7LECLbyPA2x65Kgf80OJFdroafW" crossorigin="anonymous"></script>　　←☆追記する <script src="https://cdn.jsdelivr.net/npm/@popperjs/core@2.5.4/dist/umd/popper.min.js" integrity="sha384-q2kxQ16AaE6UbzuKqyBE9/u/KzioAlnx2maXQHiDX9d4/zp8Ok3f+M7DPm+Ib6IU" crossorigin="anonymous"></script>　　←☆追記する <script src="https://cdn.jsdelivr.net/npm/bootstrap@5.0.0-beta1/dist/js/bootstrap.min.js" integrity="sha384-pQQkAEnwaBkjpqZ8RU1fF1AKtTcHJwFl3pblpTlHXybJjHpMYo79HY3hIi4NKxyj" crossorigin="anonymous"></script>　　←☆追記する </body> </html> friendslist/index.html # 見た目作成 {% extends 'friendslist/base.html' %} {% block content %} <nav class="navbar navbar-expand-lg navbar-light bg-warning mb-5"> <div class="container-fluid"> <a class="navbar-brand" href="#">{{ title }}</a> <button class="navbar-toggler" type="button" data-bs-toggle="collapse" data-bs-target="#navbarNav" aria-controls="navbarNav" aria-expanded="false" aria-label="Toggle navigation"> <span class="navbar-toggler-icon"></span> </button> <div class="collapse navbar-collapse" id="navbarNav"> <ul class="navbar-nav"> <li class="nav-item"> <a class="nav-link active" aria-current="page" href="#">友達一覧</a> </li> <li class="nav-item"> <a class="nav-link" href="#">友達登録</a> </li> <li class="nav-item"> <a class="nav-link" href="#">カテゴリ一覧</a> </li> <li class="nav-item"> <a class="nav-link disabled" href="#" tabindex="-1" aria-disabled="true">カテゴリ登録</a> </li> </ul> </div> </div> </nav> <main class="container mb-5"> <div class="row"> <div class="col col-md-4"> <h2>カテゴリ一覧</h2> <ul class="list-group"> <li class="list-group-item active" aria-current="true">仕事</li> <li class="list-group-item">小学校</li> <li class="list-group-item">中学校</li> <li class="list-group-item">高校</li> <li class="list-group-item">大学</li> </ul> </div> <div class="column col-md-8"> <h2>友達一覧</h2> <ul class="list-group"> <li class="list-group-item active" aria-current="true">山田太郎</li> <li class="list-group-item">山田花子</li> <li class="list-group-item">山田一郎</li> <li class="list-group-item">山口二郎</li> <li class="list-group-item">山口三郎</li> </ul> </div> </div> </main> <div class="container mb-5"> <div class="row"> <div class="col"> <a class="btn btn-danger" href="#" role="button">友達一覧</a> </div> <div class="col"> <a class="btn btn-success" href="#" role="button">友達登録</a> </div> <div class="col"> <a class="btn btn-warning" href="#" role="button">カテゴリ一覧</a> </div> <div class="col"> <a class="btn btn-info" href="#" role="button">カテゴリ登録</a> </div> </div> </div> <footer class="border-top py-4 text-center bg-light fixed-bottom"> <p> <a href="/">トップに戻る</a> </p> </footer> {% endblock %} 共通部品でリファクタリング friendslist/snippets/header.html <nav class="navbar navbar-expand-lg navbar-light bg-warning mb-5"> <div class="container-fluid"> <a class="navbar-brand" href="#">{{ title }}</a> <button class="navbar-toggler" type="button" data-bs-toggle="collapse" data-bs-target="#navbarNav" aria-controls="navbarNav" aria-expanded="false" aria-label="Toggle navigation"> <span class="navbar-toggler-icon"></span> </button> <div class="collapse navbar-collapse" id="navbarNav"> <ul class="navbar-nav"> <li class="nav-item"> <a class="nav-link active" aria-current="page" href="#">友達一覧</a> </li> <li class="nav-item"> <a class="nav-link" href="#">友達登録</a> </li> <li class="nav-item"> <a class="nav-link" href="#">カテゴリ一覧</a> </li> <li class="nav-item"> <a class="nav-link disabled" href="#" tabindex="-1" aria-disabled="true">カテゴリ登録</a> </li> </ul> </div> </div> </nav> friendslist/snippets/footer.html <footer class="border-top py-4 text-center bg-light fixed-bottom"> <p> <a href="/">トップに戻る</a> </p> </footer> friendslist/base.html {% include 'friendslist/snippets/header.html' %} {% include 'friendslist/snippets/footer.html' %} 管理機能実装 friendslist/models.py # カスタムユーザーモデル作成(定型文) from django.db import models from django.contrib.auth.models import BaseUserManager, AbstractBaseUser class UserManager(BaseUserManager): def create_user(self, email, password=None): if not email: raise ValueError('Users must have an email address') user = self.model( email=self.normalize_email(email), ) user.set_password(password) user.save(using=self._db) return user def create_superuser(self, email, password=None): user = self.create_user( email, password=password, ) user.is_admin = True user.save(using=self._db) return user class User(AbstractBaseUser): email = models.EmailField( max_length=255, unique=True, ) is_active = models.BooleanField(default=True) is_admin = models.BooleanField(default=False) objects = UserManager() USERNAME_FIELD = 'email' REQUIRED_FIELDS = [] def __str__(self): return self.email def has_perm(self, perm, obj=None): "Does the user have a specific permission?" # Simplest possible answer: Yes, always return True def has_module_perms(self, app_label): "Does the user have permissions to view the app `app_label`?" # Simplest possible answer: Yes, always return True @property def is_staff(self): "Is the user a member of staff?" # Simplest possible answer: All admins are staff return self.is_admin config/settings.py # 管理者ユーザー設定を追記する AUTH_USER_MODEL = 'friendslist.User' friendslist/admin.py # 管理画面の表示を編集する(定型文) from django.contrib import admin from django.contrib.auth.admin import UserAdmin from django.contrib.auth.models import Group from friendslist.models import User class CustomUserAdmin(UserAdmin): fieldsets = ( (None, { 'fields': ( 'email', 'password', ) }), (None, { 'fields': ( 'is_active', 'is_admin', ) }) ) list_display = ('email', 'is_active') list_filter = () ordering = () filter_horizontal = () add_fieldsets = ( (None, { 'fields': ('email', 'password',), }), ) admin.site.unregister(Group) admin.site.register(User, CustomUserAdmin) # マイグレーションをする $ python manage.py makemigrations $ python manage.py migrate config/urls.py # 管理画面にパスを通す(元から記述されている) path('admin/', admin.site.urls), # superuserを作る $ python manage.py createsuperuser さいごに 今回は プロジェクト作成 アプリ作成 テンプレート作成 ビュー設定 管理画面設定

はじめに 昔作成したチートシートを転記したもので、参照元は不明です。 また、私の資料整理のために記載しています。ご容赦ください。 決定係数 決定係数とは、回帰分析によって求められた目的変数の予測値が、実際の目的変数の値とどのくらい一致しているかを表している指標。１に近いほど分析の精度が高いことを表す。 「（自由度調整済み）決定係数さえ高ければなんでもよい」という考え方はもちろん不適切。 分析を行うからにはその目的があり、検証したい仮説があるはずです。こうした目的や仮説とは一切関係ない変数を、決定係数が上がるからという理由だけで無秩序に分析に加えてしまうと、そのモデルは解釈が難しくなり当初の目的を果たせないでしょう。 また決定係数はあくまで「予測の当てはまりの良さ」を表す指標です。分析の目的が「ある変数の値を予測したい」の場合には適切な指標ですが、「雨が降った日数がコンビニの月間の売上に影響があるかどうか知りたい」のように、ある変数の影響の有無が主眼であり予測は重視しない場合には、決定係数に注目することはあまり意味がないと言えます。この場合には「雨が降った日数」という説明変数の回帰係数や、その係数に対する検定の結果（有意かどうか）にまず注目すべきでしょう。 決定係数は分かりやすく便利な指標であるため、分析の際にはどうしてもここだけに目が行きがちです。しかしながら、本当に価値のある分析をするためには、数字の意味をしっかり理解した上で、分析の目的と照らし合わせて正しく使うことが求められます。 線形回帰とロジスティック回帰 異なる点 線形回帰は実際の値を予測するのに対して、ロジスティック回帰は発生確率（0~1の値）を予測する つまり握力の例で言えば、線形回帰は「年齢から平均握力そのものを予測する」のに対して、ロジスティック回帰は「（例えば）30歳での平均握力が30kgとなる確率」を予測します。 実際のデータ分析での流れ では、実際のデータ分析ではどのような流れで処理されることが多いのでしょうか。 機械学習をする場合、元データが存在した場合には、それをまずはPandasのDataFrameに落とし込んで加工します。そして、DataFrame内の必要なデータをNumPyのndarrayに変換し、機械学習の演算や高度なアルゴリズムによる処理を行う、というのがよくある方法だと思います。 機械学習や画像処理を含む高度なアルゴリズムはDataFrameではなく、ndarrayを引数にとることが多いです。NumPyに落とし込むのは最後の最後の段階で、NumPyが処理するのが得意な形に持っていくようPandasで色々前処理をする役割になることが多いです。 データ分析では作業時間の5割以上は前処理に割くといっても過言ではありません。そう考えると、Pandasを使った処理を覚えておくと役立つ場面が多いでしょう。 Feature Scalingとは 特徴量の取りうる値の範囲（スケール）を変えることです。 データセットの特徴量間でスケールが異なることは多々あります。例えば、体重と身長、家の価格と部屋数では、その単位と値の範囲が異なります。 そのような特徴量間で異なるスケールのデータセットをモデルで学習させた場合、うまく学習できないということがおこるのです。 そのため、学習前の前処理で、特徴量間のスケールを揃える必要があります。 Feature Scalingの種類 主なFeature Scalingの種類として、以下の2つがあります。 正規化(normalization) 標準化(standardization) 正規化(normalization)とは 正規化とは、特徴量の値の範囲を一定の範囲におさめる変換になります。主に[0, 1]か、[-1, 1]の範囲内におさめることが多いです。 標準化(standardization)とは 標準化とは、特徴量の平均を0、分散を1にする変換になります。 正規化と標準化の使い分け 基本は標準化を用います。理由は、正規化の場合、外れ値が大きく影響してしまうためです。 ただし、画像データの場合は学習コストを下げるため、[0,1]の範囲に収まるよう255.0で割ることで正規化するのが一般的です。 正規化 使うとき: - 画像処理におけるRGBの強さ[0,255] - sigmoid, tanhなどの活性化関数を用いる、NNのいくつかのモデル 標準化 使うとき: - ロジスティック回帰、SVM、NNなど勾配法を用いたモデル - kNN, k-meansなどの距離を用いるモデル - PCA, LDA(潜在的ディリクレ配分法), kernel PCA などのfeature extractionの手法 使わないとき: 決定木、ランダムフォレスト 機械学習コード Scikit-Learn # 基本のライブラリインポート手順 # アルゴリズム、データセット、などのライブラリをインポート from sklearn import neighbors, datasets, preprocessing # トレーニング用とテスト用に分けてくれる便利な関数 from sklearn.model_selection import train_test_split # 予測に対する正解率（accuracy）を出すために必要 from sklearn.metrics import accuracy_score # 基本の学習方法手順 # データセットからアヤメの花のデータを取り出す。 iris = datasets.load_iris() # データの列0番目から1番目までの全行をXに入れて、 # データの目的変数をyに入れる。 X, y = iris.data[:, :2], iris.target # トレーニング用とテスト用に分ける。 X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=33) # StandardScalerを使ってデータセットを標準化する。 scaler = preprocessing.StandardScaler().fit(X_train) # 実際の変換作業。 X_train = scaler.transform(X_train) X_test = scaler.transform(X_test) # K近傍法を使ってモデルを作成。 knn = neighbors.KNeighborsClassifier(n_neighbors=5) # トレーニング用データで学習をさせる。 knn.fit(X_train, y_train) # テスト用データで予測値を出す。 y_pred = knn.predict(X_test) # 正解率をテスト用目的変数と予測値を使って計算する。 accuracy_score(y_test, y_pred) Linear Regression（線形回帰） from sklearn.linear_model import LinearRegression linear_model = linear_model.LinearRegression() linear_model.fit(X_train,y_train) # R2(決定係数) linear_model.score(X_test,y_test) # coefficient(偏回帰係数) print('偏回帰係数: ', linear_model.coef_) # indercept(切片) print('切片: ', linear_model.intercept_) #予測値を出す。 predicted = linear_model.predict(X_test) LogisticRegression（ロジスティック回帰） 線形回帰と似てますが、目的変数が０か１、YES・NOなどの時に使います。 from sklearn.linear_model import LogisticRegression model = model.LogisticRegression() model.fit(X_train,y_train) # R2(決定係数) model.score(X_test,y_test) # coefficient(偏回帰係数) print('偏回帰係数: ', model.coef_) # indercept(切片) print('切片: ', model.intercept_) #予測値を出す。 predicted = model.predict(X_test) Decision Tree（決定木） from sklearn import tree # 不純度の指標をジニ係数にする。エンロトピーと使い分ける必要がある。 model = tree.DecisionTreeClassifier(criterion='gini') model.fit(X_train, y_train) model.score(X_test, y_test) predicted = model.predict(X_test) SVM (サポートベクターマシン) from sklearn import svm # クラス分類問題において、データの数がそこまで大きくない場合は # SVC(Linear SVM)を使います。 model = svm.svc() model.fit(X_train, y_train) model.score(X_test, y_test) predicted= model.predict(X_test) Naive Bayes (ナイーブベイズ) # 数種類あるうちのガウシアン分布（正規分布） # を使ったベイズ分類器をインポート from sklearn.naive_bayes import GaussianNB model.fit(X_train, y_train) predicted= model.predict(X_test) kNN (K近傍法) from sklearn.neighbors import KNeighborsClassifier # n_neighbors:Kの数はデフォルトでは5に設定されている。 KNeighborsClassifier(n_neighbors=6) model.fit(X_train, y_train) predicted= model.predict(X_test) K-Means（K平均法） from sklearn.cluster import KMeans # n_clusters:何個にクラス分けするか、random_stateは0に設定。 k_means = KMeans(n_clusters=3, random_state=0) model.fit(X) predicted= model.predict(x_test) Random Forest（ランダムフォレスト） from sklearn.ensemble import RandomForestClassifier # 決定木の深さを100に設定 model= RandomForestClassifier(max_depth=100) model.fit(X, y) predicted= model.predict(x_test) # 深さ100の時の各特徴量の重要性。 print(model.feature_importances_) Dimensionality Reduction Algorithms（次元削減） # PCA (線形アルゴリズム)を使用 from sklearn.decomposition import PCA # n_components:2次元に次元を削減する pca = PCA(n_components=2) # トレーニング用のデータセットの次元をPCAを用いて削減する。 train_reduced = pca.fit_transform(train) # テスト用のデータセットの次元をPCAを用いて削減する。 test_reduced = pca.transform(test) Gradient Boosting Algorithms（勾配ブースティング） from sklearn.ensemble import GradientBoostingClassifier # n_estimators: 決定木の数 # learning_rate: 結果に対する各決定木の影響度合。小さい値が良いとされているが、 # 大きい数値で試した後小さい数値にして試すのがセオリー。 # max_depth: 決定木の深さ。 model= GradientBoostingClassifier(n_estimators=100, learning_rate=1.0, max_depth=1, random_state=0) model.fit(X, y) predicted= model.predict(x_test) XGBoost XGBoostとは、先程のGradient BoostingとRandom Forestsを組み合わせたアンサンブル学習である。 アンサンブル学習とは、複数のモデル（学習器）を使い、1つの学習モデルを作成する手法。 pythonパッケージのインストールは以下のURLを参考に。 from xgboost import XGBClassifier from sklearn.model_selection import train_test_split from sklearn.metrics import accuracy_score X = dataset[:,0:10] Y = dataset[:,10:] seed = 1 X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size=0.33, random_state=seed) model = XGBClassifier() model.fit(X_train, y_train) y_pred = model.predict(X_test) Scikit-Learn目的別 データ前処理 # Standardization(標準化) from sklearn.preprocessing import StandardScaler scaler = StandardScaler().fit(X_train) standardized_X = scaler.transform(X_train) standardized_X_test = scaler.transform(X_test) # Normalization(正規化) from sklearn.preprocessing import Normalizer scaler = Normalizer().fit(X_train) normalized_X = scaler.transform(X_train) normalized_X_test = scaler.transform(X_test) # Binarization（二値化） from sklearn.preprocessing import Binarizer binarizer = Binarizer(threshold=0.0).fit(X) binary_X = binarizer.transform(X) # Encoding Categorical Features（文字列要素を数値コードに変換） from sklearn.preprocessing import LabelEncoder enc = LabelEncoder() y = enc.fit_transform(y) # Imputing Missing Values（欠損値の穴埋め。） from sklearn.preprocessing import Imputer # strategyで平均か中央値か最頻値を選択。 imp = Imputer(missing_values=0, strategy='mean', axis=0) imp.fit_transform(X_train) # Generating Polynomial Features (多項式の特徴量を作る) from sklearn.preprocessing import PolynomialFeatures poly = PolynomialFeatures(5) poly.fit_transform(X) モデルの成果を評価する 分類結果を評価 # Accuracy Score (正解率) knn.score(X_test, y_test) from sklearn.metrics import accuracy_score accuracy_score(y_test, y_pred) # Classification Report (結果の詳細なレポートを表示) from sklearn.metrics import classification_report print(classification_report(y_test, y_pred)) # Confusion Matrix (混同行列) # 正しく識別できた件数、誤って識別した件数を比較 from sklearn.metrics import confusion_matrix print(confusion_matrix(y_test, y_pred)) 回帰結果を評価 # Mean Absolute Error (平均絶対誤差) from sklearn.metrics import mean_absolute_error y_true = [5, -0.25, 3] mean_absolute_error(y_true, y_pred) # Mean Squared Error (平均二乗誤差) from sklearn.metrics import mean_squared_error mean_squared_error(y_test, y_pred) # R² Score (決定係数) from sklearn.metrics import r2_score r2_score(y_true, y_pred) クラスタリング結果を評価 # Adjusted Rand Index (ランド指数) from sklearn.metrics import adjusted_rand_score adjusted_rand_score(y_true, y_pred) # Homogeneity (均一性) from sklearn.metrics import homogeneity_score homogeneity_score(y_true, y_pred) # V-measure from sklearn.metrics import v_measure_score metrics.v_measure_score(y_true, y_pred) 交差検証 from sklearn.cross_validation import cross_val_score print(cross_val_score(knn, X_train, y_train, cv=4)) print(cross_val_score(lr, X, y, cv=2)) 学習モデルの調整 グリッドサーチ from sklearn.grid_search import GridSearchCV # 最適化したいパラメータを設定 params = {"n_neighbors": np.arange(1,3), "metric": ["euclidean", "cityblock"]} # パラメータを最適化 grid = GridSearchCV(estimator=knn, param_grid=params) grid.fit(X_train, y_train) print(grid.best_score_) print(grid.best_estimator_.n_neighbors ランダムサーチ from sklearn.grid_search import RandomizedSearchCV # 最適化したいパラメータを設定 params = {"n_neighbors": range(1,5), "weights": ["uniform", "distance"]} rsearch = RandomizedSearchCV(estimator=knn, param_distributions=params, cv=4, n_iter=8, random_state=5) rsearch.fit(X_train, y_train) print(rsearch.best_score_) Numpy # Numpy配列を作り、同時に2行4列にする。 np.array([1,2,3,5,6,7,8,9]).reshape(2,4) # ゼロ埋め np.zeros((3,4)) # 1で埋める np.ones((2,3,4),dtype=np.int16) # ランダムな要素をもつ4行5列の配列を作る np.empty((4, 5)) # 10から25までの整数を5間隔の配列にする np.arange(10,25,5) # 0から2までを等間隔で9に分割 np.linspace(0,2,9) # 正規分布の乱数を100個生成 np.random.randn(100) # 20から100までの整数の中から一つだけランダムで選ぶ np.random.randint(20,100) # 転置 # bは3行5列の配列 t = np.T(b) # tは5行3列の配列に入れ替わる # 行列の掛け算 np.dot(t,b) # 配列の連結 # axis=1は列方向への結合 # axis=0は行方向への結合 a=np.array([[1,2,3,4,5,6],[1,2,3,4,5,6]]) d=np.array([[1,2,3,4,5,6],[1,2,3,4,5,6]]) np.concatenate((a,d),axis=1) np.concatenate((a,d),axis=0) # 行方向への結合 np.vstack((a,b)) # 列方向への結合 np.hstack((e,f)) # 列方向に指定した数で等分する。 np.hsplit(matrix_arr,3) # 行方向に指定した数で等分する。 np.vsplit(matrix_arr,2) pandas # データフレームを作る columns = ['A','B','C'] rows = [[1,2,3],[4,5,6],[7,8,9]] df = pd.DataFrame(rows,columns=columns) # ファイルを読み込む。 df = pd.read_csv('file.csv') df = pd.read_json('file.json') df = pd.read_excel('file.xlsx') # ファイルを書き込む。 df.to_csv('file.csv') df.to_json('file.json') df.to_excel('file.xlsx') # URLからでも読み込める df = pd.read_csv('http://example/bacon/spam/ham/egg/file.csv') # エクセルシートで指定してデータフレームを作る。 xlsx = pd.ExcelFile('file.xlsx') df = pd.read_excel(xlsx, 'Sheet1') # 行と列を整数で指定 df.iloc[[0],[0]] # 列を文字列で指定 df.loc[[0], ['Country']] # 条件を指定して行や列の値を取る df[df['A']>2] df[df['B']==1] # 行や列に関数を働かせる df['A'].apply(str) df['A'].apply(lambda x: str(x) if x>2 else x) # 欠損値を見つけ、欠損値を削除 df.isnull() # すべてが欠損値の行であれば、その行は削除される。 df.dropna(how='all') # axis=1なら行方向が削除される。 # つまり、すべてが欠損値の列であれば、その列は削除される。 df.dropna(how='all', axis=1) # howにanyを指定すると、 # 欠損値が一つでも含まれる行が削除される。 df2.dropna(how='any') # axis=1にすると列に適用。つまり行方向に適用。 # 欠損値が一つでも含まれる列が削除される。 df2.dropna(how='any', axis=1) # 新しい列を加える。 # 目的変数や分析結果を加えるさいに便利。 df['new_columns'] = new_row # 基本的だが本当によく使う # 行列数を表示 df.shape # インデックスの情報 df.index # 列の情報を表示 df.columns # 欠損値やオブジェクトタイプを調べるさいに便利 df.info() # ユニークなオブジェクトの個数をカウント df.count() # axisなどで行方向や列方向も指定できるものがある。 # すべての列を足す df.sum() # 一つ上の行から順番に足して表示 df.cumsum() # 最小値、最大値を表示 df.min()/df.max() # 最大値・最小値の要素のインデックス（ラベル）が取得できる。 df.idxmin()/df.idxmax() # 平均値や偏差や最大値などの数値情報がまとめて表示される。 df.describe() # 平均値を求める。 df.mean() # 中央値を取得する df.median()

この記事では Montage の使い方について紹介します。 Montage のインストール http://montage.ipac.caltech.edu/docs/download.html からドキュメントに従って自身の環境に合わせてインストールします。 そしてパスを通し、ターミナルでどこからでも「mImgtbl」などのコマンドが使える状態であることを確認します。 データの準備 今回は例として、Spitzer 24 µm のアーカイブデータを使用します。 https://irsa.ipac.caltech.edu/data/SPITZER/MIPSGAL/ から、W43 付近のデータをダウンロードします。 または、 https://irsa.ipac.caltech.edu/data/SPITZER/MIPSGAL/images/mosaics24/ からダウンロードすることも可能です。 MG0290n005_024.fits MG0290p005_024.fits MG0300n005_024.fits MG0300p005_024.fits MG0310n005_024.fits MG0310p005_024.fits の計６つの fits をつなぎ合わせたいと思います。 Montage 使用方法 (2D) まず、「imagedir」というディレクトリを作成し、その中に fits を全て入れます。 ターミナルで fits のあるディレクトリに移動し、 mkdir imagedir mkdir resultdir mv MG*_024.fits imagedir/ と実行します。 そして mImgtbl imagedir images.tbl を実行します。images.tbl という名前のファイルが生成されます。ファイルのリストです。 次に、 mMakeHdr images.tbl template.hdr GAL を実行します。template.hdr という名前のファイルが生成されます。新しいヘッダーです。「GAL」を後ろにつけることで座標系を銀経銀緯にできます。 次に、 mProjExec -p imagedir images.tbl template.hdr resultdir stats.tbl を実行します。resultdir の中にreprojectされた fits が生成されます。少し計算に時間がかかります。 次に、 mImgtbl resultdir/ resultimages.tbl を実行します。resultimages.tbl という名前のファイルが生成されます。 次に、 mAdd -p resultdir/ resultimages.tbl template.hdr W43_24um.fits を実行します。resultdir の中の fits が足し合わされます。W43_24um.fits は保存する名前です。少し計算に時間がかかります。W43_24um_area.fits というものも生成されますが、これは中間ファイルのようなものなので消して大丈夫です。 生成された fits の値はfloat64で書かれています。python の astropy.io.fits などで開いてfloat32にすることで約半分まで軽くすることができます。 最後に、生成された fits を ds9 で確認します。 Montage 使用方法 (3D) 3D の fits (cube) 同士を足し合わせたいこともあるかと思います。 一応、Montage は cube にも対応していますが、速度軸がちょっとおかしな挙動をすることがあるので、できれば casa の「imregrid」などで速度軸を揃えてから足し合わせることをお勧めします。 casa の使い方については別の記事を参照してください。 mImgtbl imagedir/ images.tbl mMakeHdr images.tbl template.hdr for infile in raw/*.fits; do mProjectCube $infile resultdir/${infile##*/}_proj.fits template.hdr ; done mImgtbl resultdir/ resultimages.tbl mAddCube -p resultdir resultimages.tbl template.hdr result.fits 以上です。 リンク 目次

まずはデータを目で見るところからスタートしましょう。 今回は例として Spitzer のアーカイブデータを使用します。 https://sha.ipac.caltech.edu/applications/Spitzer/SHA/ から、自分の好きな領域のデータを検索してダウンロードします。 銀河面であれば、 https://irsa.ipac.caltech.edu/data/SPITZER/GLIMPSE/images/I/1.2_mosaics_v3.5/ からダウンロードすることも可能です。 例えば、 GLON_30-53/ → GLM_03000+0000_mosaic_I1.fits GLM_03000+0000_mosaic_I4.fits をクリックして、W43 付近の 3.6 μm (IRAC1) と 8.0 μm (IRAC4) のデータをダウンロードします。 ついでに 24 µm のデータもこちらからダウンロードします。 https://irsa.ipac.caltech.edu/data/SPITZER/MIPSGAL/ ここでは、「W43」で検索して出てきた MG0310p005_024.fits をダウンロードします (左にダウンロードボタンがあります)。 https://irsa.ipac.caltech.edu/data/SPITZER/MIPSGAL/images/mosaics24/ からもダウンロードすること可能です。 fits 画像を開くには、専用のソフトが必要です。色々ありますが、一番おすすめなのが、ds9 です。 ds9 をダウンロード&インストールして開きます。 ds9の基本的な使い方 fits の読み込み 「ファイル」→「開く」からfitsを選択します。 ds9の場所にパスが通っている場合は、ターミナルから ds9 GLM_03000+0000_mosaic_I1.fits とするとds9の起動とファイルのopenを同時に行えます。 ds9 GLM_03000+0000_mosaic_I1.fits GLM_03000+0000_mosaic_I4.fits で複数のファイルを開くこともできます。 GLM_03000+0000_mosaic_I1.fitsを開くとこのような画像が出てきます。 ヘッダーの確認 真ん中のパネルの「ファイル」から「ヘッダ」を押すと見ることができます。 拡大縮小・移動 マウスのホイールもしくは Mac であれば二本指の上下で拡大縮小ができます。 「ズーム」からでも変更可能です。 右上の画面の水色の四角を移動させると中心の位置を移動できます。 または、真ん中のパネルの「編集」から「パン」を選ぶとクリックしたところが真ん中にくるモードになります。 カラースケールの変更 場合によっては、logスケールの方が良いこともあります。 上の「スケール」から「対数」を選択しましょう。 また、カラーの最小値最大値は「スケール」から一番下の「スケールの設定」から行えます。 カラーマップの変更 上の「色」から変更可能です。個人的なおすすめは sls です。 座標系の変更 デフォルトでは fk5 (赤経赤緯) です。今回は銀河面のデータなので、上の「WCS」から「銀河座標」を選びましょう。 フレームの追加 フレームを追加することで複数の fits を同時に見ることができます。 真ん中のパネルの「フレーム」から「新規」をクリックします。 同じく「フレーム」から「並べて表示 (英語版だと tile)」を選択します。 フレーム間での座標等の固定 上の「フレーム」から「固定」→「フレーム」→「WCS」で全フレームの表示している座標が fix されます。 このような画面になったと思います。 他にも ds9 にはコントアを引いたり、マーカー (領域) を打ったり様々なことができます。 RGB 合成図の表示 真ん中のパネルの「フレーム」から「rgb色」をクリックします。そうすると小さなウィンドウが出てくると思います。 「選択中」を赤色にした状態で、「ファイル」→「開く」で MG0310p005_024.fits を選択します。 次に、「選択中」を緑色にした状態で、「ファイル」→「開く」で GLM_03000+0000_mosaic_I4.fits を選択します。 次に、「選択中」を青色にした状態で、「ファイル」→「開く」で GLM_03000+0000_mosaic_I1.fits を選択します。 各色のスケールを、「スケール」から一番下の「スケールの設定」で好みに合わせて調節すると、以下のような map が描けます。 基本的な使い方は以上となります。 今回の例だと、赤色 (24 µm) の fits だけ領域が狭いです。 隣り合う領域の fits もダウンロードし結合して1つの fits にしたい場合が多いと思います。そのような際には Montage が非常に便利で信頼性も高いです。 リンク 目次

はじめに この天文データ解析入門シリーズでは、天文に関する一般的な可視化ソフトの使い方から機械学習を用いた解析まで幅広く取り扱っていきます。内容は電波天文学が中心となりますが、可視や赤外線のデータなども取り扱います。 準備するもの 本シリーズを進めていく上で必要なツール・あると良いツールを紹介します。 必須 ds9 軽くて圧倒的に便利な無料のfitsの可視化ソフトです。公式のページから自分の環境にあったものをダウンロードしてインストールしてください。 https://sites.google.com/cfa.harvard.edu/saoimageds9 Python3系 本シリーズで基本とする言語です。整っていない方は方法をgoogle検索して是非整えましょう。 pip Python の module (パッケージ) を簡単にインストールするものです。整っていない方は方法をgoogle検索して是非整えましょう。 Jupyter notebook 簡単にいうとPythonのコードを書いて実行する非常に便利なノートです。整っていない方は方法をgoogle検索して是非整えましょう。 Python modules numpy (標準でインストール済み) scipy (標準でインストール済み) matplotlib (標準でインストール済み) astropy aplpy pandas Python module は、pipが正常にインストールされていれば、ターミナルから pip install astropy pip install aplpy で簡単にインストールできます。 何らかの解析サーバーを使用していて権限がない場合は、 pip install --user <パッケージ名> で大抵はインストールできます。 その他あると良いもの casa 生データからのreductionや、regridなど様々な場面で有効です。 https://casa.nrao.edu/ Montage fits の合成に使用します。 http://montage.ipac.caltech.edu/ miriad fits の回転などに便利です。 https://www.atnf.csiro.au/computing/software/miriad/ Python modules astrodendro (天体の構造同定) pycupid (天体の構造同定) 目次 天文データ解析入門 その1 (ds9の使い方) 天文データ解析入門 その2 (Montageの使い方) 天文データ解析入門 その3 (astropyでのfitsの開き方, plot方法) 天文データ解析入門 その4 (積分強度図の作り方)

SDF Editorのソースコードが公開された! プロトタイプとして公開されたSDF Editorの出力をPythonでGLSLにコンバートするコードを書いていたら、 なんと、5月30日にSDF Editorのソースコードが公開された。 https://joetech.itch.io/sdf-editor これによって、外部のPythonプログラムを使うことなく、直接GLSLコードを出力するようにソースコードを修正すればいいと考え、作業を開始したら一日ほどで動作するようになった。 これで、SDF Editorでモデリングし、プロジェクトをセーブしたら、<ファイル名>.fragを出力し、 それをVisualStudio Codeで読み込んだら、すぐに表示できるようになり、また、 htmlとしてセーブするとそのまま、WEBに置けば3DCGを表示できる。 SDF Editorでデザインし出力したコードをShadertoy.comで公開->https://www.shadertoy.com/view/Nts3W2 SDF EditorのソースコードはGODOT Engine SDF EditorのソースコードはGODOT Engineで書かれている。その本質はC++だが、Pythonライクなコーディング ができる。（ただし、Pyhtonとフルコンパチではない。）Pythonライクなので、今までPythonで書いたコンバータ がほぼそのまま使えて楽だろうと思ったら、件のPythonのコードはほとんど使わずに書けたので拍子抜けした。 今回は、改変したソースコードはデバッグ中なのでそのすべてを示さないが、いくつかポイントを紹介するので 誰でも独自にソースコードを変更することができるだろう。 大まかな方針は以下の通りである。 step1 SDFモデルや、tranformやrotateなどのmodifireの変数値をShadertoyのコードに取り込む step2 SDF関数をShadertoyのコードに取り込む 　SDF　EditorのソースコードのプロジェクトをGODOTに読み込み、 　シーン -> Saveのスクリプトをを表示し、以下のようにjson_saver.save(path)の下に emit_signal("my_signal",path)を追加する。 　func _on_file_selected(path): json_saver.save(path) emit_signal("my_signal",path) 次に、右端ペインのノードタブをタップし、my_signa()をマウス右ボタンクリックで、接続を選択し、 Scene -> Addに接続する。 sdf関数のための変数とコードを取得する また、SceneAdd.gdの先頭を次の二行をつけ加える。 var compiler = preload("res://SDF Scripts/Compilers/ShaderCompiler.gd").new() var uniforms = [] 変数の具体的数値を得る ここで、func _on_Save_my_signal(path): で、 var result = compiler.compile() var shdercode = result[0] uniforms = result[1] # sdf関数を書き込む var stemp = shdercode とすると、stempにuniformの変数名と、sdf関数の本体を文字列として取り込むことができる。 たとえば、球を一つ追加すると stempは下記のような内容になる。 uniform float SpRa_1; float sdf(vec3 p0) { float d1; d1 = pSphere(SpRa_1, p0); return d1; } ここで、SpRa_1を具体的数値に置き換え(球の半径）、このコードをShadertoyのテンプレートに読み込むと、 そのコードはShadertoyに投稿しても、VisiualStuido Code上からでも結果を表示できる。 変数の具体的数値は、以下のように取得できる。 var sdffnc = '' var j = 0 for u in uniforms: var value if u[0].has_method(u[1]): value = u[0].call(u[1]) else: value = u[0].get(u[1]) # var temp = lines[j].replace('uniform ','') このようにして取得した数値を変数名のところに置き換え、たとえば、 球を一個追加する例では、 float SpRa_1　＝1.0と書き換える。 float sdf(vec3 p0) { float d1; d1 = pSphere(SpRa_1, p0); return d1; } この方針で、すべての変数を具体的数値に変更し、Shadertoyのテンプレートコードに挿入して セーブすると、そのコードが実行し、表示できる。

ターミナルを開けると、ある日こんな警告が出ていました。 WARNING: `pyenv init -` no longer sets PATH. Run `pyenv init` to see the necessary changes to make to your configuration. 結論: https://github.com/pyenv/pyenv/blob/master/README.md に従って環境設定し直しましょう。 例えば Mac + ZSH + Homebrew の組み合わせの場合は以下の実行が必要です。 echo 'eval "$(pyenv init --path)"' >> ~/.zprofile うっかり警告どおり pyenv init を実行してしまうとこんな表示が出るのですが: # (The below instructions are intended for common # shell setups. See the README for more guidance # if they don't apply and/or don't work for you.) # Add pyenv executable to PATH and # enable shims by adding the following # to ~/.profile and ~/.zprofile: export PYENV_ROOT="$HOME/.pyenv" export PATH="$PYENV_ROOT/bin:$PATH" eval "$(pyenv init --path)" # Load pyenv into the shell by adding # the following to ~/.zshrc: eval "$(pyenv init -)" # Make sure to restart your entire logon session # for changes to profile files to take effect. あれ、私の環境には ~/.profile も ~/.zprofile も無いのにどうしたら良いんだ！と混乱する羽目になってしまうので、余計な事考えずに README.md を読むべきだという結論です。

解決すべき問題 $ git tag 0.1.0 こういうコマンドを実行すればGitレポジトリのなかでタグが作られる。Gitのタグとは何か？「サル先生のGit入門 タグ 」によればこうだ。 タグとは、コミットを参照しやすくするために、わかりやすい名前を付けるものです。 説明はこれで尽きる。でもタグを図に描いたらどうなるんだろう。複数のブランチを作ってそれぞれにたくさんタグを打ったら全体としてどう見えるだろう。 解決方法 Pythonでツール kazurayam/visualize_git_repository.py を開発した。これを使えばいま自分の手元にあるプロジェクトの .git ディレクトリのなかにあるオブジェクト群の実物を読み出し、Graphvizでグラフを生成してPNG画像ファイルを出力することができる。 説明 デモ用にディレクトリを作りGitレポジトリを作ろう。masterブランチとdevelopブランチでファイルを追加・変更などして、developをmasterにマージしよう。途中適宜にタグをつけていこう。Gitレポジトリの中がどのように形を変化させていくか、とくにタグがどんなふうに見えるか、図を示しつつ説明しよう。 ステップ1 最初のコミットをする 新しいGitレポジトリを作ろう。いくつかファイルを作ってコミットしよう。 作業用ディレクトリを作った git initした ファイルを３つ(README.md、.gitignore、src/greeting.pl)作った。 git add .した git commit -m "initial commit'した この時点ではまだタグがありません。この時点でvisualize_git_repositoryツールを実行したら次のグラフが生成された。 master develop まだ無い ステップ2 最初のコミットにタグ 0.1.0 を打つ タグ 0.1.0 を打ちましょう。 % git tag 0.1.0 コマンドを実行した。 この時点でvisualize_git_repositoryツールを実行したら次のグラフが生成された。 master develop まだ無い このグラフから次のことが読みとれる。 git tag XXXXコマンドはgit tag XXXX HEADと同じ意味だ。つまり特殊な参照名 HEAD が指すところのcommitオブジェクトにたいしてXXXXというタグを付ける。HEADが指すcommitオブジェクトとはつまり最後に実行したgit commit -m "....."コマンドによって作られた最新のcommitオブジェクトだ。 ステップ3 developブランチでタグ 0.2.0 を打つ 新しくdevelopブランチを作ろう。developブランチでファイルdoc/TODO.txtを追加してコミットしよう。そしてタグ 0.2.0 を打とう。 git branch --show-currentコマンドでいまカレントブランチがmasterであることを確認してから git branch develop を実行した ファイル doc/TODO.txt を追加した git add . を実行した git commit -m "add doc" を実行した git tag 0.2.0 を実行した この時点でvisualize_git_repositoryツールを実行したら次のグラフが生成された。 master develop 前と同じ このグラフから次のことが読みとれる。 いまカレント・ブランチがdevelopなのだが、masterブランチにおいて打ったタグ 0.1.0 が見えている。えっ！そうなんですか？「タグ 0.1.0 はmasterブランチで打ったタグなのだから、developブランチではタグ 0.1.0 が見えなくなるだろう」とわたしは思っていたが、そうじゃないんだ。 タグが個々のcommitオブジェクトに付けられた名札のようなものという基本に戻ろう。タグを作る時にカレント・ブランチが何だったかはタグの性質に影響しない。この絵をみてわたしははじめて納得しました。 ステップ4 masterブランチでタグ 0.1.1 を打つ masterブランチに戻って、masterブランチでもタグをひとつ打ってみよう。どう 見えるか？ git checkout master を実行した ファイル src/greeting.pl を修正した git add . を実行した git commit -m "modify src/greeting.pl" を実行した git tag 0.1.1 を実行した この時点でvisualize_git_repositoryツールを実行したら次のグラフが生成された。 master develop 前と同じ このグラフから次のことが読みとれる。 さきほどdevelopブランチでタグ 0.2.0 を打ったが、masterブランチではそれが見えない。 developブランチをまだmasterにマージしていないのだから、masterからdevelopの変更点が見えないのは当然だ。 ステップ5 developブランチをmasterにマージする developブランチをmasterにマージしよう。 git merge develop を実行した この時点でvisualize_git_repositoryツールを実行したら次のグラフが生成された。 master develop 前と同じ このグラフから次のことが読みとれる。 タグ 0.2.0 が付いたcommitオブジェクトがマージされて、masterブランチにおいて見えるようになった。 タグ 0.2.0 が付いたコミットとタグ 0.1.1 が付いたコミットがグラフのなかで上下逆順に並んでいるのを見ておや？と迷うかもしれない。気にしないでほしい。タグの文字つまり0.2.0と0.1.1を数として読んだときにわれわれが感じる意味的順序どおりにcommitノードをGraphvizが配置するようにツールを作ってはいない。必要なことではないと思うから。 ステップ6 masterブランチでタグ 0.2.1 を打つ さいごに、masterブランチでREADME.mdファイルを修正してコミットし、タグ 0.2.1 を打とう。 ファイル README.md を修正した git add . を実行した git commit -m "modified README.md" を実行した git tag 0.2.1 を実行した この時点でvisualize_git_repositoryツールを実行したら次のグラフが生成された。 master develop 前と同じ はい。もう付け足すことはありません。 ツールについて 本稿で示したPNG画像は自作のツール visualize_git_repository で描画した。このツールはPython言語で開発した。ソースコードは下記のGitHubレポジトリにある。 https://github.com/kazurayam/visualizing-git-repository このツールは下記2つのライブラリを利用している。 pytest python graphviz PNG画像を生成するにはコマンドラインで下記の操作をする。 $ cd $visualize_git_repository $ pytest -s kazurayam/visualize_git_repository.py::test_3_tags 上記の例を作るのにどういうgitコマンドを実行したのかを知りたいならプログラムのソースコードを読み解いてください。下記を入り口として解読してください。 kazurayam/visualize_git_repository.py まとめ 「タグとは、コミットを参照しやすくするために、わかりやすい名前を付けるものです」という説明をクリアに図示することができたとおもいます。 author: kazurayam date: June, 2021 連作の目次 第1回 commitとtreeとblob 第2回 ブランチとマージ 第3回 タグ

LINE Notify は、トークンを発行してコピペするだけでLINEと一瞬で連携できてとても簡単で簡単なので、PythonからLINEに通知するための爆速備忘録として掲載します。 LINE側 ログイン ページの右上のログインボタンから、メールアドレスとパスワードを入力してログイン。 マイページに移動 右上がログインから、LINEアカウントのユーザー名に変わっているので、クリックしてマイページを選択。 トークンを発行 一番下に「トークンを発行する」というボタンがあるのでクリック。 トークン名と送り先の設定 トークン名を適当に書いて、「1:1でLINE Notifyから通知を受け取る」を選択して「発行する」をクリック。 発行したトークンをコピー 発行したアクセストークンは、閉じると二度と見れなくなってしまうので、どこかメモ帳などで残しておくと良いですね。 Python側 下記のコードの、textに送りたいメッセージ内容と、tokenの中身を先ほど発行したアクセストークンに書き換えて実行します。 linenotify.py # ライブラリ import requests # 送りたい内容 text = 'こんにちは〜！' # アクセストークン -> LINE Notify ( https://notify-bot.line.me/ja/ ) のサイトからトークンを発行 token = '*******************************************' # 送信 url = "https://notify-api.line.me/api/notify" headers = {'Authorization': 'Bearer ' + token} payload = {'message': text} requests.post(url, headers=headers, params=payload) LINEに通知が来たら、成功です。

はじめに IoTが叫ばれている昨今（あれっ！？もうIoTって聞かなくなってきたような。死語かな…）、 手持ちのラズパイを使って何か出来ないかと考えていたところ、GoogleAPIを使うとスプレッドシートにデータをアップロード出来ることを発見。 備忘録として、それまでの手順をまとめます。 更に、Notionでグラフ表示できることも分かって、スゲー！ってなったので、それについても記載しています。 なので、基板周りはサラッと流させて下さい。 私としてはスプレッドシートにアップできたところ、Notionでグラフ表示できたところに感動があったため、その辺りにウェイトを置いた記事となります。 なお、言語はpythonで記載します。 ラズパイと温度センサの接続準備 ラズパイはRaspberry Pi 4Bを使用しています。 温度センサはアナログデバイセズ社のADT7310を使用しています。 秋月電子で売っているのを使いました（こちら)。 ラズパイとの接続はSPI通信方式となります。以下、参考接続図です。 （ADT7310がFritzingに見当たらなかったため別のセンサでSPI接続図を書いています。すみません。ホント、参考程度に） ADT7310の温度取得 SPI通信を使用するため、まずはラズパイの設定でSPI通信を使用許可して下さい。 raspi-configで設定しても良いですし、GUIであればスタートボタン→設定→Raspberry Piの設定からでも設定できます。 Raspberry Pi設定画面 pythonでSPI通信を使うために、spidevライブラリを使います。 インストールされていない場合はインストールしましょう（参考）。 ソースコード adt7310.py import spidev import time class ADT7310: def __init__(self) -> None: print('adt7310 init') self.spi = spidev.SpiDev() spi_ch = 0 # use cs0 pin self.spi.open(0, spi_ch) self.spi.mode = 0x03 self.spi.max_speed_hz = 1000000 time.sleep(0.25) self.spi.xfer2([0xff, 0xff, 0xff, 0xff]) def get_temp_adt7310(self): self.spi.xfer2([0x54]) # time.sleep(0.5) time.sleep(0.25) adc = self.spi.xfer2([0xff, 0xff]) temp = (adc[0] << 8) | adc[1] temp = temp >> 3 if(temp >= 4096): temp = temp - 8192 temp = temp / 16.0 return temp def close(self): self.spi.close() Googleスプレッドシートにアップするために さて、この記事のメインです。スプレッドシートに書き込みましょう！ 今回の方法ではGoogle Cloud PlatformにてGoogle Drive APIとGoogle Sheet APIを有効化し、サービスアカウント作成とjsonキーの取得を行います。 やっていきましょう。 Google Cloud Platformの設定1 まずはGoogle Cloud Platformにアクセス。 以下の画面になるのでそのまま「続行」しましょう。 何もせず、「認証情報に進む」ボタンをクリックします。 次の画面では、完了せずに「キャンセル」をクリックします。 APIの有効化 キャンセルした後の画面で、左メニューから「ダッシュボード」をクリック、さらに「＋APIとサービスを有効化」をクリックします。 APIライブラリでGoogle Drive APIを検索し、有効化します。 次は、Google Sheet APIを検索し、同様に有効化します。 サービスアカウントの作成とキー取得 APIとサービス画面から、「認証情報」をクリックし、「サービスアカウントの管理」をクリックします。 次の画面で「+サービスアカウントの作成」をクリックします。 サービスアカウントの作成画面に遷移するので、サービスアカウント名を記載して作成をクリックします。 今回はサービスアカウント名を「raspi_temp_sheet」としました。 サービスアカウントは自動で作成されますので特に記載する必要はないです。 作成後、②と③は省略可とあるため、そのまま完了してしまいましょう。 作成が完了したら、作成したアカウントのメールをクリックします。 「キー」タブをクリック後、鍵を追加をクリック。「新しい鍵を作成」を選択します。 鍵はjsonにしましょう。 jsonファイルがローカルにダウンロードされたと思います！ 後で使用しますのでなくさないようにしましょう。 スプレッドシートの作成と設定 ラズパイから書き込むためのスプレッドシートを作成します。 作成したら、ファイル名、シート名は各自で記載しましょう。私はシート名は「temp」にしました（あとで使います）。 また、1行目に 「date」「cpu temp」「atmos temp」の3つを記載しました。 dateはラズパイで温度をセンシングした時間 cpu tempはラズパイのCPUの温度 atmos tempは先ほどラズパイで接続したADT7310でセンシングした温度 とします。 このスプレッドシートを共有します。 共有相手は、先ほど作成したサービスアカウントのID（メールアドレス）となります。 これにてサービスアカウントとスプレッドシートとの連携が取れたかと思います。 次はラズパイからこのスプレッドシートにデータを書き込んでいきましょう。 ラズパイの設定とコード まずはライブラリをインストールしましょう。 pip install --upgrade google-auth google-api-python-client google-auth-httplib2 google-auth-oauthlib その後、以下のコードにて時刻・CPU温度・センサ温度をスプレッドシートに書き込みます。 先ほどのjsonファイルはpythonファイルと同じ階層に置いています。 各自、jsonファイル名や階層などによってソースコードを変更してください。 また、SPREAD_SHEET_IDは先ほど作成したスプレッドシートのURLの/d/の後から/editの前のURLとなります。 例えば、https://docs.google.com/spreadsheets/d/sample/edit#gid=0というURLの場合は sampleがSPREAD_SHEET_IDに該当します。 なお、ADT7310のソースコードも以下のpythonファイルの直下に置いてください。 main.py import datetime import sys import time from googleapiclient.discovery import build from googleapiclient.errors import HttpError as googleHttpError from google.oauth2 import service_account from gpiozero import CPUTemperature from adt7310 import ADT7310 # your create json file JSON_FILE = './key.json' # edit your spread sheet id SPREAD_SHEET_ID = 'your_sheet_id' # sheet name SHEET_NAME = 'temp' service = None def login_service(json_file, sheet_id): credentials = service_account.Credentials.from_service_account_file(json_file) scoped_credentials = credentials.with_scopes( [ 'https://spreadsheets.google.com/feeds', 'https://www.googleapis.com/auth/drive' ]) print('credentials: ', scoped_credentials) service = build('sheets', 'v4', credentials=scoped_credentials) sheet = service.spreadsheets() try: result = sheet.get(spreadsheetId=sheet_id).execute() except googleHttpError as e: print(e) for sheet in result['sheets']: print('Index %s: Sheet name = %s' % ( sheet['properties']['index'], sheet['properties']['title']) ) return service def append_sheet(service, sheet_id, sheet_name, val): values = [ val ] body = { 'values': values } ret = service.spreadsheets().values().append(spreadsheetId=sheet_id, range=sheet_name, valueInputOption='USER_ENTERED', body=body).execute() print('{0} cells appended.'.format(ret.get('updates').get('updatedCells'))) if __name__ == '__main__': my_adt7310 = ADT7310() time.sleep(0.1) while True: try: # Login if necessary. if service is None: service = login_service(JSON_FILE, SPREAD_SHEET_ID) # get cpu temp at gpiozero tmp = CPUTemperature() cpu_temp_gpio = tmp.temperature # get atmos temp atmos_temp = my_adt7310.get_temp_adt7310() # get date now_date = datetime.datetime.now().isoformat(' ') # upload to spreadsheet append_sheet(service, SPREAD_SHEET_ID, SHEET_NAME, [now_date, cpu_temp_gpio, atmos_temp]) print("{0}__cpu_temp: {1}'C, atmos_temp: {2}'C".format( now_date, cpu_temp_gpio, atmos_temp)) except KeyboardInterrupt: my_adt7310.close() sys.exit(0) except: # Error appending data, most likely because credentials are stale. # Null out the worksheet so a login is performed at the top of the loop. print('Append error, logging in again') worksheet = None # time.sleep(1) continue finally: # wait... # time.sleep(0.1/2) time.sleep(60) スプレッドシートに書き込めましたか？ 書き込めましたね！ やった！ Notionにシートのグラフを表示する ここまで結構な道のりでした。お疲れ様です。 最後にNotionで先ほどのスプレッドシートの情報をグラフ表示してみましょう。 Googleスプレッドシートの設定 共有設定にて、「リンクを知る人全員が閲覧可能」の設定に変更します。 Notionに表示する NotionにはNotion Chartsというサービスがあり、こちらを利用してNotion上にグラフ表示を行います。 まずは以下をクリックしてください。 こちら さて、一つずつ記載していきましょう。 GOOGLE SHEETS DOCUMENT ID : pythonファイルにも記載した、自身のスプレッドシートのリンクの一部です SHEET NAME : スプレッドシートのシート名（私の場合はtemp） DATA RANGE : グラフ表示したいセルの領域 ここでの変更は、「Craft your pallet.」の「Color2」を追加したのみです。 今回のグラフはCPU温度と気温の2つを表示したいため、Color2を追加しました。 最後に、作成です（Make Magic）！ 作成されたリンクをNotionに「Create embed」で貼り付けましょう。 すると、、、 表示が出来ましたね！ さらに、DATA RANGEを広めに指定しておけば、Notionのグラフが自動で更新されていくのがわかると思います。 まとめ ラズパイで得たデータをGoogleスプレッドシートに書き込んで、Notionでグラフ表示しました。 IoTやった感がありますね。 また、Google APIはスプレッドシート以外にもたくさんのサービスを利用できます。 気が向いたら他のサービスも触ってみたいと思いますし、皆さんもぜひ色んなサービスを使ってみてはいかがでしょうか。 最後に、Notionを使用し始めて3週間程度ですが、Notion凄い便利ですね。 凄い。スプレッドシートのグラフ表示が出来ることがわかり、より一層Notionが好きになりました。 参考リンク Google Sheet API python quickstart read write cell

※画像が多いです 今回はこれと これを組み合わせて これを作ります。(日本の形が変なのは許してください。) python、pillow、opencv、そしてnumpyを使います。 def mixTest(): material = cv2.imread("./resources/weather_material_FF0000.png") bgrLower = np.array([0, 0, 255]) bgrUpper = np.array([0, 0, 255]) # 抜き出す色は白に、それ以外は黒になる。 MonoMask = cv2.inRange(material, bgrLower, bgrUpper) cv2.imwrite("./resources/mask.png", MonoMask) orImg = Image.open("./resources/weather_report_for_test.png").convert("RGBA") mrImg = Image.open("./resources/weather_material_FF0000.png").convert("RGBA") msImg = Image.open("./resources/mask.png").convert("L") im = Image.composite(orImg, mrImg, msImg) im.show()

はじめに ワルラス法則とミクロ主体予算制約式1を用いて、インフレ目標達成への検討を行なった。 レポジトリーはこちらです。 ワルラス法則とミクロ主体予算制約式1 P1(t)*(y(t) - c(t)) + P2(t)*(b(t-1) - b(t)) +P2(t)*(bc(t-1) - bc(t)) + W(t)*(ls(t) - Id(t)) = MD(t) - MS(t) P1:財価格   y:財供給   c:財消費   P2:債券価格   b:家計保有債券   bc:政府保有債券   W:賃金   ls:労働供給   ld:労働需給   マネー需要:MD   MS:マネー供給 上式の意味は以下の通りである。 マネーの超過供給量は、通貨発酵益を経由して、財の超過需要を生み出している。1 本論では、以下の政策を実施する。 政府政策 ・インフレ目標2% ・GDP600兆円を目指す 日銀政策 ・フォワードガイダンス ・量的緩和(毎年80兆円の国債買い入れ、内民間と統合政府それぞれの保有残高増しは25兆円とする) また、上式の制約として、 一般の政府と中央銀行の統合政府を考える。ただし、簡単のために一般の政府は重要な役割を果たさず、統合政府は中央銀行とする。1 次回以降、 ・イールドカーブコントロール ・為替とGDPの関係(消費、投資、政府支出、輸出入) の検討とそれらを加味したモデルの検討を行いたい。 フィリップス曲線 インフレ率と失業率には強い逆相関(フィリップス曲線)があることが知られている。そこで、フィリップス曲線を用いて、インフレ率から失業率を導出する。 インポート import numpy as np import pandas as pd import matplotlib.pyplot as plt from google.colab import drive import time from sklearn.linear_model import LinearRegression drive.mount('/content/drive/') nov_dir = 'Colab Notebooks/dataset/phillips_curve/' work_dir = 'Colab Notebooks/workspace/export/' uer_path = '/content/drive/My Drive/' + nov_dir + 'UnemploymentRate.csv' cpi_path = '/content/drive/My Drive/' + nov_dir + 'corecoreCPI.csv' save_path = '/content/drive/My Drive/' + work_dir df_uer = pd.read_csv(uer_path, index_col='DATE', parse_dates=True) df_cpi = pd.read_csv(cpi_path, index_col='DATE', parse_dates=True) df_cpi = df_cpi[df_cpi.index != '2021-04-01'] print(df_uer.shape) print(df_cpi.shape) print(df_uer.head()) print(df_cpi.head()) Mounted at /content/drive/ (603, 1) (603, 1) UnemploymentRate DATE 1971-01-01 1.1 1971-02-01 1.2 1971-03-01 1.2 1971-04-01 1.2 1971-05-01 1.2 corecoreCPI DATE 1971-01-01 6.5 1971-02-01 6.2 1971-03-01 6.2 1971-04-01 6.8 1971-05-01 7.3 描写 plt.scatter(df_uer['UnemploymentRate'], df_cpi['corecoreCPI']) plt.title('Unemployment Rate & corecoreCPI') plt.xlabel('Unemployment Rate(%)') plt.ylabel('Core Core CPI(%)') plt.grid() plt.savefig(save_path+'PhillipsCurvePlot.png') plt.show() 教育 model = LinearRegression() X = df_uer[['UnemploymentRate']].values Y = df_cpi['corecoreCPI'].values model.fit(X, Y) a=model.intercept_ b=model.coef_ r2=model.score(X,Y) p_y=a/X+b 描写と保存 plt.scatter(X, Y) plt.plot(X, p_y, color='red', label="Phillips Curve") plt.title('Phillips Curve') plt.xlabel('Unemployment Rate(%)') plt.ylabel('Core Core CPI(%)') plt.grid() plt.text(3, 8, 'y = %s x %s' % (round(a,3),round(b[0],3))) plt.text(3, 6, 'R^2 = %s' % round(r2,3)) plt.legend() plt.savefig(save_path+'PhillipsCurve.png') plt.show() フィリップス曲線の結果より、 (インフレ率,失業率) = (2,2.567) となることがわかった。 乗数効果 内閣府の調査によれば、以下の乗数があることがわかった。(ともに2015年)2 名称 乗数 公共投資 1.14 減税 0.30 以上より以下のことがわかる。 公共投資の乗数効果は、短期目的には減税のそれよりも大きい。2 また、乗数効果の式は以下の式である。 \Delta Y = \alpha \Delta G \\ \Delta Y:国民所得,\alpha:乗数,G:政府支出 ここで、政府支出は一般会計の国費と限定する。 令和3年の12ヶ月会計 名称 額面(兆円) 国費 23.8 信用乗数(貨幣乗数)の推定 一般にマネタリーベースとマネーストックとの関係がある。 \Delta M = m \Delta H \\ \Delta M: マネタリーベース,m: 信用乗数,\Delta H: マネーストック しかし、景気循環にタイムラグがあり、定数が一定でないことが知られている。 そこで、関係性を調査し、可能であれば、推定を行いたい。 マネタリーベースのデータ加工 df = pd.read_excel(mbase_path, sheet_name=excel_sheet_name) df = df.drop(range(6)) df = pd.concat([df['Unnamed: 1'], df['Unnamed: 7']], axis=1) df = df.rename(columns={'Unnamed: 1': 'DATE', 'Unnamed: 7': 'MONETTALYBASE'}) df['DATE'] = pd.to_datetime(df['DATE'], format='%Y%m%d') df = df.reset_index(drop=True) drop_index = df.index[(df.index >= 612)] df = df.drop(drop_index) drop_index = df.index[(df.index <= 551)] df = df.drop(drop_index) df = df.set_index('DATE') df.index = df.index.strftime('%Y/%m') df.tail() マネーストックのデータ加工 df2 = pd.read_csv(mstock_path,encoding='cp932') df2 = df2[df2.index != 0] df2 = df2.rename(columns={'データコード': 'DATE','MD02\'MAM1NAM2M2MO': 'MONEYSTOCK'}) df2 = df2.drop(['MD02\'MAM1NAM3CCMO','MD02\'MAM1NAM3DMMO'], axis=1) df2 = df2.reset_index(drop=True) drop_index = df2.index[(df2.index >= 60)] df2 = df2.drop(drop_index) df2 = df2.set_index('DATE') df2.tail() 信用乗数 df3 = pd.concat([df['MONETTALYBASE'], df2['MONEYSTOCK']], axis=1) df3['MONEYSTOCK'] = df3['MONEYSTOCK'].astype(int) df3['MoneyMultiplier'] = df3['MONEYSTOCK']/df3['MONETTALYBASE'] df3 = df3.drop(['MONEYSTOCK','MONETTALYBASE'], axis=1) serial_num = pd.RangeIndex(start=1, stop=len(df.index) + 1, step=1) df3['No'] = serial_num df3.tail() 描画 plt.plot(df3.index,df3['MoneyMultiplier'], label="Plot of Days vs MoneyMultiplier") plt.title('Plot of Days vs MoneyMultiplier') plt.ylabel('MoneyMultiplier [-]') plt.xlabel('Days [Days]') plt.xticks(['2016/01','2017/01', '2018/01','2019/01', '2020/01']) plt.legend() plt.show() 教育 model = LinearRegression() X = df3[['No']].values Y = df3['MoneyMultiplier'].values model.fit(X, Y) a=model.coef_[0] b=model.intercept_ r2=model.score(X,Y) print('coefficient = ', model.coef_[0]) print('intercept = ', model.intercept_) coefficient = -0.008187783320053916 intercept = 2.3618821612964167 描画と保存 fig = plt.figure() plt.plot(df3.index,df3['MoneyMultiplier'], color = 'blue', label="Plot of Days vs MoneyMultiplier") plt.plot(X, model.predict(X), color = 'red', label="Estimated Money multiplier") plt.title('Estimated Money multiplier') plt.ylabel('MoneyMultiplier [-]') plt.xlabel('Days [Days]') plt.xticks(['2016/01','2017/01', '2018/01','2019/01', '2020/01']) plt.text(24, 2.3, 'MoneyMultiplier = MONEYSTOCK / \n MONETTALYBASE') plt.text(24, 2.25, 'y = %s x + %s (x: Days=>times)' % (round(a,3),round(b,3))) plt.text(24, 2.2, 'R^2 = %s' % round(r2,3)) plt.legend() fig.savefig(save_path) plt.show() 以上より、信用乗数を2.0と推定する。 マネタリーベースの平均変化率 マネタリーベースの定義は以下の通りである。 マネタリーベース＝「日本銀行券発行高」＋「貨幣流通高」＋「日銀当座預金」3 EXCELデータの読み込み加工 df = pd.read_excel(mbase_path, sheet_name=excel_sheet_name) df = df.drop(range(6)) df = pd.concat([df['Unnamed: 1'], df['Unnamed: 7']], axis=1) df = df.rename(columns={'Unnamed: 1': 'DATE', 'Unnamed: 7': 'MONETTALYBASE'}) df['DATE'] = pd.to_datetime(df['DATE'], format='%Y%m%d') df = df.reset_index(drop=True) drop_index = df.index[(df.index >= 612)] df = df.drop(drop_index) df = df.set_index('DATE') df.head() 年の増加率 df2 = df[(df.index.month == 1)] df2.index = df2.index.year df2 = df2.rename(columns={'MONETTALYBASE':'JAN'}) df3 = df[(df.index.month == 12)] df3.index = df3.index.year df3 = df3.rename(columns={'MONETTALYBASE':'DEC'}) dfy = pd.concat([df2['JAN'], df3['DEC']], axis=1) dfy['MONETTALYBASE'] = dfy['DEC'] - dfy['JAN'] dfy = dfy.drop(['DEC', 'JAN'], axis=1) dfy.to_csv(save_path, sep=",") dfy.tail() 描画 fig = plt.figure() plt.plot(dfy.index,dfy['MONETTALYBASE'], label="Plot of Years vs Average rate of change") plt.title('Plot of Days vs MoneyMultiplier') plt.ylabel('Average rate of change [-]') plt.xlabel('Years [Years]') plt.legend() fig.savefig(save_path4) plt.show() 期待GDPとマネタリーベースの算出 GDP=23.8*1.14 print('GDP =' + str(round(GDP,2))) now_GDP=df2['REAL'][0] print(now_GDP) for_GDP=now_GDP+GDP print('for_GDP =' + str(for_GDP)) all_wage = df2['REAL'][2] * df2['REAL'][3]*0.000000000001 print('all_wage =' + str(all_wage)) for_mometqary_base=(50+all_wage+GDP)/2 print('for_mometqary_base ='+str(round(for_mometqary_base,2))) now_GDP =27.13 535.7 for_GDP =562.832 all_wage =9.3075 for_mometqary_base =43.22 以上より、 完全雇用(インフレ率2%)を達成した時の 予想名目GDPは、562兆円 予想マネタリーベースは43兆円 と妥当性が高く、可能性は高い。 GDP600兆円は達成できないと思われる。 岩田喜久男 編、まずデフレを止めよ、日本経済新聞社、2003年 ↩ 社会資本整備等の現状 ↩ 「マネタリーベース」とは何ですか？ ↩

Pythonのコマンドラインツール作成方法 Pythonでは、標準ライブラリーのargparseを使って、コマンドラインツールを作ることができます。 また、サードパーティのライブラリーが、Awesome Python - Command-line Interface Development に紹介されています。 この中でGitHub上でStar数が最も多いのはPython Fireです。 Python Fireについて Python Fireは、Google製のライブラリーです。pip install fireでインストールすると使えます。 Python Fireを使うと、「コマンドラインからオブジェクトを自由に操作」できるようになります。 そのために必要なことは、オブジェクトをFireに登録することだけです。簡単に使い始められるというのがメリットです。 オブジェクトには、関数やクラスも含めて任意のオブジェクトが扱えます。 ドキュメントは、The Python Fire Guide です。 以下に主な使用方法を紹介します。 関数のコマンドライン化 Python Fireを使うには、from fire import Fireして、Fire(コマンドラインで使用したい関数)とします。つまり「使用したい関数をFireに登録」だけです。 下記をsample.pyというファイル名で作成してください（以降も同名のファイルを扱います）。 from random import choice from fire import Fire Fire(choice) python sample.py [1,2,3]を実行すると、1, 2, 3のいずれかが出力されます。これは、Pythonでprint(choice([1, 2, 3]))を実行するのと同じになります。 また、python sample.py [大吉,吉,小吉]を実行すると、大吉, 吉, 小吉のいずれかが出力されます。 python sample.py ["大吉","吉","小吉"]としなくても「大吉」などを文字列と認識します。 python sample.py [大吉, 吉, 小吉]のようにスペースを入れるとエラーになるので注意してください。 Fire(choice)とするだけで、簡単に、おみくじとして使うことができます。 ヘルプの表示 以下のいずれの方法でもヘルプを表示できます。 python sample.py --help python sample.py -- --help python sample.py -h python sample.py -- -h このように、Python Fireでは、自動的に「ヘルプを表示する機能」が使えます。 --helpあるいは-hをつけることでヘルプが表示されます。 「Choose a random element from a non-empty sequence.」という説明が出てきます。 --を書くことによって、以降の引数が、sample.pyに対してはなくFireに対するものになります。今回は、「--helpや-h」がsample.pyの引数でないことが自明なので、--は不要になっています。 また、自作の関数を登録して、詳しいヘルプを出したい場合は、その関数にdocstringを記述してください。 複数の関数のコマンドライン化 複数の関数を使えるようにするには、Fireに辞書を渡します。 辞書のキーにコマンド名を、値に関数を指定します。 たとえば、Fire({"min": min, "max": max})とすることで、minとmaxの2つのコマンドが使えます。 python sample.py -hまたはpython sample.pyでヘルプを表示します。 ヘルプを見ると、python sample.py コマンドのように使えて、コマンドとしてminとmaxが選べることがわかります。 python sample.py min 1 -2 3を実行すると、min(1, -2, 3)を計算して-2と表示されます。 python sample.py max 1 -2 3を実行すると、max(1, -2, 3)を計算して3と表示されます。 引数にリストを指定できます。 python sample.py min [1,-2,3]を実行すると、min([1, -2, 3])を計算して-2と表示されます。 引数に文字列を指定できます。 python sample.py min catを実行すると、min("cat")を計算してaと表示されます。min("cat")は、min(["c", "a", "t"])と同じです。 辞書のキーがコマンドになります。 もし、Fire({"MIN": min, "MAX": max})としていれば、python sample.py MIN 1 -2 3のように使います。 辞書の要素数は、3つ以上でも可能です。 もし、Fire({"min": min, "max": max, "sum": sum})としていれば、python sample.py sum [1,-2,3]を実行すると、sum([1, -2, 3])を計算して2と表示されます。 引数の文字列はPythonのコードのように評価され、適切な型のオブジェクトになって実行されます。 python sample.py min [1,-2,3]を実行すると、min関数の引数はリストになります。 python sample.py min catを実行すると、min関数の引数は文字列になります。 Fire()のように空にした場合 Fire()のように引数を指定しないと、そのモジュールに出てくる変数をすべて登録します。 from builtins import min, max from fire import Fire as _Fire _Fire() 上記のようにすると、Fire({"min": min, "max": max})と同じになります。なお、_Fireは、アンダースコアで始まっているため無視されます。 クラスのコマンドライン化 クラスを指定することもできます。 下記をsample.pyとします。 from fire import Fire class Command: min = min max = max Fire(Command) Commandクラスは、minとmaxのスタティックなメソッドを持っています。 このとき、Fire(Command)とすることで、下記と同じように動作します。 c = Command() Fire({"min": c.min, "max": c.max}) クラスを登録した場合でも、使い方は複数の関数の場合と変わりません。 実行例 結果 python sample.py min cat a python sample.py max cat t python sample.py ヘルプを表示 python sample.py min -- --help minのヘルプを表示 python sample.py min catを実行すると、Command().min('cat')を実行し結果を表示します。 メソッドを定義したクラスのコマンドライン化 スタティックでないメソッドを定義したクラスでは、以下のようになります。 from fire import Fire class Command: def min(self, *args): """return smallest args""" return min(args) def max(self, *args): """return largest args""" return max(args) Fire(Command) Commandクラスのメソッドは、スタティックでないので、第1引数は自分自身（self）にします。 実行例は以下のようになります。 実行例 結果 python sample.py min 3 -2 1 -2 python sample.py max 3 -2 1 3 ※ python sample.py min 3 -2 1は、Command().min(3, -2, 1)に対応し、min((3, -2, 1))を計算し、-2と出力されます。 ※ python sample.py min catは、Command().min("cat")に対応し、min(("cat", ))を計算し、catと出力されることに注意してください。 __init__メソッドを定義したクラスのコマンドライン化 前節では、python sample.py minを実行すると、Command().min()が呼ばれ、min(())を実行しようとしてエラーになります。 一般に、min((), default=0)のようにdefaultを指定すると、エラーにならずに0を返すようになります。 ここでは、__init__メソッドを作成し、minのdefault引数を指定できるようにします。 以下のように、Commandクラスの__init__メソッドでdefaultを受け取れるようにし、それをminやmaxで使うようにします。 from fire import Fire class Command: def __init__(self, default=0): self.default = default def min(self, *args): """return smallest args""" return min(args, default=self.default) def max(self, *args): """return largest args""" return max(args, default=self.default) Fire(Command()) このようにすることで、引数が空でもエラーにならずに、0を返します。 引数defaultに値を指定するには、下記のように--default 値というオプションをつけます。 実行例 結果 python sample.py min 0 python sample.py max --default -1 -1 別々の引数があるクラスのコマンドライン化 __init__メソッドと、通常のメソッド（some）の両方に引数がある場合は以下のようになります。 from fire import Fire class Command: def __init__(self, arg1=1): self.arg1 = arg1 def some(self, arg2=2): return self.arg1, arg2 Fire(Command) 実行例 結果 python sample.py some [1, 2] python sample.py some --arg1 10 [10, 2] python sample.py some 20 [1, 20] python sample.py some --arg2 20 [1, 20] python sample.py some --arg1 10 20 [10, 20] python sample.py some --arg1 10 --arg2 20 [10, 20] メソッドの引数にデフォルト値があれば、コマンドラインで省略できます。 __init__メソッドと通常のメソッドのどちらの引数も、同じように指定できます。同名の引数の場合、__init__メソッドが優先されます。 --引数名をつけない場合、通常のメソッドの引数と解釈されます。 Python Fireへのオプションを指定したい場合は、python sample.py -- --help のように--の後に書きます。 補足 もし、def some(self, arg2):のように定義していれば、下記はarg2が指定されていないのでエラーになります。 python sample.py some python sample.py some --arg1 10 サブコマンドについて Python Fireでは、サブコマンドを持つツールも簡単に作成できます。 from fire import Fire class Bin: """Convert binary number from/to decimal number.""" def from_dec(self, n): """Convert binary number from decimal number.""" return bin(n) def to_dec(self, n): """Convert binary number to decimal number.""" return int(str(n), 2) class Oct: """Convert octal number from/to decimal number.""" def from_dec(self, n): """Convert octal number from decimal number.""" return oct(n) def to_dec(self, n): """Convert octal number to decimal number.""" return int(str(n), 8) class Command: bin = Bin oct = Oct Fire(Command) 上記のように記述することで、binコマンドとoctコマンドが使えます。 binコマンドの値は、Binなので、Binクラスのメソッドをサブコマンドとして使えます。 octコマンドの値は、Octなので、Octクラスのメソッドをサブコマンドとして使えます。 実行例 結果 意味 python sample.py （略） コマンドのヘルプ表示 python sample.py bin （略） binコマンドのサブコマンドのヘルプ表示 python sample.py bin from_dec 3 0b11 10進数の数字3を2進数の数字に変換 python sample.py bin to_dec 11 3 2進数の数字11を10進数の数字に変換 python sample.py oct from_dec 9 0o11 10進数の数字9を8進数の数字に変換 python sample.py oct to_dec 11 9 8進数の数字11を10進数の数字に変換 補足その1 class Command: bin = Bin oct = Oct Fire(Command) 上記の代わりに、Fire({"bin": Bin, "oct": Oct})のように辞書で書いても同じ機能になります。 ただし、Commandクラスにはdocstringを書けますが、辞書では書けません。 補足その2 class Command: bin = Bin() oct = Oct() Fire(Command) 上記のように、「クラスではなくオブジェクトを指定」しても同じように使えます。 コマンドチェーンについて Python Fireでは、コマンドをつなげていくことができます。 python sample.py xxxが返すオブジェクトがyyyという（メソッドなどの）属性を持っていれば、python sample.py xxx yyyを実行できます。 また、それが返すオブジェクトがzzzという属性を持っていれば、python sample.py xxx yyy zzzを実行できます。 これは、いくつでもつなげることができます。 下記をsample.pyとします。 import calendar from fire import Fire Fire(calendar) python sample.py TextCalendar prmonth 2021 2を実行すると下記のように、2021年2月のカレンダーが表示されます。 February 2021 Mo Tu We Th Fr Sa Su 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 Fireでは、calendarモジュールを登録しています。 calendarモジュールでは、TextCalendarクラスが使えます。 TextCalendarクラスでは、prmonthというカレンダーを表示するメソッドが使えます。 prmonthメソッドの引数は、年と月です。 以上から、python sample.py TextCalendar prmonth 2021 2を実行すると、calendar.TextCalendar().prmonth(2021, 2)に対応し、カレンダーが表示されます。 補足 Fireは、オプションで実行状態のトレースを出力できます。デバッグ時に役立ちます。 python sample.py TextCalendar prmonth 2021 2 -- --traceとすると下記のように出力されます。 Fire trace: 1. Initial component 2. Accessed property "TextCalendar" (中略/calendar.py:293) 3. Instantiated class "TextCalendar" (中略/calendar.py:293) 4. Accessed property "prmonth" (中略/calendar.py:346) 5. Called routine "prmonth" (中略/calendar.py:346) クラス指定、オブジェクト指定の違い クラスを指定する例 class Command: def __init__(self, hello='Hello'): self.hello = hello def say(self) return self.hello Fire(Command) # クラスを指定 python sample.py --hello Hi sayでHiと出力されます。このように、実行時に__init__メソッドの引数を指定できます。 オブジェクトを指定する例 class Command: def __init__(self, hello='Hello'): self.hello = hello def say(self): return self.hello Fire(Command('Hi')) # オブジェクトを指定 python sample.py sayでHiと出力されます。このように、実行時に指定していない__init__メソッドの引数が使われます。__init__メソッドの引数は変更できません。 クラス指定、オブジェクト指定の違いのまとめ クラスを指定 オブジェクトを指定 ヘルプ中の文言「sample.py COMMAND」 ヘルプ中の文言「sample.py GROUP」 オブジェクトを自動で生成 オブジェクトを明示的に生成 __init__メソッドの引数を指定可能 __init__メソッドの引数は指定不可 副作用のあるメソッドのコマンドチェーン ここでは、print()で結果を出力するメソッドをつなげる例を紹介します。 import random from fire import Fire class Command: def __init__(self, *words): self._words = list(words) def s(self): """shuffle""" random.shuffle(self._words) print(" ".join(self._words)) return self def r(self): """reverse""" self._words.reverse() print(" ".join(self._words)) return self def __str__(self): return "" Fire(Command) 実行例 python __N__.py A B C D E - s：A〜Eをシャッフルして出力 python __N__.py A B C D E - s s：A〜Eをシャッフルして出力を繰り返す python __N__.py A B C D E - s s r：A〜Eをシャッフルして出力を繰り返し、反転し出力 コードの説明 -引数は、現在のメソッド（__init__メソッド）の引数の終わりを指示します。 s()は、self._wordsをシャッフルして出力します。 r()は、self._wordsを反転して出力します。 s()やr()はselfを返すので、最後に__str__メソッドが呼ばれて空文字列を出力します。 _wordsという代わりにwordsという名前のデータ属性にすると、wordsというコマンドがヘルプに表示されます。 _wordsのようにアンダースコアで始まる名前にすると、ヘルプには表示されませんが、利用はできます。 ブール型の引数 ブール型の引数は値を省略できます。 def show_args(args): print(args) Fire(show_args) 実行例 コマンド 出力 python sample.py --args True python sample.py --noargs False python sample.py --args True True python sample.py --args False False 引数の値を省略すると、Trueを指定したとみなされます。 引数の名前にnoをつけて、値を省略すると、Falseを指定したとみなされます。

今回やる事 競技プログラミングをやっていたら、再帰下降構文解析を応用したようなアルゴリズムの問題がけっこう出て解けなかったので、まずは基本から理解しようと思い、自分なりにまとめてみました。 ※間違っていたらご指摘いただけると喜びます。 環境 python 3.8.5 jupyter notebook 6.2.0 構文解析とは 文字列を解析して、その文字の構成要素の意味を解析する事。 自然言語における解析と、コンパイラやプログラムなどでの使用を目的とした解析では手法が異なる。 今回は、競技プログラミングなどのプログラムおける使用を前提としているので、後者の方を追っていきます。 再帰下降構文解析とは 再帰により構文を小さな単位に分割し、解析するもので、比較的簡単に実装出来て性能はそこそこいいため競技プログラミングなどに頻出する。 処理の流れ的には、入力文字列を左から構文解析していって末尾から値を確定させていく。 Wikipediaより引用 再帰下降構文解析は相互再帰型(※1)の手続き（あるいは再帰型でない同等の手続き）で構成されるようなLL法のトップダウン構文解析(※2)であり、各プロシージャが文法の各生成規則を実装する事が多い。従って、生成されるプログラムの構造はほぼ正確にその文法を反映したものとなる。そのような実装の構文解析器を再帰下降パーサと呼ぶ。 1. 相互再帰型 再帰の種類の1つで、それ自身をループさせるのではなく、複数の関数の間で相互に再帰になっているもの。 2. LL法のトップダウン構文解析 構文解析のアルゴリズムを大きく分けると、 トップダウン型（下降型解析） LL法 再帰下降構文解析 末尾再帰構文解析 パックラット構文解析 アップダウン型（上昇型解析） LR法 SLR法 LALR法 CLR法 GLR法 アーリー法 CYK法 などがあり、LL法というのは、トップダウン型の構文解析のアルゴリズムの1つ。 文脈自由文法というものの一種で、文字列を左（Left）から解析していき、左端導出（Leftmost Derivation）という左から非終端記号（他と置換が出来るもの）を終端記号（置換が出来ないもの）に置き換えていく手法を取る。 また、左から解析していき、右端導出を行うものはLR法という。 文脈自由文法については、こちらのサイトが分かりやすくまとめてありました。 うさぎでもわかるオートマトンと言語理論　第07羽　文脈自由文法 つまり、再帰下降構文解析も文脈自由文法に対する構文解析手法の一種みたいです。 どんな事が出来るのか？ 四則演算や、何らかの規則によって成り立っている文字列の構造を解析出来たりする。 このアルゴリズムの例として、四則演算がよく挙げられているが、競技プログラミングとかでは四則演算のような仕組みを応用し、生成した文字列を解析する問題がけっこう存在する。 四則演算で理解する 例として、2(2+3*(6+1)+5)といった簡単な四則演算を実装する事を考えます。 ここで重要になるのは、計算の優先順位の実装をする事で ①. 【括弧の中】 ②. 【 * 】 or 【 / 】 or 【 ( 】※ ③. 【 + 】 or 【 - 】 の順で計算させる事です。 括弧の中に括弧があったら一番内側を優先して計算させる必要もあります。 ※*が省略されている場合は【（ 】の次に数字が来る。 上記の実装はEBNFで表すとこうなります。 expr = term, {("+", term) | ("-", term)} term = factor, {("*", factor) | ("/", factor) | ("(", factor)} factor = ("(", expr, ")") | number number = 1つ以上の数字 EBNFはBNFの拡張版で、構文定義を書いたりするメタ言語です。ここではプログラムの仕様を定義します。 EBNFでは左辺が非終端記号、右辺が終端記号になっており、左辺を右辺に置き換えて構文を解析していきます。 つまり、exprで言うとexprはterm（項）と、+か-だったらterm（項）に置き換えられるという事みたいです。 2(3+4)+5という式を例にして大まかに考え方を見る ①. exprは項（+か-だったら、さらに項）に置き換えられる。 まず先頭から探索し、2をtermとして置き換える。 ②. termはfactor（* か / か ( だったら更にfactor）に置き換えられる 2をfactorとして置き換える。 ③. factorはexprかnumberに置き換えられる 2は括弧始まりではないため、number（1つ以上の数字）に置き換えられる。 このようなexpr→term→factorという再帰を繰り返していき計算を進めていきます。 EBNFを元に実装してみる test.py def countup(): global i i += 1 # number = 1つ以上の数字 def number(s): global i res = '' while i < len(s) and s[i].isdigit(): res += s[i] countup() return int(res) # expr = term, {("+", term) | ("-", term)} def expr(s): global i res = term(s) while i < len(s): if s[i] == '+': countup() res += term(s) continue elif s[i] == '-': countup() res -= term(s) continue return res return res # term = factor, {("*", factor) | ("/", factor) | ("(", factor)} def term(s): global i res = factor(s) while i < len(s): if s[i] == '*': countup() res *= factor(s) continue elif s[i] == '/': countup() res /= factor(s) continue elif s[i] == '(' and str(res).isdigit(): res *= factor(s) continue else: break return res # factor = ("(", expr, ")") | number def factor(s): global i if s[i] == '(': countup() res = expr(s) if s[i] == ')': countup() return res return number(s) i = 0 s = '2(3+4)+5' print(expr(s)) 実行結果 19 処理の流れを見る。 実装内容 expr termを呼び結果を取得。 + or - だったらtermを呼ぶ。 term factorを呼び結果を取得。 * or / or ( の次数字 だったらfactorを呼ぶ。 factor ( だったらexprを呼び結果を取得。 ( じゃなかったらnumberを呼ぶ。 処理の流れ ①. expr まずtermを呼ぶ。 ②. term factorを呼ぶ。 ③. factor 括弧始まりじゃないためnumberを呼ぶ。 ④. number 先頭から数字じゃなくなるまで検索し、2を返す。 ※これから検索した位置は保持しておく。 ⑤. term factorから2が返り、先を検索すると(でその次が数字のためfactorを呼ぶ。 　 ⑥. factor (は括弧始まりのため次に進めexprを呼び結果を取得。 ※括弧が終わるまで再帰的に探索される。 ⑦. expr termを呼ぶ。 ⑧. term factorを呼ぶ。 ⑨. factor 3なのでnumberによって3が返される。 ⑩. term factorから3が返されるが *, /, ( ではないので3を返す。 ⑪. expr 次は+なので次に進めtermを呼ぶ。 ⑫. term factorを呼ぶ。 ⑬. factor 4を返す。 ⑭. expr 3+4が処理され、次に進む。 ⑮. factor ⑥の結果が返り、括弧が終わると、⑤に結果が返り2*7が処理される。 その後再度exprが呼ばれる。 ⑯. expr +のため、次に進みtermを呼ぶ。 ⑰. term factorを呼ぶ。 ⑱. factor 5のためnumberを返す。 ⑲. expr numberから返った5と14が足し合わされ、末尾まで検索したため結果が返る。 感想 構文解析系のアルゴリズムは全く触れたことが無かったので、すごい勉強になりますね。 そもそもの基礎的な知識が全く足りてない事に気づいたのでそちらも勉強していきたいと思いました。 簡単な四則演算だけじゃなく、文字列が入った計算とかも実装してみます。 参考にさせていただいた記事 Java 再帰下降構文解析 超入門

Rやpythonで作図したグラフを論文で使用する際に、フォント周りで苦労しないように設定方法をメモしておく。 Python matplotlib 3.3.4 スタイルシートmatplotlibrcに以下の設定を加える .matplotlibrc font.family : sans-serif font.sans-serif : Helvetica #(or your favorite) pdf.fonttype : 42 # 同様にrcParamsを変更すれば自分のデフォルトテーマを設定可能 matplotlibの初期フォント設定は、"DejaVu Serif"でType 3のフォント形式となっている。フォントの形式の詳細に関しては勉強不足なので割愛。 スタイルシートの書き方は下記参考： https://note.nkmk.me/python-matplotlib-matplotlibrc-stylesheet/ R ggplot2 3.3.3 themeのbase_familyで指定する。 Rが起動する時に読み込まれる.Rprofileに好きなthemeを定義し、theme_set()で設定しておくと楽。 .Rprofile # dplyrのselectやfilterがマスクされないように後からtidyverseを読み込む options(defaultPackages = c(getOption('defaultPackages'), "MASS", 'tidyverse')) # 無いとtheme_setが見つからないとエラーが出た library(ggplot2) # マイテーマの記述　フォント以外の設定もお好みで my_theme <- function(TXTSIZE = 24, LEDSIZE = 18, LINEWIDTH = 2) { theme_bw(base_family="Helvetica") + theme( axis.text=element_text(color="black",size= TXTSIZE), panel.border = element_rect(color = "black", fill = NA, size = LINEWIDTH), legend.text = element_text(size= LEDSIZE) # ... ) } # マイテーマの設定 theme_set(my_theme()) Rprofileの書き方は下記参考： https://stats.biopapyrus.jp/r/devel/rprofile.html

Pythonのグローバル変数 Pythonを使っていると、モジュール間共通のグローバル変数を使いたくなるシーンがあります。 このグローバル変数がとても扱いにくく 関数内からグローバル変数を書き換える場合、 使用前にglobalで変数を宣言しなければいけません。 main.py global val1 def func1(): global val1 val1 = 100 def main(): global val1 func1() print(val1) 扱う変数の量が増えると、行頭のglobal宣言だけで数行埋まってしまいます。 プログラムをクラス化すれば良いのですが、Pythonではクラス変数の前にself.をつけなければいけないという制約があるため、コードが冗長になってしまいます。 また、モジュール外から他のモジュールのグローバル変数にアクセスしようとしても、思うようにいきません。 my_function.py def func1(): global val1 val1 = 100 main.py from my_function import * global val1 def main(): func1() print(val1) このように書いてもグローバル変数を別モジュールから書き換えることが出来ません。 共通のファイルをインポートしてas gとする これらの問題を解決する、とっておきの方法を思いつきました。 global_value.pyという中身空っぽのファイルを作成します(名前は何でも良いです) global_value.py #中身は空っぽ 使う側で、それをgとしてimportします そして、グローバル変数をg.val1のように使用します。 my_function.py import global_value as g def func1(): g.val1 = 1 main.py import my_function import global_value as g def main(): func1() print(g.val1) こうすることでglobalの宣言は一切必要なくなり g.をつけた変数はグローバル変数として扱えるようになります。 この変数は、他のモジュールからでも、クラス内からでも自由にアクセスできます。 変数名はやや冗長になりますが、ローカルとグローバルを明示する前置詞として働くため、プラスの効果も大きいと思います。 グローバル変数のメリット Pythonはスクリプト言語です。 スクリプト言語は遅いだの、レベルが低いだの難癖を付けられがちですが 他の言語にはない強力な機能が用意されています。 それが、文字列をコードとして扱えるeval関数とexec関数です。 exec("g.val=100") とすれば、g.valの値を100に書き換えることが出来ます。 これがtkinterと変数の同期処理で輝きます。 tkinterウィジェットを管理する辞書のキー値に変数名を用いることで決まったループ処理１つで変数とGUIの同期ができます。 my_gui.py widgets = dict() widgets["g.val1"] = widget1 widgets["g.val2"] = widget2 widgets["g.val3"] = widget3 def UpdateValue(): for key in widgets.keys(): exec(key+"="+str(widgets[key].data.get())) 実際はウィジェットクラスに適切な文字列を返すget_data_textメソッドを用意したほうが良いです。 この壁になっていたのが、グローバル変数問題でした。 どうにかしてウィジェットクラス内部から、メインプログラムのグローバル変数を書き換えられないかと悩んでいたら、この方法を思いつきました。 これが、誰かのお役に立てば幸いです。

Owaspのやられサイトで勉強してみた。 KBID 1 - Path traversal (LFI) やられサイトの構築 requirements.txtでパッケージを一括でインストールする。 $ pip3 install -r requirements.txt pythonでLFIをたちあげる。 $ python3 LFI.py Reconnaissance Intro, Chapter 1とChapter 2を選択するプルダウンメニューがある。 valueはintro.txtファイルから取得されている。 filenameパラメータはシステム内のファイルからコンテンツを取得するために利用されている。 filenameをユーザーが操作できる状態にあるため、想定外の機密情報を盗まれる可能性がある。 LFI.py @app.route("/home", methods=['POST']) def home(): filename = request.form['filename'] if filename == "": filename = "text/default.txt" f = open(filename,'r') read = f.read() return render_template("index.html",read = read) Exploitation text/default.txtを../../../../../etc/passwdに書き換える。 ペネトレーションテストではLFIの証明によく/etc/passwdを取得しようとする。らしい。 Submit buttonをクリックすると/etc/passwd情報が表示される。成功。

出発地 Python実践データ分析100本ノック の「ノック86:　画像内の人がどこに顔を向けているのかを検出してみよう」で用いられているライブラリ dlib のインストールで迷ったので記録しておきます。 環境によって対処方法が異なると思いますが、私の場合はこうでした。 環境 OS: windows 10.0 python: 3.7.3 Anaconda: conda 4.10.1 実施日 2021/06/09 利用ツール Anaconda Prompt を管理者として実行 実施コマンド conda install -y -c conda-forge dlib 到着地 dlib: 19.22.0 conda list # packages in environment at C:\ProgramData\Anaconda3: # # Name Version Build Channel （中略） dlib 19.22.0 py37h38cb4a7_0 conda-forge （後略）

私が最初に線形代数を学んだのはこの本でした。 マンガ　線形代数入門　はじめての人でも楽しく学べる しばらくはこの本に従って解説を進めていきたいと思います。 第一章「行列って何？　線形代数って何？」 まずは直線性(Linearity=グラフとして表した時に直線となるような数学的関係)とは何か説明する為、以下の中学校で習う公式から出発します。 ax+bx=c 確率の期待値(平均)の場合 E(aX+bY)=aE(X)+bE(Y) 微分の場合 (af(x)+bg(x))'=af'(x)+bg'(x) ax+bx=cの場合、代数方程式(Algebraic Formula)の表現(Expression)では以下となります。 【数理考古学】代数方程式について。 ax+bx-c=0 これを行列(Matrix)で表してみましょう。 行列の行とは横の事、列とは縦の事。合わせると縦と横に数字を綺麗に並べたものといえる。日本では数学用語Matrixを和訳する際にこれを「行と列の組み合わせ」と捉えたのでかかる表現が採用されたが、英語では窓口に並ぶ行列の事をQueueという。 マトリックス - Wikipedia 本来は「子宮」を意味するラテン語（< Mater母+ix）に由来するMatrixの音写で（英語では「メイトリクス」）、そこから何かを生み出すものを意味する。この「生み出す機能」に着目して命名されることが多い。また、子宮状の形状・状態に着目して命名される場合もある。 日本語にあえて翻訳する場合は「基盤」「基質」などの訳語が当てられることがあるが、原語で強く感じられる「生み出す機能」や「形状」が伝わりにくく、必ずしも評判が良くない。例えば「母体」あるいは「子宮体」ならばニュアンスも伝わるのだろうが、このような訳語はほとんど採用されていない。結局、カタカナで表記されることが多い。 生物学や医学では「間質」という言葉も使われている(例：細胞間質)。 材料学では「母材」、鉱物学では「母岩」という訳語も使われている(複合材料の母材)。 SFドラマ「ドクター・フー（1963年-）」88話「Deadly Assassin（1976年）」に登場する用語。知識が集積された仮想空間のことを「マトリックス」と呼んでいる。 聖悠紀の漫画作品「超人ロック」シリーズ(1967年~)において、超能力者が他人に変身するときに、他人の身体からコピーして自分の体を作り変える遺伝子情報を遺伝子マトリクスと呼んでいる。 ウィリアム・ギブスンのサイバーパンクSF小説「ニューロマンサー（1984年）」「カウント・ゼロ（1986年）」「モナリザ・オーバードライブ（1988年）」「記憶屋ジョニィ（1981年）」「クローム襲撃（1982年）」、あるいは「記憶屋ジョニィ」を基にした映画「ＪＭ(Johnny Mnemonic,1995年）での用語。それらの作品において、コンピューター・ネットワーク上のサイバースペース（電脳空間）に築かれた「仮想現実空間」、人類の全コンピューター・システムから引き出されたデータの「視覚的再現」、「共感覚幻想」のことを「マトリックス」と呼んでいる。 上掲のサイバーパンク文学上の文脈から映画「マトリックス三部作(1999年~2003年)」でも「コンピュータの作り出した仮想現実」をそう呼んでいる。ちなみに「MATRIX」という言葉自体はボードリヤールの著書「シミュラークルとシミュレーション(Simulacres et simulation,1981年)」の中にも掲げられており、そこに登場する「シミュラークル=不換紙幣(本位貨幣すなわち正貨たる金貨や銀貨との兌換が保障されない法定紙幣（Fiat Money）で、政府の信用において流通する事から信用紙幣（英: Faith Money）とも呼ばれる)の様に、単に特定の実体と紐付けられたラベルに留まらない記号そのもの」の概念こそが出所となったという見方もある。作中ではハードカバーのボードリヤールの本が映るシーンも見られ、2作目からボードリヤール本人をアドバイザーに迎える計画があったが、断られたという。【雑想】「自己の専有と疎外の終わり」?【データベース消費】 穿り返すと案外根深い？ これが数理的にどういうニュアンスになるかというと… import sympy as sp a,b,c,x,y = sp.symbols('a,b,c,x,y') A = sp.Matrix([[a,b,c]]) X = sp.Matrix([[x],[y],[-1]]) AX=A*X expr00 =ax[0] expr01 = sp.solve(expr00,y) expr02 = sp.solve(expr00,x) expr03 = sp.solve(expr00,c) sp.init_printing() display(A) display(X) display(AX) print(sp.latex(A)+sp.latex(X)+"="+sp.latex(AX)) display(expr00) print(sp.latex(expr00)+"=0") display(expr01) print("y="+sp.latex(expr01)) display(expr02) print("x="+sp.latex(expr02)) display(expr03) print(sp.latex(expr03)+"=c") \left[\begin{matrix}a & b & c\end{matrix}\right]\left[\begin{matrix}x\\y\\-1\end{matrix}\right]=\left[\begin{matrix}a x + b y - c\end{matrix}\right] \\ a x + b y - c=0 \\ y=\left[ \frac{- a x + c}{b}\right] \\ x=\left[ \frac{- b y + c}{a}\right] \\ \left[ a x + b y\right]=c 全体像を構築しようとすると行列に「-1」なる成分を追加する必要が生じますが、これは複式簿記における「借方(左側)」と「貸方(右側)」の関係の定義に対応します。そう歴史上における複式簿記の導入とは実は経済分野における会計処理の代数計算化に他ならなかったのです!! 【数理考古学】三次方程式(cubic equation)から虚数(Imaginary Number)へ。 例えば300円のコーヒーと350円の紅茶しかメニューが存在しない喫茶店において、1日の売り上げを知りたいとすれば、その日のコーヒーの注文回数をx杯、紅茶の注文回数をy杯と置いて「300x+350y」の計算が成立します。 import sympy as sp a,b,c,x,y = sp.symbols('a,b,c,x,y') A = sp.Matrix([[a,b,c]]) X = sp.Matrix([[x],[y],[-1]]) AX=A*X expr00 =AX[0] expr01 = sp.solve(expr00,y) expr02 = sp.solve(expr00,x) expr03 = sp.solve(expr00,c) Aa=A.subs([(a,300),(b,350)]) expr00a = expr00.subs([(a,300),(b,350)]) expr01a = expr01[0].subs([(a,300),(b,350)]) expr02a = expr02[0].subs([(a,300),(b,350)]) expr03a = expr03[0].subs([(a,300),(b,350)]) sp.init_printing() display(Aa) display(X) display(Aa*X) print(sp.latex(Aa)+sp.latex(X)+"="+sp.latex(Aa*X)) display(expr00a) print(sp.latex(expr00a)+"=0") display(expr01a) print("y="+sp.latex(expr01a)) display(expr02a) print("x="+sp.latex(expr02a)) display(expr03a) print(sp.latex(expr03a)+"=c") \left[\begin{matrix}300 & 350 & c\end{matrix}\right]\left[\begin{matrix}x\\y\\-1\end{matrix}\right]=\left[\begin{matrix}- c + 300 x + 350 y\end{matrix}\right] \\ - c + 300 x + 350 y=0 \\ y=\frac{c}{350} - \frac{6 x}{7} \\ x=\frac{c}{300} - \frac{7 y}{6} \\ 300 x + 350 y=c Sympyが自動的に約分を施すせいで元データすなわち「特定の実体と紐付けられたラベル」のイメージから離れ始めてますね。それではこの喫茶店の1週間の売り上げをランダムに生成してみましょう。とりあえず売り上げ幅は0杯から50杯とします。 import sympy as sp a,b,x,y = sp.symbols('a,b,x,y') A = sp.Matrix([[a,b]]) X = sp.Matrix([[x],[y]]) AX=A*X expr00 =AX[0] Aa=A.subs([(a,300),(b,350)]) Xr=sp.randMatrix(2,7, min=0, max=50) sp.init_printing() display(Aa) display(X) display(Aa*X) print(sp.latex(Aa)+sp.latex(X)+"="+sp.latex(Aa*X)) display(Xr) display(Aa*Xr) print(sp.latex(Aa)+sp.latex(Xr)+"="+sp.latex(Aa*Xr)) \left[\begin{matrix}300 & 350\end{matrix}\right]\left[\begin{matrix}x\\y\end{matrix}\right]=\left[\begin{matrix}300 x + 350 y\end{matrix}\right] \\ \left[\begin{matrix}300 & 350\end{matrix}\right]\left[\begin{matrix}18 & 4 & 41 & 49 & 4 & 41 & 33\\10 & 5 & 6 & 23 & 9 & 6 & 35\end{matrix}\right]=\left[\begin{matrix}8900 & 2950 & 14400 & 22750 & 4350 & 14400 & 22150\end{matrix}\right] ああ、これはもう立派な統計データですね!! かくしてデータ処理の主体がnumpyを介してpadasに推移する。 Pythonで改行コードLF「\n」を削除する import sympy as sp import numpy as np import pandas as pd a,b,x,y = sp.symbols('a,b,x,y') A = sp.Matrix([[300,350]]) X = sp.Matrix([[18,4,41,49,4,41,33],[10,5,6,23,9,6,35]]) X1=X.row_insert(2, A*X) x=X1.transpose() df=pd.DataFrame(np.matrix(x),columns=['Coffe', 'Tea', 'Total']) sp.init_printing() display(A) display(X) display(A*X) print(sp.latex(Aa)+sp.latex(X)+"="+sp.latex(Aa*X)) display(x) org=df.to_html() print(org.replace('\n', '')) \left[\begin{matrix}18 & 10 & 8900\\4 & 5 & 2950\\41 & 6 & 14400\\49 & 23 & 22750\\4 & 9 & 4350\\41 & 6 & 14400\\33 & 35 & 22150\end{matrix}\right] Coffe Tea Total 0 18 10 8900 1 4 5 2950 2 41 6 14400 3 49 23 22750 4 4 9 4350 5 41 6 14400 6 33 35 22150 だから以下の様な記述統計諸元が計算可能となります。 import sympy as sp import numpy as np import pandas as pd X = sp.Matrix([[18,4,41,49,4,41,33],[10,5,6,23,9,6,35]]) df=pd.DataFrame(np.matrix(X.transpose()),columns=['Coffee', 'Tea']) sp.init_printing() display(X.transpose()) print(sp.latex(X.transpose())) org=df.to_html() print(org.replace('\n', '')) print(df.astype({'Coffee': int, 'Tea': int}).describe()) \left[\begin{matrix}18 & 10\\4 & 5\\41 & 6\\49 & 23\\4 & 9\\41 & 6\\33 & 35\end{matrix}\right] Coffee Tea 0 18 10 1 4 5 2 41 6 3 49 23 4 4 9 5 41 6 6 33 35 Coffee Tea count 7.000000 7.000000 mean 27.142857 13.428571 std 18.488091 11.326328 min 4.000000 5.000000 25% 11.000000 6.000000 50% 33.000000 9.000000 75% 41.000000 16.500000 max 49.000000 35.000000 Coffee Tea count 7 7 unique 5 6 top 4 6 freq 2 2 このうちcountは成分の個数(1週間分のデータなので7個)、uniqueはユニークな（一意な）値の要素の個数に対応します。 【初心者向け】記述統計学と代表値 要約統計量についてまとめてみた クォンタイル(パーセンタイル)/四分位(Quantile) 以下の要素から構成される要約統計量(Summary Statistics)。 最小値(Min) 第1四分位(1st Quartile=25%) 中央値(Median=50%) 第3四分位（3rd Quartile=75%） 最大値(Max) 以下の形で統計的代表値の候補となります。 散布度基準の一つである範囲(Range)は最小値(Min)と最大値（Max)から構成される(最大値-最小値)。 分布の中心にある中央値(Median)は、候補データが複数存在する場合にはその平均を取る。年収リストの様な中央値と平均が異なる分布では、概ね中央値の方が平均値を下回る。平均値(Mean)より外れ値の影響に対してロバスト(頑強)とされる。【初心者向け】平均と標準偏差、中央値と平均偏差 中央値(Median)は平均値(Mean)ほど外れ値の影響を受けない事の簡単な検証。 import sympy as sp import numpy as np import pandas as pd df1=pd.DataFrame([300,300,300,300,400,400,400,2000,9000]) df2=pd.DataFrame([300,300,300,300,400,400,400,2000,12000]) print("mean1=" + str(df1.mean())) print("median1=" + str(df1.median())) print("mean2=" + str(df2.mean())) print("median2=" + str(df2.median())) mean1=0 1488.888889 dtype: float64 median1=0 400.0 dtype: float64 mean2=0 1822.222222 dtype: float64 median2=0 400.0 dtype: float64 いわゆる「箱髭図(Box Plot)」はこのデータから起こされます。 【プログラマーのための統計学】箱ひげ図 %matplotlib inline import matplotlib.pyplot as plt # 数学の点数 math = [74, 65, 40, 62, 85, 67, 82, 71, 60, 99] # 国語の点数 literature = [81, 62, 32, 67, 41, 50, 85, 70, 67, 97] # 点数のタプル points = (math, literature) # 箱ひげ図 fig, ax = plt.subplots() bp = ax.boxplot(points) ax.set_xticklabels(['math', 'literature']) plt.title('Box plot') plt.xlabel('exams') plt.ylabel('point') # Y軸のメモリのrange plt.ylim([0,100]) plt.grid() # 描画 plt.show() それぞれの意味はこうなっています。 箱ひげ図の場合、最大値、最小値が極端であった場合、ヒゲが長くなります。外れ値を考慮する箱ひげ図の場合は、ヒゲの長さは、最大値側、最小値側でそれぞれ、箱の1.5倍以下とし、それを超えるデータは、外れ値とみなします。 pythonのmatplotlibでは、外れ値を自動で検出してくれるようです。以下のコードでは、国語の点数結果に170点、190点を追加してみました。テストは100点満点なので、この2つは外れ値になるはずです。グラフの目盛りは200までに増やしています。これでグラフを作成してみます。 %matplotlib inline import matplotlib.pyplot as plt # 国語の点数 literature = [81, 62, 32, 67, 41, 50, 85, 100, 170, 190] # 点数のタプル points = (literature) # 箱ひげ図 fig, ax = plt.subplots() bp = ax.boxplot(points) ax.set_xticklabels(['literature']) plt.title('Box plot') plt.xlabel('exams') plt.ylabel('point') # Y軸のメモリのrange plt.ylim([0,200]) plt.grid() # 描画 plt.show() グラフの上部の方に、◯が2つできました。この2つは、170点、190点が外れ値としてみなされたものです。pythonのmatplotlibでは、特に外れ値を定義しなくても、このように自動で判別してくれるようなので、非常に便利ですね。 独自に「外れ値」を除外する仕組みが備わっていたりもするのですね。ちなみに今回の「とある架空の喫茶店の7日分の売り上げ(データランダム生成)」をプロットするとこんな感じに。 %matplotlib inline import matplotlib.pyplot as plt # コーヒーの売り上げ coffee = [18,4,41,49,4,41,33] # 紅茶の売り上げ tea = [10,5,6,23,9,6,35] # 点数のタプル points = (coffee, tea) # 箱ひげ図 fig, ax = plt.subplots() bp = ax.boxplot(points) ax.set_xticklabels(['coffee', 'tea']) plt.title('Box plot') plt.xlabel('exams') plt.ylabel('point') # Y軸のメモリのrange plt.ylim([0,50]) plt.grid() # 描画 plt.show() なんとランダム生成データなのに外れ値が検出されてしまいました? 範囲(Range)と最頻値(Mode) 以下の要素から構成される要約統計量(Summary Statistics)。 範囲(Range)…散布度基準の一つだが、あまり普及してない。最大値(Max)-最小値(Min)。度数分布（Frequency Distribution）の表現形態の一つたるヒストグラム（Histogram）では「階級(Class)の最大値-階級(Class)の最小値」と計算される。 最頻値(Mode)…日本工業規格の定義では「離散分布の場合は確率関数が，連続分布の場合は密度関数が最大となる確率変数の値。分布が多峰性の場合は，それぞれの極大値を与える確率変数の値」と定めている。Pandasのdescribe関数ではTopが最頻値、freqがその出現回数を表す。 そう、連続分布を度数分布表(Histogram)で離散的に表す過程はしばしば、「データ分布の(確率密度に対応した)クラス分けの最適化」なる分析過程を必要とするのです。 【初心者向け】度数分布と最頻値 試しに「とある架空の喫茶店の7日分の売り上げ(データランダム生成)」をプロットすると… %matplotlib inline import numpy as np import matplotlib.pyplot as plt # コーヒーの売り上げ coffee = [18,4,41,49,4,41,33] # 紅茶の売り上げ tea = [10,5,6,23,9,6,35] plt.hist([coffee, tea], stacked=False) # 描画 plt.show() まぁ当然の帰結として見るからに有意性の感じられない結果を迎えました。これではクラス分けを最適化する意欲も起こりません。 平均(Mean)と偏差（Deviation）と分散（Dispersion） ここで共通して使用する要約統計量(Summary Statistics)は以下です。 平均(Mean) とりあえず「とある架空の喫茶店の7日分の売り上げ(データランダム生成)」における偏差(Deviation)、すなわち「データの各数値より、その平均を引いた残り」概念から出発します。数理的にいうと「平均値(mean)を原点とする分布への射影結果」といったイメージですね。 import sympy as sp import numpy as np import pandas as pd x=np.matrix([8900,2950,14400,22750,4350,14400,22150]) y=pd.DataFrame([8900,2950,14400,22750,4350,14400,22150]) org=y.to_html() print(org.replace('\n', '')) print("元データ:" + str(x)) print("平均:" + str(np.round(np.mean(x)))) print("偏差:" + str(np.round(x-np.mean(x)))) print("標本分散:" + str(np.round(np.var(x)))) print("不偏分散:" + str(np.round(y.var()))) print("標本偏差:" + str(np.round(np.std(x)))) print("標準偏差:" + str(np.round(y.std()))) print("偏差の絶対値:" + str(np.round(abs(x-np.mean(x))))) print("平均偏差=偏差の絶対値の平均:" + str(np.round(np.sum(abs(x-np.mean(x))))/len(x))) 0 0 8900 1 2950 2 14400 3 22750 4 4350 5 14400 6 22150 元データ:[[ 8900 2950 14400 22750 4350 14400 22150]] 平均:12843.0 偏差:[[-3943. -9893. 1557. 9907. -8493. 1557. 9307.]] 標本分散:53595306.0 不偏分散:0 62527857.0 dtype: float64 標本偏差:7321.0 標準偏差:0 7907.0 dtype: float64 偏差の絶対値:[[3943. 9893. 1557. 9907. 8493. 1557. 9307.]] 平均偏差=偏差の絶対値の平均:44657.0 何やらややこしい計算を重ねた末に「標本分散(Sample Dispersion)/標本偏差(Sample Dispersion)」「不偏分散(Unbiased Dispersion)/標準偏差(Standard Deviation)」「平均偏差(Mean Deviation)」の3組に到達します。どういう事なのでしょうか? 順を追って見ていく事にしましょう。 偏差(Deviation)の合計自体は0 統計学的画像再構成法であるOSEMアルゴリズムの基礎論 まずは(母平均μ（ミュー）を中心位置、母標準偏差σ（シグマ）をばらつき具合とする未知の分布によって構成される)母集合(Parameter=パラメータ)からn個の標本集合(Sample Set)$[x_1,…,x_n]$を抽出し、その「標本平均(Sample Mean)」を求めます。 \hat{μ}=\overline{x} = \frac{1}{n}\sum_{ i = 1 }^{ n } x_i この標本集合の偏差(Deviation)すなわち各データから平均値を引いた結果の合計は必ず0となります。 \sum_{ i = 1 }^{ n } x_i-\overline{x}=(x_1-\overline{x})+…+(x_n-\overline{x})=0 試しに「とある架空の喫茶店の7日分の売り上げ(データランダム生成)」で計算すると… import sympy as sp import numpy as np import pandas as pd from scipy import stats a=[8900,2950,14400,22750,4350,14400,22150] print("リスト表現") print(a) print("平均") print(np.round(np.mean(a))) print("偏差") print(a-np.round(np.mean(a))) print("偏差計") print(sum(a-np.round(np.mean(a)))) リスト表現 [8900, 2950, 14400, 22750, 4350, 14400, 22150] 平均 12843.0 偏差 [-3943. -9893. 1557. 9907. -8493. 1557. 9307.] 偏差計 -1.0 四捨五入による誤差の結果、中心が「ちょうど0」にはなりませんが、イメージ的には大体こんな感じ。群論でいうと「加法単位元と加法逆元が設定された加法群」に射影された状態に該当する様です。統計学における標準化(Standardization)にはそういうニュアンスもあったのですね。 ただし、この考え方ではデータのばらつき具合を把握出来ないので、さらに考え方を発展させなくてはいけません。 標本分散(Sample Dispersion)/標本偏差(Sample Dispersion) numpyがデフォルトで立拠する分散と偏差の概念です。 標本分散(Sample Dispersion)…偏差^2の合計/標本数 s^2 = \frac{1}{n}\sum_{i=1}^{n} (x_i – \bar{x})^2 標本偏差(Sample Dispersion)…標本分散の平方根 s = \sqrt{\frac{1}{n}\displaystyle \sum_{ i = 1 }^{ n } (x_i-\overline{x})^2} 上掲の「架空の喫茶店NO枚日の売り上げ」の例でいうと、その開店が1週間限定だったら全データが列記されている事になり、こちらの計算が採用されます。 不偏分散(Sample Dispersion)/標準偏差(Sample Dispersion) padasが立拠する分散と偏差の概念です。何やら出自が胡散臭い「平均への回帰」理論を使いますが、確かに計算上はこちらの方が正解に近づくのです。 平均への回帰 -Wikipedia 【初心者向け】標本分散と不偏分散 不偏分散(Unbiased Dispersion)…偏差^2の合計/（標本数-1） s^2 = \frac{1}{n-1}\sum_{i=1}^{n} (x_i – \bar{x})^2 標準偏差(Standard Deviation)…不偏分散の平方根 S = \sqrt{\frac{1}{n-1}\displaystyle \sum_{ i = 1 }^{ n } (x_i-\overline{x})^2} 上掲の「架空の喫茶店NO枚日の売り上げ」の例でいうと、そのデータは連日の売り上げ集合から7日分だけ抽出した結果なら「母分散は別に実際する」と想定せざるを得なくなり、こちらの計算が採用されます。確かにまぁ、普通はそう考えるものなんです? 「統計学における自由度」の考え方自体はこちらを参照して下さい。 自由度 | 統計学の時間 | 統計WEB - BellCurve ちなみにいわゆる「正規分布(Normal Distribution)」が平均(中央値0)同様にパラメーターとして使用する分散(中央値1)は、この不偏分散の方です。【初心者向け】正規分布とは何か？ ちなみにnumpyのstd関数も引数としてddof=1を与える事で不偏分散/標準偏差を扱う様になります。ddofはDelta Degrees of Freedomの略との事。numpy.varやnumpy.stdの自由度 import sympy as sp import numpy as np import pandas as pd from scipy import stats a=[8900,2950,14400,22750,4350,14400,22150] print("リスト表現") print(a) print("np var ddof=0") print(np.var(a,ddof=0)) print("np std ddof=0") print(np.std(a,ddof=0)) print("np var ddof=1") print(np.var(a,ddof=1)) print("np std ddof=1") print(np.std(a,ddof=1)) x=pd.DataFrame(a,columns=["Daily Total"]) print("Dataframe表現") print(x) print("np var ddof=0") print(np.round(x.var(),decimals=2)) print("np var ddof=1") print(np.round(x.std(),decimals=2)) リスト表現 [8900, 2950, 14400, 22750, 4350, 14400, 22150] np var ddof=0 53595306.12244898 np std ddof=0 7320.88151266287 np var ddof=1 62527857.14285714 np std ddof=1 7907.455794555993 Dataframe表現 Daily Total 0 8900 1 2950 2 14400 3 22750 4 4350 5 14400 6 22150 np var ddof=0 Daily Total 62527857.14 dtype: float64 np var ddof=1 Daily Total 7907.46 dtype: float64 ところで「原点ゼロからの自乗を取る」とは行列演算の二次形式で表現すると以下となります。 【Python演算処理】行列演算の基本③パスカルの三角形から二次形式へ import sympy as sp a,b = sp.symbols('a,b') #Rows x=sp.Matrix([[a,b]]) #Columns y=sp.Matrix([[a],[b]]) a=sp.Matrix([[0,0],[0,1]]) sp.init_printing() display(x) display(a) display(y) display(sp.simplify(x*a*y)) print(sp.latex(x)+sp.latex(a)+sp.latex(y)+"="+sp.latex(sp.simplify(x*a*y))) \left[\begin{matrix}a & b\end{matrix}\right]\left[\begin{matrix}0 & 0\\0 & 1\end{matrix}\right]\left[\begin{matrix}a\\b\end{matrix}\right]=\left[\begin{matrix}b^{2}\end{matrix}\right] ユークリッド距離演算$\sqrt{0+x^2}$によって等長性が保証された空間への射影とも見て取れる訳ですね。そもそもいわゆる正規分布のBellCurve自体がこの方法によって(極座標系における)円状分布から平面に射影された結果である事に対応する様なのです。 【初心者向け】正規分布(Normal Distribution)とは何か？ ただしこの空間自体は(中心を定める為の)スカラー演算と(そうして検出された中心からの距離を扱う)加減算の概念こそ存在するものの(ベクトル空間ではある)、まだ角度の概念を備えていない（内積が計算可能なベクトル計量空間ではない）。 平均偏差(Mean Deviation) 標準偏差も標本偏差も平均(Mean)の概念を代表値に選ぶ点は同じですが、平均偏差は中央値(Median)を代表値とする要約統計量(Summary Statistics)です。 【初心者向け】平均と標準偏差、中央値と平均偏差 平均偏差(Mean Deviation) s = \frac{1}{n}\sum_{i=1}^{n} |x_i – \bar{x}| 試しに「とある架空の喫茶店の7日分の売り上げ(データランダム生成)」で計算すると… import sympy as sp import numpy as np import pandas as pd from scipy import stats a=[8900,2950,14400,22750,4350,14400,22150] print("リスト表現") print(a) print("平均") print(np.round(np.mean(a))) print("偏差") print(a-np.round(np.mean(a))) print("偏差の絶対値") print(abs(a-np.round(np.mean(a)))) print("偏差の絶対値計") print(sum(abs(a-np.round(np.mean(a))))) リスト表現 [8900, 2950, 14400, 22750, 4350, 14400, 22150] 平均 12843.0 偏差 [-3943. -9893. 1557. 9907. -8493. 1557. 9307.] 偏差の絶対値 [3943. 9893. 1557. 9907. 8493. 1557. 9307.] 偏差の絶対値計 44657.0 そういえば正規化(Feature Scaling)には「最小値0、最大値1への正規化(Min-Max Normalization)」なんて考え方もあるのです。 x'=\frac{x-min(x)}{max(x)-min(x)} import sympy as sp import numpy as np import pandas as pd from scipy import stats a=[8900,2950,14400,22750,4350,14400,22150] print("リスト表現") print(a) print("Min(0)-Max(1) Normalization") print((a-np.round(min(a)))/(max(a)-min(a))) リスト表現 [8900, 2950, 14400, 22750, 4350, 14400, 22150] Min(0)-Max(1) Normalization [0.30050505 0. 0.57828283 1. 0.07070707 0.57828283 0.96969697] 群論的志向様式に寄せるなら「最小値-1、最大値1への正規化(Min-Max Normalization)」なんて考え方も出来るかもしれません。 x'=2\frac{x-min(x)}{max(x)-min(x)}-1 import sympy as sp import numpy as np import pandas as pd from scipy import stats a=[8900,2950,14400,22750,4350,14400,22150] print("リスト表現") print(a) print("Min(-1)-Max(1) Normalization") print(((a-np.round(min(a)))/(max(a)-min(a)))*2-1) リスト表現 [8900, 2950, 14400, 22750, 4350, 14400, 22150] Min(-1)-Max(1) Normalization [-0.3989899 -1. 0.15656566 1. -0.85858586 0.15656566 0.93939394] 本来ならこれらは状況によって使い分けるのが正しいのでしょうが、結果として一時期(特に20世紀後半)不偏分散/標準偏差を絶対視するイデオロギーが広まりました。この辺り本当に注意が必要です。 標準得点（z得点）への変換と偏差値（Deviation Value） Z値とは Z値とは、標準偏差の単位で観測統計量とその仮説母集団パラメータの差を測定するZ検定の統計量です。 観測値をZ値に変換することを標準化と呼びます。母集団の観測値を標準化するには、対象の観測値から母集団平均を引き、その結果を母集団の標準偏差で除算します。この計算結果が、対象の観測値に関連付けられるZ値です。 Z値を使用して、帰無仮説を棄却するかどうかを判断できます。帰無仮説を棄却するかどうかを判断するには、Z値を棄却値と比較します。これは、ほとんどの統計の教科書の標準正規表に示されています。棄却値は、両側検定の場合はZ1-α/2、片側検定の場合はZ1-αです。Z値の絶対値が棄却値より大きい場合、帰無仮説を棄却します。そうでない場合、帰無仮説を棄却できません。 Z得点(Z Value)…$\frac{偏差}{標準偏差}$ 偏差値(Deviation Value)…Z得点*10+50 【python】pandasでデータを標準得点（z得点）に変換 データの正規化（標準化）とは、データを分散1、平均0に変換する操作である。 Pythonで正規化・標準化（リスト、NumPy配列、pandas.DataFrame） ①z得点と偏差値の計算(デフォルト) import sympy as sp import numpy as np import pandas as pd from scipy import stats a=[8900,2950,14400,22750,4350,14400,22150] A1=stats.zscore(a) print("z得点") print(A1) A2=A1*10+50 print("偏差値") print(A2) x=pd.DataFrame(a,columns=["Daily Total"]) X1=x.apply(stats.zscore) print("z得点") print(X1) X2=X1*10+50 print("偏差値") print( X2) z得点 \\ [-0.53857683 -1.35132048 0.21269882 1.35327185 -1.16008668 0.21269882 1.27131451] \\ 偏差値 \\ [44.61423172 36.48679517 52.1269882 63.53271848 38.39913318 52.1269882 62.71314505] z得点 Daily Total 0 -0.538577 1 -1.351320 2 0.212699 3 1.353272 4 -1.160087 5 0.212699 6 1.271315 偏差値 Daily Total 0 44.614232 1 36.486795 2 52.126988 3 63.532718 4 38.399133 5 52.126988 6 62.713145 ②z得点と偏差値の計算(ddof=0)=$\frac{偏差}{標本偏差}$(*10+50) import sympy as sp import numpy as np import pandas as pd from scipy import stats a=[8900,2950,14400,22750,4350,14400,22150] A1=stats.zscore(a,ddof=0) print("z得点") print(A1) A2=A1*10+50 print("偏差値") print(A2) x=pd.DataFrame(a,columns=["Daily Total"]) X1=x.apply(stats.zscore,ddof=0) print("z得点") print(X1) X2=X1*10+50 print("偏差値") print( X2) z得点 \\ [-0.53857683 -1.35132048 0.21269882 1.35327185 -1.16008668 0.21269882 1.27131451] \\ 偏差値 \\ [44.61423172 36.48679517 52.1269882 63.53271848 38.39913318 52.1269882 62.71314505] z得点 Daily Total 0 -0.538577 1 -1.351320 2 0.212699 3 1.353272 4 -1.160087 5 0.212699 6 1.271315 偏差値 Daily Total 0 44.614232 1 36.486795 2 52.126988 3 63.532718 4 38.399133 5 52.126988 6 62.713145 ③z得点と偏差値の計算(ddof=1)=$\frac{偏差}{標準偏差}$(*10+50) import sympy as sp import numpy as np import pandas as pd from scipy import stats a=[8900,2950,14400,22750,4350,14400,22150] A1=stats.zscore(a,ddof=1) print("z得点") print(A1) A2=A1*10+50 print("偏差値") print(A2) x=pd.DataFrame(a,columns=["Daily Total"]) X1=x.apply(stats.zscore,ddof=1) print("z得点") print(X1) X2=X1*10+50 print("偏差値") print( X2) z得点 \\ [-0.49862525 -1.25107966 0.19692084 1.25288628 -1.07403157 0.19692084 1.17700852] \\ 偏差値 \\ [45.01374747 37.48920335 51.96920843 62.52886278 39.25968432 51.96920843 61.77008522] z得点 Daily Total 0 -0.498625 1 -1.251080 2 0.196921 3 1.252886 4 -1.074032 5 0.196921 6 1.177009 偏差値 Daily Total 0 45.013747 1 37.489203 2 51.969208 3 62.528863 4 39.259684 5 51.969208 6 61.770085 母集団からのランダムサンプリングの場合は③が使われます。偏差値が百点満点のテスト(最低得点0点,最高得点100点)のイメージに立脚する評価空間ながら、大元のZ得点の特徴から0点以下や100点以上も扱える点に注意して下さい。しかしまぁ大抵は例えデータとして現れても外れ値(Outliers)として切り捨てられてしまう事でしょう。 言葉としての「偏差値」が20世紀の受験戦争以降一般にも定着しましたが、その過程でその内容についての誤解も広まりました。この辺りもまた注意が必要です。【無限遠点を巡る数理】無限遠点としての正規分布と分散概念の歴史 こうして全体像を俯瞰してみると統計データにおける標準化(スカラー倍率を抽出して増減範囲を1にまとめ行列演算的に対角化する)と行列演算における標準化(同じくスカラー倍率を抽出して値1の直行行列概念を援用)にある種の相似性が見て取れるのです。２次形式の標準化 正規分布(Normal Distribution)を母集合分布として想定した展開 正規分布を母集合分布として想定すると何が嬉しいって、以下の様な計算が可能となるのです。 標準誤差（SE：Standard Srror） 標準偏差と標準誤差 推定量の標準偏差であり、標本から得られる推定量そのもののバラつき（＝精度）を表します。一般に「標本平均の標準偏差」を意味し、その計算に「中心極限定理」を使います。 「中心極限定理」は、平均μ、分散σに従う母集団からサンプルサイズnの標本を抽出する時、その平均値の分布はnが大きくなるにつれ正規分布$N(μ,\frac{σ^2}{N})$に近づくと考えます。すなわち標本平均の標準偏差は以下に近づくのです。 \sqrt{\frac{σ^2}{N}}=\frac{σ}{\sqrt{N}} ただし、標本の分散は$σ^2$ではなく不偏分散$s^2$を用いることから、標本平均の標準偏差（＝標準誤差SE）は標準偏差sを用いて次の式から計算できます。 SE=\frac{s}{\sqrt{N}}=\frac{\sqrt{\frac{1}{n-1}\displaystyle \sum_{ i = 1 }^{ n } (x_i-\overline{x})^2}}{\sqrt{N}} 標準誤差は、母集団から抽出された標本から標本平均を求める場合、標本平均の値が母平均に対してどの程度ばらついているかを表すものです。サンプルサイズが大きくなると標準誤差は小さくなります。 「統計備忘録」 標準誤差 標準誤差はｎの二乗根に反比例することになりますから、サンプルサイズを４倍にすれば標準誤差は半分になります。 import sympy as sp import numpy as np import pandas as pd from scipy import stats a=[126,224,34,25,199,89,178,14,38,11] A1=np.std(a,ddof=1) print("標準偏差") print(A1) A2=A1/np.sqrt(len(a)) print("標準誤差") print(A2) 標準偏差 82.202730422346 標準誤差 25.994785801942832 統計的仮説検定(Testing of Statistical Hypothes) Pythonで統計学を学ぶ(4) 以下の手順に従う。 ①母集団に関する帰無仮説と対立仮説を設定する。 「帰無仮説」とは提案する手法が従来の手法と「差がない」、すなわち提案する手法は「効果がない」という本来主張したいこととは逆の仮説である。 この仮説が棄却されることを目標として仮説検定を行う。具体的には、母平均μ=0(母平均は0である),母相関係数ρ=0(相関がない),母平均の差μ1−μ2=0(差がない)というような仮説を指す。 一方「対立仮説」とは帰無仮説が棄却されたときに採択される逆の仮説であり、実験などで示したい・主張したいことを表したもの。具体的には、母平均μ≠0(母平均は0でない),母相関係数ρ≠0(相関がある),母平均の差μ1−μ2≠0(差がある)というような仮説を指す。 対立仮説の設定により、検定は次のどちらかで行う(両側検定の方がより厳しい条件であり、普通は両側検定で行う)。 片側検定：対立仮説が、母平均μ＞0(もしくはμ＜0) 、母相関係数ρ＞0(もしくはρ＜0)、母平均の差μ1＞μ2(もしくは μ1＜μ2)の場合 両側検定：対立仮説が、母平均μ≠0、母相関係数ρ≠0、母平均の差μ1-μ2≠0の場合。要するに、両側検定では例えば母平均μ≠0を調べるには、母平均μ>0とμ<0の両方を調べなければならない。 帰無仮説が正しいものとして分析を行う。 実際に得られたデータから計算された検定統計量の値によって採択を判断する。帰無仮説が正しいとしたとき、検定統計量が、ほぼ起こり得ない値(それほど極端な値)であれば、帰無仮説を棄却する(つまり、本来の主張を表す対立仮説が採択される)。そうでなければ(確率的に十分起こりうるような値であれば、帰無仮説を採択する(この場合は、本来主張したかった対立仮説が棄却されてしまう)。 ②検定統計量を選択する。 「検定統計量」とは統計的仮説検定のために用いられる標本統計量のこと。代表的な検定統計量の例:t,χ2、F ③有意水準αの値を決定する。 「有意水準」とは対立仮説を採択するか決定するときに基準である(αで表されます)。 有意水準は5％または1％(α=0.05またはα=0.01)に設定することが多い(つまり、標本が100回に5回(5%の場合)以下にしか現れないデータであった---こんなことは偶然では起こりえない---、だから帰無仮説が成り立たないと考えて良いのではないか、という判断基準) 帰無仮説が正しいものとして考えた時の標本分布を「帰無分布」という---帰無分布に基づいて確率が計算される。 棄却域以外の部分すなわち「帰無仮説が採択される領域」を「採択域」という。 棄却域と採択域の境目の値を「臨界値」という。 帰無仮説のもとで、非常に生じにくい検定統計量の値の範囲を「棄却域」という---帰無仮説が棄却される領域(だから、この範囲に入るのが『望ましい』)。 ④データを収集した後、データから検定統計量の実現値を求める。 「検定統計量の実現値」とは実際のデータ（手に入った標本）を基に計算してえられる具体的な値のことで、対立仮説に合うほど 0から離れた値を示す。 ⑤検定統計量の実現値が棄却域に入れば帰無仮説を棄却し、対立仮説を採択する。そうでなければ、帰無仮説を採択する。 帰無仮説が正しいという仮定のもとで、標本から計算した検定統計量の実現値以上の値が得られる確率を「p値」という。p値が有意水準より小さい時に帰無仮説を棄却する。 p値の大きさこそが対立仮説を採択する(帰無仮説を棄却する)決め手となる。p値が小さいということは 『帰無仮説が正しいとすると』確率的にほとんど起こりえないことが起きた(有意水準が5%なら100回中5回以下、1%なら100回中1回以下)ということを意味する。逆にp値が大きいということは、確率的にはよくあることが起きた(だから、この結果では差があるとはいえない)ということになる。 検定統計量の実現値が棄却域に入った場合、「差がない」という帰無仮説を棄却し、 「差がある」という対立仮説を採択する。この時は以下の様に記述する。 [検定結果は5%(または1%)水準で有意である」 「p<.05(またはp<.01)で有意差が見られた」 。 帰無仮説が棄却できない場合は、「検定の結果、差が有意でなかった」または「有意差が認められなかった」と書く。 import sympy as sp import numpy as np import pandas as pd from scipy import stats #「心理学テスト」が N(12,10) の正規分布に従うものとする。 #次のデータ(「指導法データ」と呼ぶ)はこの母集団から #無作為抽出した標本と考えてよいかどうかを判定せよ SampleData = np.array([13,14,7,12,10,6,8,15,4,14,9,6,10,12,5,12,8,8,12,15]) print("判定データ") print(SampleData) #帰無仮説「無作為抽出した標本と考えて良い」つまり μ=12 #対立仮説「無作為抽出した標本ではない」つまり μ≠12 #検定統計量の選択: 標本データを標準化した値(Zで表す) #有意水準の決定: 両側検定で、有意水準 5%、つまり α=0.05 z = (np.mean(SampleData) - 12) / (10.0/len(SampleData))**0.5 # 標準化 print("標準化データZ") print(z) #帰無仮説によればこの標本は正規分布に従う。 print("下側確率") print(stats.norm.ppf(0.025)) # 下側確率0.05/2 = 0.025となるzの値を求める # 下側確率であるから、この値よりもZ値が小さければ棄却される print("上側確率") print(stats.norm.ppf(0.975)) # 上側確率 1 - 0.05/2 = 0.975となるzの値を求める # 上側確率であるから、この値よりもZ値が大きければ棄却される #棄却域は Z<−1.959964 または Z>1.959964 #従ってZ の値は棄却域に入る。 #関数cdfを用いて、直接p値を求めた場合 print("下側確率(p値)") print(stats.norm.cdf(-2.828427)) # 下側確率とすれば、p値は0.0023という小さな値(< 0.05) print("上側確率(p値)") print(stats.norm.cdf(2.828427)) # 両側検定なので2倍する print("両側検定(p値)") print(2*stats.norm.cdf(-2.828427)) # 両側検定であるから2倍したp値は0.0047という小さな値(< 0.05) #よって、結論は以下となる。 #「有意水準5%において、指導法データは心理学テスト #(という母集団)から無作為抽出した標本とはいえない」。 判定データ [13 14 7 12 10 6 8 15 4 14 9 6 10 12 5 12 8 8 12 15] 標準化データZ -2.82842712474619 下側確率 -1.9599639845400545 上側確率 1.959963984540054 下側確率(p値) 0.0023388684020295768 上側確率(p値) 0.9976611315979704 両側検定(p値) 0.0046777368040591535 こういう定義もありますが… 統計科学事典，清水良一訳(https://www.ism.ac.jp/~ayaka/2017_gairon_1.pdf) 記述統計学とはデータのもっている主要な特性をより鮮明に表現するために，データを要約したり作表をしたりすること一般を指する。 まぁこの辺りが「とりあえず手元データ自体を母集合と考える」記述統計学(Descriptive Statistics)と「手元データは母集合から抽出した部分標本に過ぎない」と考える推計統計学(Inferential Statistics)の基本的態度の境界線となってくる訳で、そこに「大数の法則」概念より発展した正規分布の概念が密接に絡んでくる訳です。 記述統計学と推計統計学の違い 相関係数(Correlation Coefficient)から単回帰分析(Simple Regression Analysis)へ。 相関係数の意味と求め方 - 公式と計算例 \begin{align*} r &= \frac{s_{xy}}{s_xs_y} \\[5pt] &= \frac{\frac{1}{n}\sum_{i=1}^n(x_i-\overline{x})(y_i-\overline{y})}{\sqrt{\frac{1}{n}\sum_{i=1}^n(x_i-\overline{x})^2}\sqrt{\frac{1}{n}\sum_{i=1}^n(y_i-\overline{y})^2}} \end{align*}=\frac{\sum_{i=1}^n(x_i-\overline{x})(y_i-\overline{y})}{\sqrt{\sum_{i=1}^n(x_i-\overline{x})^2}\sqrt{\sum_{i=1}^n(y_i-\overline{y})^2}} 相関係数を求めるには、共分散をそれぞれの変数の標準偏差で割ります。 具体的には次の6つのステップで求めます。 ①それぞれの変数の平均値を求める。 ②それぞれの変数の偏差（数値－平均値）を求める。 ③それぞれの変数の分散（偏差の二乗平均）を求める。 ④それぞれの変数の標準偏差（分散の正の平方根）を求める。 ⑤共分散（偏差の積の平均）を求める。 ⑥共分散を2つの変数の標準偏差の積で割って相関係数を得る。 式にもあります様に($\frac{1}{n}$項だろうが$\frac{1}{n-1}$項だろうが最終的に分母と分子で打ち消し合うので)標本相関係数(Sample Correlation Coefficient)と不偏相関係数(Unbiased Correlation Coefficient)は一致します。 import sympy as sp import numpy as np import pandas as pd from scipy import stats #元データ #A~Eさんの英語(x)と数学(y)のテスト結果。 x=[50,60,70,80,90] y=[40,70,90,60,100] print("データx") print(x) print("データy") print(y) #①それぞれの変数の平均値を求める。 x_mean=np.mean(x) print("x平均") print(x_mean) y_mean=np.mean(y) print("y平均") print(y_mean) #②それぞれの変数の偏差（数値－平均値）を求める。 x_deviation=x-x_mean print("x偏差") print(x_deviation) y_deviation=y-y_mean print("y偏差") print(y_deviation) #③それぞれの変数の分散（偏差の二乗平均）を求める。 x_sample_var=np.var(x,ddof=0) x_unbiased_var=np.var(x,ddof=1) print("x標本分散") print(x_sample_var) print("x不偏分散") print(x_unbiased_var) y_sample_var=np.var(y,ddof=0) y_unbiased_var=np.var(y,ddof=1) print("y標本分散") print(y_sample_var) print("y不偏分散") print(y_unbiased_var) #④それぞれの変数の標準偏差（分散の正の平方根）を求める。 x_sample_std=np.std(x,ddof=0) x_unbiased_std=np.std(x,ddof=1) print("x標本偏差") print(x_sample_std) print("x不偏偏差") print(x_unbiased_std) y_sample_std=np.std(y,ddof=0) y_unbiased_std=np.std(y,ddof=1) print("y標本偏差") print(y_sample_std) print("y不偏偏差") print(y_unbiased_std) #⑤**共分散**（偏差の積の平均）を求める。 xy_sample_cov=sum(x_deviation*y_deviation)/len(x) xy_unbiased_cov=sum(x_deviation*y_deviation)/(len(x)-1) print("xy標本共分散") print(xy_sample_cov) print(np.cov(x,y,bias=True)) print("x不偏共分散") print(xy_unbiased_cov) print(np.cov(x,y,bias=False)) #⑥共分散を2つの変数の標準偏差の積で割って相関係数を得る。 xy_sample_cor=sum(x_deviation*y_deviation)/len(x) xy_unbiased_cor=sum(x_deviation*y_deviation)/len(x) print("xy共分散") print(np.corrcoef(x,y)) データx [50, 60, 70, 80, 90] データy [40, 70, 90, 60, 100] x平均 70.0 y平均 72.0 x偏差 [-20. -10. 0. 10. 20.] y偏差 [-32. -2. 18. -12. 28.] x標本分散 200.0 x不偏分散 250.0 y標本分散 456.0 y不偏分散 570.0 x標本偏差 14.142135623730951 x不偏偏差 15.811388300841896 y標本偏差 21.354156504062622 y不偏偏差 23.874672772626646 xy標本共分散 220.0 [[200. 220.] [220. 456.]] x不偏共分散 275.0 [[250. 275.] [275. 570.]] xy共分散 [[1. 0.7284928] [0.7284928 1. ]] numpy.cov()は共分散行列、numpy.corrcoef()は相関行列を返しますが、これはこう見ます。 [[1つ目のデータと1つ目のデータ, 1つ目のデータと2つ目のデータ], [2つ目のデータと1つ目のデータ, 2つ目のデータと2つ目のデータ]] NumPyの共分散を求める関数np.cov関数の使い方 Pythonで相関係数を求める方法を現役エンジニアが解説【初心者向け】 相関係数は−1から1までの値をとり、1であれば2つの変数の値は完全に同期している事になります。対象によって相関係数の意味はかなり変わってきますが、一般的見方は以下。 −1〜−0.7…強い負の相関 −0.7〜−0.4…かなりの負の相関 −0.4〜−0.2…やや負の相関 −0.2〜0.2…ほとんど相関なし 0.2〜0.4…やや正の相関 0.4〜0.7…かなりの正の相関 0.7〜1…強い正の相関 このデータでは「強い正の相関」が検出されたので、さらに最小二乗法(Least Squares Method)による単回帰分析(Simple Regression Analysis)へと進みます。 Numpyだけで回帰分析その4。polyfit()について。 import numpy as np import matplotlib.pyplot as plt from scipy import stats #元データ #A~Eさんの英語(x)と数学(y)のテスト結果。 xo=[50,60,70,80,90] yo=[40,70,90,60,100] print("データx") print(xo) print("データy") print(yo) #Z得点化 x=stats.zscore(x,ddof=1) y=stats.zscore(y,ddof=1) print("データx(標準化)") print(np.round(x,decimals=3)) print("データy(標準化)") print(np.round(y,decimals=3)) # 近似線の係数を算出 k = np.polyfit(x, y, 1) a,b=k #適当にx軸の両端の値 x_reg = np.array([-1.5, 1.5]) #直線の式を作って配列xを代入 y_reg = a * x_reg + b plt.scatter(x,y) plt.plot(x_reg, y_reg, color="red") plt.show() データx [50, 60, 70, 80, 90] データy [40, 70, 90, 60, 100] データx(標準化) [-1.265 -0.632 -0. 0.632 1.265] データy(標準化) [-1.34 -0.084 0.754 -0.503 1.173] 21世紀は数理モデルの時代? 冒頭に述べた様に「行列」の原語たる「Matrix」には「不換紙幣の如く特定の実体の代替物を超えた何か」なる深淵なニュアンスまで存在するものの、かかる考え方を20世紀末に広めたフランスの思想家ボードリヤールなどは、それほどポジティブな意味では捉えていなかった様なのです。この辺り、同時期広まった「仮想(Virtual)」概念に通じるものがあります。当時はまだまだ「人間の知性は数理を超越して実存する」なる信念が健在だったのです。 仮想とは - IT用語辞典 e-Words 実際、当時その「シニフィエ(Signifié=意味されるもの。当時の流行語の一つ)」の領域は随分と被っていました。惜しくも1980年代末から1990年代にかけては著者高齢化に伴うハイファンタジー小説やTV系サイバーパンク文学の衰退、AI冬の時代の到来を迎え、これに有名なウォシャウスキー姉妹(当時は兄弟)監督映画「マトリックス(Matrix)三部作(1999年~2003年)」の成功が続いた訳ですが、かかる時代をDeep Blue(IBMが開発したチェス専用スーパーコンピュータ)が乗り越え(人間の知性の慎ましい模倣ではなく)ただひたすら「数理モデルの洗練」を目指す第三次AIブームが始まる訳です。 シニフィアンとシニフィエ - Wikipedia 「TV系サイバーパンク」という不思議なジャンル AI冬の時代とは!？AIの過去と今後について ディープ・ブルー (コンピュータ) - Wikipedia ディープラーニング - Wikipedia ようやくこうした景色の一端まで辿り着く事ができました。次回以降はさらにこの考え方を広げていきたいと思います。