Pythonメモ

全くのところ、自分用のメモです。

始めの宣言

%matplotlib widget
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()

区切りの良いカラー番号

Pythonとは関係ありませんが、0から255まで等間隔に区切りたいだけです。17の倍数です。255を素因数分解すると、

255 = 17 * 5 * 3

となる訳で。

# 0から255まで4段階に分ける
0, 85, 170, 255,
# 0から255まで6段階に分ける
0, 51, 102, 153, 204, 255,
# 0から255まで16段階に分ける
0, 17, 34, 51, 68, 85, 102, 119, 136, 153, 170, 187, 204, 221, 238, 255

matplotlibのデフォルトカラーコード

'tab:blue', 'tab:orange', 'tab:green', 'tab:red', 'tab:purple', 'tab:brown', 'tab:pink', 'tab:gray', 'tab:olive', 'tab:cyan'

グラフの描画

# seaborn-paper, seaborn-notebook, seaborn-talk, seaborn-posterの順で文字が大きくなる
with plt.style.context('seaborn-poster'):
    fig, axes = plt.subplots(2, 2, figsize=(8,6))
    ax = np.rabel(axes)
    ....
    plt.tight_layout()
    plt.show()

インタラクティブなmatplotlibグラフ表示時のショートカット

Command	Shortcut(s)
Home/Reset	'h', 'r', 'home'
Back	'left', 'c', 'backspace'
Forward	'right', 'v'
Pan/Zoom	'p'
Zoom-to-rect	'o'
Save	's', 'ctrl+s'
Toggle fullscreen	'f', 'ctrl+f'
Toggle major grids	'g'
Toggle minor grids	'G'
Toggle x axis scale (log/linear)	'k', 'L'
Toggle y axis scale (log/linear)	'l'
Close Figure	'ctrl+w', 'cmd+w', 'q'
Constrain pan/zoom to x axis	hold x when panning/zooming with mouse
Constrain pan/zoom to y axis	hold y when panning/zooming with mouse
Preserve aspect ratio	hold CONTROL when panning/zooming with mouse

好みのグラフにする

# グラフのサイズ
plt.rcParams['figure.figsize'] = (12, 5.5)
# 外枠を消す
#plt.rcParams['axes.spines.left'] = False
plt.rcParams['axes.spines.right'] = False
#plt.rcParams['axes.spines.bottom'] = False
plt.rcParams['axes.spines.top'] = False
# 左の目盛りを消す
plt.rcParams['ytick.left'] = False
# 水平グリッドを表示する
plt.rcParams['axes.grid.axis'] = 'y'
plt.rcParams['axes.grid'] = True
plt.rcParams['grid.linestyle'] = 'dotted'
# カラーマップを再定義する
colors = plt.cm.get_cmap('tab20c')(np.arange(0, 1, 0.05))
plt.rcParams['axes.prop_cycle'] = plt.cycler(color=colors)

2020-08-26

Pythonによる分散分析

統計

分散分析（Analysis Of VAriance; ANOVA）は3標本以上の差の検定です。

Pythonの統計パッケージには、scipy.statsやstatsmodelsがあります。最近はpingouinというものもあります。いずれもANOVAをサポートしています。それぞれのパッケージでの使用方法を説明します。

scipy.statsによるANOVA

from scipy import stats

x = [43, 55, 57, 72, 51, 52, 48, 46, 58, 54]
y = [53, 44, 54, 51, 68, 64, 45, 67, 49, 56]
z = [77, 55, 67, 54, 46, 75, 65, 57, 49, 61]
print(stats.f_oneway(x, y, z))

結果

F_onewayResult(statistic=1.6463168290164742, pvalue=0.21152459574247603)

statsmodelsによるANOVA

import pandas as pd
import statsmodels.api as sm
from statsmodels.formula.api import ols

x = [43, 55, 57, 72, 51, 52, 48, 46, 58, 54]
y = [53, 44, 54, 51, 68, 64, 45, 67, 49, 56]
z = [77, 55, 67, 54, 46, 75, 65, 57, 49, 61]

values = x + y + z
groups = ['x'] * len(x) + ['y'] * len(y) + ['z'] * len(z)
data = pd.DataFrame({'values': values, 'groups': groups})

lm = ols('values ~ groups', data=data).fit()
sm.stats.anova_lm(lm, typ=2) # Type 2 Anova DataFrame

結果

               sum_sq    df         F    PR(>F)
groups     271.666667   2.0  1.646317  0.211525
Residual  2227.700000  27.0       NaN       NaN

pingouinによるANOVA

import pandas as pd
import pingouin as pg

x = [43, 55, 57, 72, 51, 52, 48, 46, 58, 54]
y = [53, 44, 54, 51, 68, 64, 45, 67, 49, 56]
z = [77, 55, 67, 54, 46, 75, 65, 57, 49, 61]

values = x + y + z
groups = ['x'] * len(x) + ['y'] * len(y) + ['z'] * len(z)
data = pd.DataFrame({'values': values, 'groups': groups})

print(pg.anova(data, dv='values', between='groups'))

結果

   Source  ddof1  ddof2         F     p-unc       np2
0  groups      2     27  1.646317  0.211525  0.108694

2020-08-15

Pythonによる級内相関係数の計算

統計

検者内・検者間の信頼性を表す指標として級内相関係数（ICC; Intraclass correlation)というのがあります[1, 2]。Pythonの統計パッケージにはscipy.statsやstatsmodelsがありますが、ICCは実装されていない様です。

最近の統計パッケージpingouinにはありますので紹介します。

pingouinはpipでインストールできます。

pip install pingouin

サンプルコードは以下のとおりです。

import pandas as pd
import pingouin as pg

A = [1,1,3,1,1,2,1,2,1,1]
B = [2,1,3,1,3,2,1,3,3,3]
C = [2,3,3,1,1,1,1,2,3,3]
D = [2,3,3,1,1,2,1,2,3,1]
E = [2,3,3,3,3,2,1,2,3,1]
ratings = A + B + C + D + E
raters = ['A'] * len(A) + ['B'] * len(B) + ['C'] * len(C) + \
         ['D'] * len(D) + ['E'] * len(E)
targets = list('abcdefghij') * 5
data = pd.DataFrame({'targets':targets, 'raters':raters, 'ratings':ratings})
icc = pg.intraclass_corr(data=data, targets='targets', 
                         raters='raters', ratings='ratings')
print(icc.set_index('Type'))

出力結果は以下のとおりです。

                   Description       ICC         F  df1  df2      pval  \
Type                                                                     
ICC1    Single raters absolute  0.266854  2.819923    9   40  0.011519   
ICC2      Single random raters  0.280000  3.221007    9   36  0.005798   
ICC3       Single fixed raters  0.307576  3.221007    9   36  0.005798   
ICC1k  Average raters absolute  0.645380  2.819923    9   40  0.011519   
ICC2k    Average random raters  0.660377  3.221007    9   36  0.005798   
ICC3k     Average fixed raters  0.689538  3.221007    9   36  0.005798   

              CI95%  
Type                 
ICC1   [0.03, 0.64]  
ICC2   [0.05, 0.64]  
ICC3   [0.06, 0.67]  
ICC1k   [0.13, 0.9]  
ICC2k   [0.21, 0.9]  
ICC3k  [0.23, 0.91]

[1] Shrout, P. E., & Fleiss, J. L. (1979). Intraclass correlations: uses in assessing rater reliability. Psychological bulletin, 86(2), 420. https://pdfs.semanticscholar.org/b8d4/7b0c0b12dd77543e82e6bf6636ddd335cfea.pdf

[2] 医療系データのとり方・まとめ方、対馬栄輝・石田水里、東京図書 (2013年). https://www.amazon.co.jp/dp/4489021437

2020-04-26

M5StickCでパルスオキシメーターを作る

プログラミング

M5StickCでパルスオキシメーターを作りました。センサにはM5Stack用心拍セEンサユニット（MAX30100）を使用しました。新型コロナウイルス感染症(COVID-19)では自覚症状がないうちに肺炎が悪化するとの報告もありますので、日頃から血中酸素飽和度(SpO2)を測定するのもいいかもしれません。

使用したパーツは以下のとおりです。

M5StickC - スイッチサイエンス
https://www.switch-science.com/catalog/5517/
M5Stack用心拍センサユニット - スイッチサイエンス
https://www.switch-science.com/catalog/5695/

合計金額3,000円以下です。

M5StickCとM5Stack用心拍センサユニットの組み合わせ。Grove用４ピンケーブルで接続する。ケーブルは心拍センサユニットに付属する。

ArduinoのスケッチはMAX30100libライブラリに付属するサンプルスケッチMAX30100_Minimal.icoを改変しました。

MAX30100lib - Arduino Libraries
https://www.arduinolibraries.info/libraries/max30100lib

変更した点は、以下のとおりです。

M5StickCのGroveポートをI2Cとして用いる事を宣言するコマンドを追加しました。これをしないと通信ができませんでした。

Wire.begin(32, 33);

LEDの電流量を50mAから7.6mAに変更しました。これをしないと測定値が飽和してしまい正しいデータが得られません。下記のサイトを参考にしました。

pox.setIRLedCurrent(MAX30100_LED_CURR_7_6MA);

迷走の果て・Tiny Objects 脈拍センサー MAX30100を試す(3)正常動作？？
https://edycube.blog.fc2.com/blog-entry-982.html

M5StickCのLCD表示に関わる部分を追加しました。また、表示の周期を２秒にしました。

M5.Lcd.fillScreen(BLACK);

M5.Lcd.setTextSize(2);
M5.Lcd.setCursor(10, 10);
M5.Lcd.print("SpO2%");
M5.Lcd.setCursor(10, 100);
M5.Lcd.print("PRbpm");
M5.Lcd.setTextSize(5);
M5.Lcd.setCursor(15, 40);
M5.Lcd.printf("%2d", spo2);
M5.Lcd.setTextSize(3);
M5.Lcd.setCursor(20, 130);
M5.Lcd.printf("%3.0f", pulse_rate);