python でパワースペクトルを計算する方法

matplotlib と spicy で FFT を計算する方法

matplotlib と spipy のどちらもパワースペクトルを計算する方法があるが、両者でデフォルトで持っている関数と定義が違う。大概の使い分けは次の通り。

• 単純にパワースペクトルを計算する場合、mlab.psd が一番簡単。
• scipy.fftpack を使うと実部と虚部が得られる。位相の計算など細かいことが可能だか、規格化は自分で行う必要あり。

どちらでやっても、パーセバルの定理から全パワーは同じになるが、両者の規格化の違いを考えないとパワースペクトルの縦軸が完璧に同じにならない。

ここでは単純な三角関数を両方の方法でフーリエ変換して、完璧にmatplotlibとscipyのFFTで結果が同じになる例を示す。おまけで、パワースペクトルから三角波の振幅を求めて、それがもとの三角関数の振幅と同じになること、フィルターを入れるとどういう影響があるかを紹介する。

問題設定

フーリエ変換する前の関数

縦軸は特に意味はないが、電圧V(t)という信号で単一の三角関数を考える。最も簡単な表式として$$y=V(t)=\sin(t)$$を用いる。周期 $f$ の三角関数は、$$\sin(2\pi ft)$$と書けるので、周期 $$f=1/2\pi$$ という設定になる。振幅は1.0である。

これを mlab.psdとscipy.fftpackの両方でフーリエ変換してパワーを計算する。それぞれ、hanning filter をかけた場合も計算し、違いを見てみる。

時間分解能 dt は0.1秒とする。ナイキスト周波数 f_n は、その逆数の半分で、5Hzとなる。

サンプルコード

qiita_fft_matplotlib_scipy_comp.py

#!/usr/bin/env python

"""
#2013-03-12 ; ver 1.0; First version, using Sawada-kun's code
#2013-08-28 ; ver 1.1; refactoring
#2020-03-31 ; ver 1.2; updated for python3
"""
__version__=  '1.2'

import numpy as np
import scipy as sp
import scipy.fftpack as sf
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import optparse

def scipy_fft(inputarray,dt,filtname=None, pltfrag=False):

    bin = len(inputarray)
    if filtname == None:
        filt = np.ones(len(inputarray))
    elif filtname == 'hanning':
#        filt = sp.hanning(len(inputarray))
        filt = np.hanning(len(inputarray))                
    else:
        print('No support for filtname=%s' % filtname)
        exit()

    freq = sf.fftfreq(bin,dt)[0:int(bin/2)]
    df = freq[1]-freq[0]
    fft = sf.fft(inputarray*filt,bin)[0:int(bin/2)]

    real = fft.real
    imag = fft.imag
    psd = np.abs(fft)/np.sqrt(df)/(bin/np.sqrt(2.))

    if pltfrag: 
        binarray = range(0,bin,1)
        F = plt.figure(figsize=(10,8.))
        plt.subplots_adjust(wspace = 0.1, hspace = 0.3, top=0.9, bottom = 0.08, right=0.92, left = 0.1)
        ax = plt.subplot(1,1,1)
        plt.title("check fft of scipy")
        plt.xscale('linear')
        plt.grid(True)
        plt.xlabel(r'number of bins')
        plt.ylabel('FFT input')
        plt.errorbar(binarray, inputarray, fmt='r', label="raw data")
        plt.errorbar(binarray, filt, fmt='b--', label="filter")
        plt.errorbar(binarray, inputarray * filt, fmt='g', label="filtered raw data")
        plt.legend(numpoints=1, frameon=False, loc=1)
        plt.savefig("scipy_rawdata_beforefft.png")
        plt.show()

    return(freq,real,imag,psd)


usage = u'%prog [-t 100] [-d TRUE]'
version = __version__
parser = optparse.OptionParser(usage=usage, version=version)
parser.add_option('-f', '--outputfilename', action='store',type='string',help='output name',metavar='OUTPUTFILENAME', default='normcheck')
parser.add_option('-d', '--debug', action='store_true', help='The flag to show detailed information', metavar='DEBUG', default=False)
parser.add_option('-t', '--timelength', action='store', type='float', help='Time length of input data', metavar='TIMELENGTH', default=20.)
options, args = parser.parse_args()
argc = len(args)

outputfilename =  options.outputfilename
debug =  options.debug
timelength =  options.timelength

timebin = 0.1 # timing resolution
# Define Sine curve
inputf = 1/ (2 * np.pi) # f = 1/2pi -> sin (t)
inputt = 1 / inputf # T = 1/f 
t = np.arange(0.0, timelength*np.pi, timebin)
c = np.sin(t)
fNum = len(c)
Nyquist = 0.5 / timebin
Fbin = Nyquist / (0.5 * fNum)
fftnorm =  1. / (np.sqrt(2) * np.sqrt(Fbin))

print("................. Time domain ................................")   
print("timebin = " + str(timebin) + " (sec) ")
print("frequency = " + str(inputf) + " (Hz) ")
print("period = " + str(inputt) + " (s) ")
print("Nyquist (Hz) = " + str(Nyquist) + " (Hz) ")
print("Frequency bin size (Hz) = " + str(Fbin) + " (Hz) ")
print("FFT norm  = " + str(fftnorm) + "\n")

# Do FFT.
# Hanning Filtered
psd2, freqlist = mlab.psd(c,
                         fNum,
                         1./timebin,
                         window=mlab.window_hanning,
                         sides='onesided',
                         scale_by_freq=True
                         )
psd = np.sqrt(psd2)    

spfreq,spreal,spimag,sppower = scipy_fft(c,timebin,'hanning',True)

# No Filter
psd2_nofil, freqlist_nofil = mlab.psd(c,
                                     fNum,
                                     1./timebin,
                                     window=mlab.window_none,
                                     sides='onesided',
                                     scale_by_freq=True
                                     )
psd_nofil = np.sqrt(psd2_nofil)    
spfreq_nf,spreal_nf,spimag_nf,sppower_nf = scipy_fft(c,timebin,None)

# Get input norm from FFT intergral
print("................. FFT results ................................")       
amp = np.sqrt(2) * np.sqrt(np.sum(psd * psd) * Fbin)
print("Amp  = " + str(amp))
peakval = amp / (np.sqrt(2) * np.sqrt(Fbin) )
print("Peakval  = " + str(peakval))

# (1) Plot Time vs. Raw value
F = plt.figure(figsize=(10,8.))
plt.subplots_adjust(wspace = 0.1, hspace = 0.3, top=0.9, bottom = 0.08, right=0.92, left = 0.1)

ax = plt.subplot(3,1,1)
plt.title("FFT Norm Check")
plt.xscale('linear')
plt.xlabel(r'Time (s)')
plt.ylabel('FFT input (V)')
plt.errorbar(t, c, fmt='r', label="raw data")
plt.errorbar(t, amp * np.ones(len(c)), fmt='b--', label="Amplitude from FFT", alpha=0.5, lw=1)
plt.ylim(-2,2)
plt.legend(numpoints=1, frameon=True, loc="best", fontsize=8)

# (2a) Plot Freq vs. Power (log-log)
ax = plt.subplot(3,1,2)
plt.xscale('log')
plt.yscale('log')
plt.xlabel(r'Frequency (Hz)')
plt.ylabel(r'FFT Output (V/$\sqrt{Hz}$)')
plt.errorbar(freqlist_nofil , psd_nofil, fmt='r-', label="Window None, matplotlib", alpha=0.8)    
plt.errorbar(freqlist, psd, fmt='b-', label="Window Hanning, matplotlib", alpha=0.8)
plt.errorbar(spfreq_nf, sppower_nf, fmt='mo', label="Window None, scipy", ms = 3, alpha=0.8)    
plt.errorbar(spfreq, sppower, fmt='co', label="Window Hanning, scipy", ms = 3, alpha=0.8)
plt.errorbar(freqlist, fftnorm * np.ones(len(freqlist)), fmt='r--', label=r"1./($\sqrt{2} \sqrt{Fbin}$)", alpha=0.5, lw=1)
plt.legend(numpoints=1, frameon=True, loc="best", fontsize=8)

# (2b) Plot Freq vs. Power (lin-lin)
ax = plt.subplot(3,1,3)
plt.xscale('linear')
plt.xlim(0.1,0.2)
plt.yscale('linear')
plt.xlabel(r'Frequency (Hz)')
plt.ylabel(r'FFT Output (V/$\sqrt{Hz}$)')    
plt.errorbar(freqlist_nofil , psd_nofil, fmt='r-', label="Window None, matplotlib", alpha=0.8)    
plt.errorbar(freqlist, psd, fmt='b-', label="Window Hanning, matplotlib", alpha=0.8)
plt.errorbar(spfreq_nf, sppower_nf, fmt='mo', label="Window None, scipy", ms = 3, alpha=0.8)    
plt.errorbar(spfreq, sppower, fmt='co', label="Window Hanning, scipy", ms = 3, alpha=0.8)
plt.errorbar(freqlist, fftnorm * np.ones(len(freqlist)), fmt='r--', label=r"1./($\sqrt{2} \sqrt{Fbin}$)", alpha=0.5, lw=1)
plt.legend(numpoints=1, frameon=True, loc="best", fontsize=8)

plt.savefig("fft_norm_check.png")
plt.show()

実行方法は、単に何も指定しないか、-t で時間幅を変えるくらい。

$ python qiita_fft_matplotlib_scipy_comp.py -h
Usage: qiita_fft_matplotlib_scipy_comp.py [-t 100] [-d TRUE]

Options:
  --version             show program's version number and exit
  -h, --help            show this help message and exit
  -f OUTPUTFILENAME, --outputfilename=OUTPUTFILENAME
                        output name
  -d, --debug           The flag to show detailed information
  -t TIMELENGTH, --timelength=TIMELENGTH
                        Time length of input data

フーリエ変換する前の時系列データ

何もオプションを指定しないで実行すると、次のような時系列データのフーリエ変換を実行する。

hanningフィルタがかかる場合、元の時系列データ(赤)に対して、hanning フィルタ(青)がかかって、緑のデータをフーリエ変換することになる。フィルターを使わない場合は、赤を周期境界条件でフーリエ変換することになる。フィルターをかける場合は、オフセット(周波数=0の成分)が乗っている場合は、フィルターによる人為的な低周波ノイズがのる場合があるので、フィルターを使う場合はDC成分は引いてから使うほうがよい。

パワースペクトルの計算結果

(上段) フーリエ変換する時系列データ。パワースペクトルから計算した時系列データの振幅を重ねて表示した。(中段)赤がフィルターなしの場合、青がhanning filter ありの場合。実線がmatplotlibで計算した場合で、ドットはscipyで計算したもの。赤い薄い線は時系列データから推定したパワーの値。(下段)中段の縦軸リニアバージョン。

パーセバルの定理の確認

パーセバルの定理より、周波数空間と時空間の積分値のパワーは一致する。三角関数の場合は、変動のパワー(RMS)と振幅は$\sqrt{2}$で結ばれているので、振幅とRMSは相互に変換できる。ここでは、周波数空間を全積分したパワーから時系列データの振幅(この例では1)を計算し、それが一致することを上の絵の青点線で示した。

周波数空間で積分している箇所は、コードの、

amp = np.sqrt(2) * np.sqrt(np.sum(psd * psd) * Fbin)
print("Amp  = " + str(amp))

に該当する。$$ amplitude = \sqrt{2} \sqrt{ \Sigma_i psd(f_i)^2 \Delta f }$$ で計算している。この値がこの例では1になっていて、

plt.errorbar(t, amp * np.ones(len(c)), fmt='b--', label="Amplitude from FFT", alpha=0.5, lw=1)

で上段に青い点線でプロットしている。

hanning フィルターの効果

フィルタがなければ、mlab.psdとscipy.fftpack は完璧に一致する。hanning filterを使うとパワーはある割合で下がるのでその分を補正する必要があるのだが、matplotlibとscipyで定義が違うのか(深追いしてない)、デフォルトのまま使うと両者の結果は一致はしない。

コードの補足

mlab.psdの使い方

# Do FFT.
# Hanning Filtered
psd2, freqlist = mlab.psd(c,
                         fNum,
                         1./timebin,
                         window=mlab.window_hanning,
                         sides='onesided',
                         scale_by_freq=True
                         )
psd = np.sqrt(psd2)    

ここで、c はFFTしたい同時間感覚幅で取得されたデータ。FFTはサンプルが全く同じ時間間隔でされたことを前提としている。(もし、時間間隔が変わる場合にはダイレクトに離散フーリエ変換をするが、計算時間はかかる。精密な実験などではそうする場合もある。) fNumはサンプル数、1/timebin で時間幅の逆数を指定する、window は None とすれば何もなし、mlabのオプションにいろんなフィルターがあるが、ここでは hanningの例を示した。sidesは onesided で片面(周波数が正の値)で積分したときにトータルのパワーが時系列のパワーと一致するように規格化するという意味。scale_by_freq で、周波数あたりに直す。これにより、返り値はpower^2/Hzの単位となる。ここでは最後に、psd = np.sqrt(psd2) をすることで、power/sqrt(Hz)の単位にしている。

scipy.fftpack の使い方

    freq = sf.fftfreq(bin,dt)[0:int(bin/2)]
    df = freq[1]-freq[0]
    fft = sf.fft(inputarray*filt,bin)[0:int(bin/2)]

fftfreq という関数を使うことで、周波数を計算してくれる。これは実はとっても大事な機能である。なぜかというと、ここscipyのfftでは複素フーリエ変換を実行するのであるが、返り値は実部と虚部の係数であり、それをどういう周波数の順番でユーザーに返すのかは自明ではないので、対応のとれた周波数を返してくれる関数が必要になる。

    real = fft.real
    imag = fft.imag
    psd = np.abs(fft)/np.sqrt(df)/(bin/np.sqrt(2.))

で、実部を虚部を取り出し、片側(周波数が正)で規格化して帳尻があるように規格化を行う。

ずっと同じカラースキームでターミナルの画面を見ていると飽きるし、日中はライトテーマ・夕方はダークテーマにして気分をリフレッシュしたいと思いませんか？ ですが、毎日いちいち手動で朝晩にテーマを切り替えるのは面倒なのでPythonにやらせます。

Python API

2019年7月にリリースされたiTerm2バージョン3.3.0からPython APIが提供されるようになりました（それ以前はAppleScriptというAPIがあったらしい？）。pip install iterm2 すればiTerm2外部からもAPIを叩けます。

• iTerm2 Python API 公式ドキュメント

APIを使用するには Preferences > General > Magic にある Enable Python API にチェックを入れます。

自前のスクリプトを用意したい場合、Scripts > Manage > New Python Script を選ぶとガイドが出てきて雛形を作ってくれます。
作成したPythonのコードは ~/Library/Application Support/iTerm2/Scripts に保存するようになっていますが、アプリ立ち上げ時に自動で読み込ませたい場合は上記のScripts傘下にAutoLaunchフォルダを作り、そこに保管します。

実はカラースキームの自動切り替えスクリプトはフルスクラッチする必要がなく、公式の Example Scripts にほぼ完成形が用意してあります。

• iTerm2: Change Color Preset by Time of Day

元のサンプルコードだとカラースキームの切り替え時にiTerm2が起動している必要があったので、アプリ起動時にも現在の時刻を参考に切り替えてくれるよう変更してあります。以下のコードを~/Library/Application Support/iTerm2/Scripts/AutoLaunch/change_color.pyとして保存してください。

• Gist: https://gist.github.com/elpnt/b1976b090a0e629bc838eb28fa7fec71

change_color.py

#!/usr/bin/env python3.7

import asyncio
import datetime
import iterm2

# Clock time to change colors.
LIGHT_TIME = (7, 0)
DARK_TIME = (17, 0)

# Color presets to use
LIGHT_PRESET_NAME = "material"
DARK_PRESET_NAME = "onedark"

# Profiles to update
PROFILES = ["Default"]

def get_datetime(t, time):
    return datetime.datetime(t.year, t.month, t.day, time[0], time[1])

def datetime_after(t, time):
    today = get_datetime(t, time)
    if today > t:
        return today
    # Same time tomorrow
    return today + datetime.timedelta(1)


def next_deadline_after(t):
    light_deadline = datetime_after(t, LIGHT_TIME)
    dark_deadline = datetime_after(t, DARK_TIME)
    if light_deadline < dark_deadline:
        return (LIGHT_PRESET_NAME, light_deadline)
    return (DARK_PRESET_NAME, dark_deadline)


def get_duration():
    now = datetime.datetime.now()
    preset_name, deadline = next_deadline_after(now)
    duration = (deadline - now).seconds
    print("Sleep for {} seconds until {}".format(duration, deadline))
    return duration, preset_name


async def set_colors(connection, preset_name):
    print("Change to preset {}".format(preset_name))
    preset = await iterm2.ColorPreset.async_get(connection, preset_name)
    for partial in (await iterm2.PartialProfile.async_query(connection)):
        if partial.name in PROFILES:
            await partial.async_set_color_preset(preset)


async def main(connection):
    now = datetime.datetime.now()
    begin = get_datetime(now, LIGHT_TIME)
    end = get_datetime(now, DARK_TIME)
    if (now > begin and now < end):
        await set_colors(connection, LIGHT_PRESET_NAME)
    else:
        await set_colors(connection, DARK_PRESET_NAME)
    while True:
        duration, preset_name = get_duration()
        await asyncio.sleep(duration)
        await set_colors(connection, preset_name)
        await asyncio.sleep(1)


iterm2.run_forever(main)

最初に定義した5つの定数

• LIGHT_TIME=(7, 0) : ライトテーマに切り替える時間(HH, MM)
• DARK_TIME=(16, 0) : ダークテーマに切り替える時間(HH, MM)
• LIGHT_PRESET_NAME="..." : ライトテーマ名（デフォルトテーマはLight Background）
• DARK_PRESET_NAME="..." : ダークテーマ名（デフォルトテーマはDark Background）
• PROFILES=["..."] : テーマを変更させるプロファイル名(複数選択可)

は好みに合わせて変更してください。プロファイル名は自分の作ってあるプロファイル一覧から、テーマのプリセット名は右下の Color Presets... メニューから選びます。

ちなみに print(...) の標準出力や、存在しないプリセット名を指定してしまったときなどのエラーは Scripts > Manage > Console に吐き出されます。

退屈なことPythonにやらせるのサイコ〜〜〜〜〜〜

参考

iTerm2 にステータスバーが付いた

記事の全体像

この記事は2部構成になっています。

・第一部：EDA・前処理 ← いまここ
・第二部：ベクトル化・モデル作成編(LSTM、BERT) (工事中、Coming Soon!)

第二部では第一部で前処理を行ったデータを用いる為、
記事の最下部に前処理を一括で実行できるコードを記載しているので、
第一部をスキップされたい方はそちらを実行して第二部に進んでください。
(※データはこの記事を参照して入手してください！)

コンペティションの概要

Kaggle : Real or Not? NLP with Disaster Tweets
自然言語処理の入門コンペです。
タスクは、ツイートを「災害ツイート」or「災害以外のツイート」の2つに分類することです。

最近では、災害時の救助要請などにツイッターが使われており、災害救援団体や通信社などで災害発生時のツイート自動監視の関心が高まっています。
しかし、ツイートが実際に災害を表しているかどうかを機械的に判断することは難しいです。
それは、例えば災害を明示的に表す「燃えている」という表現も、比喩表現として「空が燃えている」のように使われることがあるからです。
このコンペでは、10,000のツイートのデータセットを使用して、「災害ツイート」or「災害以外のツイート」を予測する機械学習モデルを構築します。

いざハンズオン！

Jupyter Notebookを使用してハンズオンを進めていきます。

0. データの入手

Kaggle : Real or Not? NLP with Disaster Tweets - Data
上述のページのから
・test.csv
・train.csv
の2つをダウンロードします。
以後は、分析を行うノートブックと同じ階層にこのcsvファイルを置いてください。

csvファイルの各カラムの内容は以下の通り。

ファイル名	説明
id	各ツイートの一意の識別子
text	ツイートのテキスト本文
location	ツイートが送信された場所(空欄有り)
keyword	ツイート内の特定のキーワード(空欄有り)
target	ツイートのラベル(災害ツイート=1、災害以外のツイート=0)

具体的な内容や欠損値の情報は次項で説明します。
コンペでは、test.csvのtargetカラムを予測します。

1.ライブラリの読み込み

今回のハンズオンで使用するライブラリを事前に定義します。

# データ分析
import pandas as pd
import numpy as np

# 可視化
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns

# 自然言語処理
import string # 記号(punctuation)の一覧を取得
import re
import contractions # 文書の短縮形を
from wordcloud import STOPWORDS # ストップワードのリストを取得する
from collections import defaultdict # n-gram作成時に利用

2.データの内容確認

早速データの内容を確認してみたいと思います。
まずは、データをデータフレームとして読み込みます。
データフレームの行数・列数を表示して、行をランダムで10個抽出してみます。

# 訓練データとテストデータの読み込み
df_train = pd.read_csv('train.csv', dtype={'id': np.int16, 'target': np.int8})
df_test = pd.read_csv('test.csv', dtype={'id': np.int16})

# 訓練データとテストデータの行数・列数を表示する
print('Training Set Shape = {}'.format(df_train.shape))
print('Test Set Shape = {}'.format(df_test.shape))

# 訓練データの中から10行をランダムで抽出
df_train.sample(n=10, random_state=28)

実行結果の一部

Training Set Shape = (7613, 5)
Test Set Shape = (3263, 4)

以下、実行結果の続き(ランダムで抽出された行の内容)

keywordとlocationのカラムの中には欠損値NaN がある場合もありますね。

欠損値について詳細に調べてみたいと思います。

# 訓練データ、テストデータの各カラムの欠損値率を算出
print("missing-value ratio of training data(%)")
print(df_train.isnull().sum()/df_train.shape[0]*100)
print("\nmissing-value ratio of test data(%)")
print(df_test.isnull().sum()/df_test.shape[0]*100)

実行結果

missing-value ratio of training data(%)
id           0.000000
keyword      0.801261
location    33.272035
text         0.000000
target       0.000000
dtype: float64

missing-value ratio of test data(%)
id           0.000000
keyword      0.796813
location    33.864542
text         0.000000
dtype: float64

欠損値は、訓練データとテストデータいずれにおいても keyword で0.8%、location で33~34%あることが分かります。

次に、訓練データのtargetの分布を見てみましょう。

# ターゲットの要素とその個数をプロット
target_vals = df_train.target.value_counts()
sns.barplot(target_vals.index, target_vals)
plt.gca().set_ylabel('samples')

データセット中には災害ツイート=1よりも災害以外のツイート=0の方が多いことが分かります。

次に、text、keyword、locationそれぞれのカラムのユニークな要素の個数を調べてみたいと思います。

# text、keyword、locationのユニークな要素の個数を表示
print(f'Number of unique values in text = {df_train["text"].nunique()} (Training) - {df_test["text"].nunique()} (Test)')
print(f'Number of unique values in keyword = {df_train["keyword"].nunique()} (Training) - {df_test["keyword"].nunique()} (Test)')
print(f'Number of unique values in location = {df_train["location"].nunique()} (Training) - {df_test["location"].nunique()} (Test)')

実行結果

Number of unique values in text = 7503 (Training) - 3243 (Test)
Number of unique values in keyword = 221 (Training) - 221 (Test)
Number of unique values in location = 3342 (Training) - 1603 (Test)

textとlocationは自由入力であることが分かります。
一方でkeywordはtextの中から事前に定義された221のキーワードを自動抽出していることが分かります。

3.探索的データ分析(Explanatory Data Analysis：EDA)

データの特徴を理解するための処理をいくつか行っていきたいと思います。
まずはtext内の特徴を大雑把に把握・比較したいと思います。
特徴は、単語数、ユニークな単語数、ストップワードの数、URLの数、単語文字数の平均、文字数、punctuation(※)の数、ハッシュタグ数、メンション数の9個を調べてみます。
1点補足をしておくと、ここでいうpunctuation(※)とは句読点の意味ではなく、英数字以外のアスキー文字のことを指します。
ざっくりと、string.punctuationで定義される「記号」のようなイメージを持ってもらえれば大丈夫です。
また、これらの特徴はメタ特徴量として、データフレームに結合していきます。
更に、訓練データ内で災害ツイート=1⇄災害以外のツイート=0、訓練データ⇄テストデータ
について9個の特徴の分布を比較します。
災害に関連するツイート、そうでいないツイート、訓練データ、テストデータはそれぞれデータ点数が異なるため、
分布の可視化を行う際には、カーネル密度推定を行うことによってスケールが揃えられて直感的に分布を比較することができます。
可視化のライブラリであるseabornのdistplotのデフォルト引数でkde=Trueとなっていますが、今回は明示的にkde=Trueを記載します。

# 単語数
df_train['word_count'] = df_train['text'].apply(lambda x: len(str(x).split()))
df_test['word_count'] = df_test['text'].apply(lambda x: len(str(x).split()))

# ユニークな単語数
df_train['unique_word_count'] = df_train['text'].apply(lambda x: len(set(str(x).split())))
df_test['unique_word_count'] = df_test['text'].apply(lambda x: len(set(str(x).split())))

# ストップワードの数
df_train['stop_word_count'] = df_train['text'].apply(lambda x: len([w for w in str(x).lower().split() if w in STOPWORDS]))
df_test['stop_word_count'] = df_test['text'].apply(lambda x: len([w for w in str(x).lower().split() if w in STOPWORDS]))

# URLの数
df_train['url_count'] = df_train['text'].apply(lambda x: len([w for w in str(x).lower().split() if 'http' in w or 'https' in w]))
df_test['url_count'] = df_test['text'].apply(lambda x: len([w for w in str(x).lower().split() if 'http' in w or 'https' in w]))

# 単語文字数の平均
df_train['mean_word_length'] = df_train['text'].apply(lambda x: np.mean([len(w) for w in str(x).split()]))
df_test['mean_word_length'] = df_test['text'].apply(lambda x: np.mean([len(w) for w in str(x).split()]))

# 文字数
df_train['char_count'] = df_train['text'].apply(lambda x: len(str(x)))
df_test['char_count'] = df_test['text'].apply(lambda x: len(str(x)))

# 句読点の個数
df_train['punctuation_count'] = df_train['text'].apply(lambda x: len([c for c in str(x) if c in string.punctuation]))
df_test['punctuation_count'] = df_test['text'].apply(lambda x: len([c for c in str(x) if c in string.punctuation]))

# ハッシュタグの個数
df_train['hashtag_count'] = df_train['text'].apply(lambda x: len([c for c in str(x) if c == '#']))
df_test['hashtag_count'] = df_test['text'].apply(lambda x: len([c for c in str(x) if c == '#']))

# メンションの個数
df_train['mention_count'] = df_train['text'].apply(lambda x: len([c for c in str(x) if c == '@']))
df_test['mention_count'] = df_test['text'].apply(lambda x: len([c for c in str(x) if c == '@']))

# 9つの特徴の分布を、災害ツイート=1 ⇄ 災害以外のツイート=0、訓練データ ⇄ テストデータで比較する
METAFEATURES = ['word_count', 'unique_word_count', 'stop_word_count', 'url_count', 'mean_word_length',
                'char_count', 'punctuation_count', 'hashtag_count', 'mention_count']
DISASTER_TWEETS = df_train['target'] == 1

fig, axes = plt.subplots(ncols=2, nrows=len(METAFEATURES), figsize=(20, 50), dpi=100)

for i, feature in enumerate(METAFEATURES):
    # 災害ツイート=1 ⇄ 災害以外のツイート=0の分布を比較する(カーネル密度推定を行う)
    sns.distplot(df_train.loc[~DISASTER_TWEETS][feature], label='Not Disaster', ax=axes[i][0], color='green', kde=True)
    sns.distplot(df_train.loc[DISASTER_TWEETS][feature], label='Disaster', ax=axes[i][0], color='red', kde=True)

    # 訓練データ ⇄ テストデータの分布を比較する(カーネル密度推定を行う)
    sns.distplot(df_train[feature], label='Training', ax=axes[i][1], kde=True)
    sns.distplot(df_test[feature], label='Test', ax=axes[i][1], kde=True)

    for j in range(2):
        axes[i][j].set_xlabel('')
        axes[i][j].tick_params(axis='x', labelsize=12)
        axes[i][j].tick_params(axis='y', labelsize=12)
        axes[i][j].legend()

    axes[i][0].set_title(f'{feature} Target Distribution in Training Set', fontsize=13)
    axes[i][1].set_title(f'{feature} Training & Test Set Distribution', fontsize=13)

plt.show()

災害ツイートと災害以外のツイート、訓練データとテストデータで分布に大きな違いが無いことがわかります。
URL、ハッシュタグ、メンションが含まれるツイートも、それらを含まないツイートと同程度かそれ以上あることが分かります。
URL、ハッシュタグ、メンション等の表記は災害ツイート、災害以外のツイートの判断に必要な情報でない可能性が高いので、前処理でクレンジングすると良いかもしれません。

次に、keywordの単語のうち、どれが災害ツイートで多く登場し、どれが災害以外ツイートで多く登場するかを調べます。
災害ツイートのtargetは整数で1、災害以外のツイートは整数で0のため、
keywordの単語毎のtargetの平均値を取り、1に近ければ災害ツイートに登場する傾向がある単語、
0に近ければ災害以外のツイートに登場する傾向がある単語であることが分かります。
pandasのgroupbyメソッドを使用して、keywordの単語毎にtargetの平均値を求めて、その値を訓練データ全体に付加します。
その後、災害ツイートで登場する単語から順に登場回数をプロットします。

# keywordの単語毎に、targetの平均値を求めて、その値を訓練データ全体に付加する
df_train['target_mean'] = df_train.groupby('keyword')['target'].transform('mean')

fig = plt.figure(figsize=(8, 72), dpi=100)

# keyword に含まれるラベル分布を確認
sns.countplot(y=df_train.sort_values(by='target_mean', ascending=False)['keyword'],
             hue=df_train.sort_values(by='target_mean', ascending=False)['target'])

plt.tick_params(axis='x', labelsize=15)
plt.tick_params(axis='y', labelsize=12)
plt.legend(loc=1)
plt.title('Target Distribution in Keywords')

plt.show()

# targetの値の平均値のカラムは以降使用しないので削除する
df_train.drop(columns=['target_mean'], inplace=True)

出力結果は以下のようになります。
221個の単語のラベル分布全てを、本記事に表示すると縦に長くなってしまうので、
災害ツイートに登場する傾向がある単語と、災害以外のツイートに登場する傾向がある単語、
それぞれ上位のみを表示します(実際は縦長のラベル分布の図が得られます)。

derailment(=脱線)、debris(=残骸)、wreckage(=残骸)といった災害に関連する具体的な状態を表す名詞は、災害ツイートで登場する傾向があります。
一方で、aftershock(=余震、余波)、body bags(=遺体袋)、ruin(=破滅(名詞)、台無しにする(動詞))などは一見災害に関連する単語ですが、
災害ツイートでは登場しない傾向があります。
これは、比喩表現としても使われる単語であるためだと考えられます。

次に、n-gramによる頻出単語の確認を行います。
今回は、ユニグラム(n=1)、バイグラム(n=2)、トライグラム(n=3)で頻出単語を、災害ツイート=1⇄災害以外のツイート=0でそれぞれ確認します。
まずは、n-gramのリストを生成する関数を定義します。
n-gramのリストを生成する関数の詳しい解説記事も良かったらご覧ください。

def generate_ngrams(text, n_gram=1):
    # ストップワードのリストにない単語のみをトークナイズ
    token = [token for token in text.lower().split(' ') if token != '' if token not in STOPWORDS]
    # n_gramのタプルの作成、zip(*)によってリストの先頭から同じインデックスの要素を取り出す。
    ngrams = zip(*[token[i:] for i in range(n_gram)])
    return [' '.join(ngram) for ngram in ngrams]

次に、関数を用いてユニグラム(n=1)、バイグラム(n=2)、トライグラム(n=3)とその頻度を計算します。

# ユニグラム
disaster_unigrams = defaultdict(int)
nondisaster_unigrams = defaultdict(int)

# df_trainの中の災害ツイートのユニグラムを作成する
for tweet in df_train[DISASTER_TWEETS]['text']:
    for word in generate_ngrams(tweet):
        disaster_unigrams[word] += 1

# df_trainの中の災害以外のツイートのユニグラムを作成する。
for tweet in df_train[~DISASTER_TWEETS]['text']:
    for word in generate_ngrams(tweet):
        nondisaster_unigrams[word] += 1

# 出現頻度でソートする
df_disaster_unigrams = pd.DataFrame(sorted(disaster_unigrams.items(), key=lambda x: x[1])[::-1])
df_nondisaster_unigrams = pd.DataFrame(sorted(nondisaster_unigrams.items(), key=lambda x: x[1])[::-1])

# バイグラム
disaster_bigrams = defaultdict(int)
nondisaster_bigrams = defaultdict(int)

for tweet in df_train[DISASTER_TWEETS]['text']:
    for word in generate_ngrams(tweet, n_gram=2):
        disaster_bigrams[word] += 1

for tweet in df_train[~DISASTER_TWEETS]['text']:
    for word in generate_ngrams(tweet, n_gram=2):
        nondisaster_bigrams[word] += 1

df_disaster_bigrams = pd.DataFrame(sorted(disaster_bigrams.items(), key=lambda x: x[1])[::-1])
df_nondisaster_bigrams = pd.DataFrame(sorted(nondisaster_bigrams.items(), key=lambda x: x[1])[::-1])

# トライグラム
disaster_trigrams = defaultdict(int)
nondisaster_trigrams = defaultdict(int)

for tweet in df_train[DISASTER_TWEETS]['text']:
    for word in generate_ngrams(tweet, n_gram=3):
        disaster_trigrams[word] += 1

for tweet in df_train[~DISASTER_TWEETS]['text']:
    for word in generate_ngrams(tweet, n_gram=3):
        nondisaster_trigrams[word] += 1

df_disaster_trigrams = pd.DataFrame(sorted(disaster_trigrams.items(), key=lambda x: x[1])[::-1])
df_nondisaster_trigrams = pd.DataFrame(sorted(nondisaster_trigrams.items(), key=lambda x: x[1])[::-1])

まずは、出現頻度の高い30のユニグラムを見てみましょう。

N = 30 # 上位30のユニグラムだけ表示する

fig, axes = plt.subplots(ncols=2, figsize=(15, 15), dpi=100)
plt.tight_layout()

sns.barplot(y=df_disaster_unigrams[0].values[:N], x=df_disaster_unigrams[1].values[:N], ax=axes[0], color='red')
sns.barplot(y=df_nondisaster_unigrams[0].values[:N], x=df_nondisaster_unigrams[1].values[:N], ax=axes[1], color='green')

for i in range(2):
    axes[i].spines['right'].set_visible(False)
    axes[i].set_xlabel('')
    axes[i].set_ylabel('')
    axes[i].tick_params(axis='x', labelsize=13)
    axes[i].tick_params(axis='y', labelsize=13)

axes[0].set_title(f'Top {N} most common unigrams in Disaster Tweets', fontsize=15)
axes[1].set_title(f'Top {N} most common unigrams in Non-disaster Tweets', fontsize=15)

plt.show()

災害ツイート、災害以外のツイートであっても、最頻出のユニグラムの多くが記号、取り除ききれなかったストップワード、数字であることが分かります。
これらのユニグラムはtargetの判断基準にならないので、モデリングの前に削除した方がよいでしょう。

また、災害ツイートで頻出のユニグラムは、災害に関する具体的な情報を提供していることが分かります。
一方、災害以外のツイートで頻出のユニグラムの中に、動詞が多いことが分かります。
これは、災害以外のツイートでは、ユーザーが自身もしくは何かの動作についてツイートする傾向があるからだと考えられます。

バイグラムとトライグラムについても見てみましょう。

# バイグラム
fig, axes = plt.subplots(ncols=2, figsize=(20, 15), dpi=100)
plt.subplots_adjust(wspace=0.4, hspace=0.6)

sns.barplot(y=df_disaster_bigrams[0].values[:N], x=df_disaster_bigrams[1].values[:N], ax=axes[0], color='red')
sns.barplot(y=df_nondisaster_bigrams[0].values[:N], x=df_nondisaster_bigrams[1].values[:N], ax=axes[1], color='green')

for i in range(2):
    axes[i].spines['right'].set_visible(False)
    axes[i].set_xlabel('')
    axes[i].set_ylabel('')
    axes[i].tick_params(axis='x', labelsize=20)
    axes[i].tick_params(axis='y', labelsize=20)

axes[0].set_title(f'Top {N} most common bigrams in Disaster Tweets', fontsize=20)
axes[1].set_title(f'Top {N} most common bigrams in Non-disaster Tweets', fontsize=20)

plt.show()

# トライグラム
fig, axes = plt.subplots(ncols=2, figsize=(20, 15), dpi=100)
plt.subplots_adjust(wspace=0.7, hspace=0.6)

sns.barplot(y=df_disaster_trigrams[0].values[:N], x=df_disaster_trigrams[1].values[:N], ax=axes[0], color='red')
sns.barplot(y=df_nondisaster_trigrams[0].values[:N], x=df_nondisaster_trigrams[1].values[:N], ax=axes[1], color='green')

for i in range(2):
    axes[i].spines['right'].set_visible(False)
    axes[i].set_xlabel('')
    axes[i].set_ylabel('')
    axes[i].tick_params(axis='x', labelsize=20)
    axes[i].tick_params(axis='y', labelsize=20)

axes[0].set_title(f'Top {N} most common trigrams in Disaster Tweets', fontsize=20)
axes[1].set_title(f'Top {N} most common trigrams in Non-disaster Tweets', fontsize=20)

plt.show()

災害ツイートは、バイグラムもトライグラムも共通して具体的な災害の内容を多く含んでいることが分かります。
また、災害ツイートではユニグラムにあったような記号やストップワード、数字があまり登場しないことも分かります。
一方、災害以外のツイートは、区切り文字やストップワードが登場し、redditやyoutubeといった単語も多く登場します。

4.データの前処理

探索的データ分析によって、ツイートからモデル構築に不要な情報を取り除く必要があることが分かりました。
モデル構築の前にいくつかの前処理を訓練データ、テストデータに行います。

まず、I'mや、we'veといった短縮形の単語をI amやwe haveといった形に戻します。
pythonではcontractionsというモジュールが提供されており、contractions.fix(text)で、
短縮形を元の形に戻すことができます。

def fix_contractions(text):
    return contractions.fix(text)

# 関数適応前のツイート例
print("tweet before contractions fix : ", df_train.iloc[1055]["text"])

# 関数の適用
df_train['text']=df_train['text'].apply(lambda x : fix_contractions(x))
df_test['text']=df_test['text'].apply(lambda x : fix_contractions(x))

# 関数適用後のツイート例
print("tweet after contractions fix : ", df_train.iloc[1055]["text"])

実行結果

tweet before contractions fix :  @asymbina @tithenai I'm hampered by only liking cross-body bags. I really like Ella Vickers bags: machine washable. http://t.co/YsFYEahpVg
tweet after contractions fix :  @asymbina @tithenai I am hampered by only liking cross-body bags. I really like Ella Vickers bags: machine washable. http://t.co/YsFYEahpVg

次に、URLを含むツイートから、正規表現を用いてURLのみ削除します。

def remove_URL(text):
    url = re.compile(r'https?://\S+|www\.\S+')
    return url.sub(r'',text)

# 関数適応前のツイート例
print("tweet before URL removal : ", df_train.iloc[1055]["text"])

# 関数の適用
df_train['text']=df_train['text'].apply(lambda x : remove_URL(x))
df_test['text']=df_test['text'].apply(lambda x : remove_URL(x))

# 関数適用後のツイート例
print("tweet after URL removal : ", df_train.iloc[1055]["text"])

実行結果

tweet before URL removal :  @asymbina @tithenai I am hampered by only liking cross-body bags. I really like Ella Vickers bags: machine washable. http://t.co/YsFYEahpVg
tweet after URL removal :  @asymbina @tithenai I am hampered by only liking cross-body bags. I really like Ella Vickers bags: machine washable. 

次に、記号を削除したいと思います。
これによって、ハッシュタグやメンションを表す記号を含む#@!"$%&\'()*+,-./:;<=>?[\\]^_`{|}~の記号が削除されます。
削除する記号の一覧はstring.punctuationで取得可能です。

def remove_punct(text):
    table=str.maketrans('','',string.punctuation)
    return text.translate(table)

# 関数適応前のツイート例
print("tweet before punctuation removal : ", df_train.iloc[1055]["text"])

# 関数の適用
df_train['text']=df_train['text'].apply(lambda x : remove_punct(x))
df_test['text']=df_test['text'].apply(lambda x : remove_punct(x))

# 関数適用後のツイート例
print("tweet after punctuation removal : ", df_train.iloc[1055]["text"])

実行結果

tweet before punctuation removal :  @asymbina @tithenai I am hampered by only liking cross-body bags. I really like Ella Vickers bags: machine washable. 
tweet after punctuation removal :  asymbina tithenai I am hampered by only liking crossbody bags I really like Ella Vickers bags machine washable 

上述の3つの処理によってツイート本文が綺麗になりました。
今後は、綺麗になったテキストを用いてベクトル化、モデル化を行っていきたいと思います！

ベクトル化、モデル化の記事は現在執筆中です！お待ちください！

参考にしたノートブック

本記事で参考に(or 使用)したコンペのノートブック

・https://www.kaggle.com/gunesevitan/nlp-with-disaster-tweets-eda-cleaning-and-bert
(EDA、データ前処理、ベクトル化、BERT)

・https://www.kaggle.com/shahules/basic-eda-cleaning-and-glove
(EDA、データ前処理、GloVe)

次章を進める前に行う前処理

以下のコードを実行した状態で次章に進んでください

# データ分析
import pandas as pd
import numpy as np

# 可視化
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns

# 自然言語処理
import string # 記号(punctuation)の一覧を取得
import re
import contractions # 文書の短縮形を
from wordcloud import STOPWORDS # ストップワードのリストを取得する
from collections import defaultdict # n-gram作成時に利用

# 短縮形の復元
def fix_contractions(text):
    return contractions.fix(text)

# URLの削除
def remove_URL(text):
    url = re.compile(r'https?://\S+|www\.\S+')
    return url.sub(r'',text)

# 記号の削除
def remove_punct(text):
    table=str.maketrans('','',string.punctuation)
    return text.translate(table)

# 関数の適用
df_train['text']=df_train['text'].apply(lambda x : fix_contractions(x))
df_test['text']=df_test['text'].apply(lambda x : fix_contractions(x))

df_train['text']=df_train['text'].apply(lambda x : remove_URL(x))
df_test['text']=df_test['text'].apply(lambda x : remove_URL(x))

df_train['text']=df_train['text'].apply(lambda x : remove_punct(x))
df_test['text']=df_test['text'].apply(lambda x : remove_punct(x))

自分の低クオリティツイートを抽出したかった

• Twitterには「リプライした先」のIDであるin_reply_to_status_idが存在するが「リプライを受けている」といったフラグは存在しないようだ
  • ちなみにリプライではないツイートはin_reply_to_status_idはNoneとなっていた
• 最初はいいね、RTされていないツイートを抽出して、検索を行って抽出していたが、検索APIの上限にすぐに引っかかってしまった
• そこで、以下の流れにした
  • いいね、RTされていないツイートを抽出
  • 直近のリプライを抽出
  • いいね、RTされていないツイートと直近のリプライを突合してリプライを受けるツイートを除外
  • いいね、RT、リプライを受けていないツイートのIDリストが完成

コード

import tweepy

consumer_key = '自分の'
consumer_secret = '自分の'
access_key = '自分の'
access_secret = '自分の'
auth = tweepy.OAuthHandler(consumer_key, consumer_secret)
auth.set_access_token(access_key, access_secret)
api = tweepy.API(auth)

noFavRtTweet = []
mentionTweet = []
noDeleteTweet = []

# No Favorite No RT Tweet
for tweet in tweepy.Cursor(api.user_timeline,screen_name = "自分のID",exclude_replies = True).items(200):
    if tweet.favorite_count == 0 and tweet.retweet_count == 0:
        print(tweet.id,tweet.created_at,tweet.text)
        noFavRtTweet.append(tweet.id)

# Reply Tweet
for mentions in tweepy.Cursor(api.mentions_timeline).items(200):
    print(mentions.id,mentions.created_at,mentions.text)
    mentionTweet.append(mentions.in_reply_to_status_id)

# No Favorite No RT No Reply Tweet
for tweet in noFavRtTweet:
    for reptw in mentionTweet:
        if tweet == reptw:
            print(api.get_status(tweet).id,api.get_status(tweet).created_at,api.get_status(tweet).text)
            noDeleteTweet.append(tweet)

# Extraction Delete Tweet
perfectList = set(noFavRtTweet) ^ set(noDeleteTweet)
print(list(perfectList))

# Delete Tweet
for tweet in perfectList:
    print(api.get_status(tweet).id,api.get_status(tweet).created_at,api.get_status(tweet).text)
# api.destroy_status(deltw) 自分はLambdaを使って日次で消している

使い道

自分は上記スクリプトをLambdaを使って日次で消している。次はLambdaで日次バッチ的に自分の低クオリティツイートを削除するCloudFormationを投稿したい。実際無料枠で十分に消すことが可能。

読者の想定

初心者向けとありますが、私も初心者です。
webスクレイピングの簡単なサンプルコードの理解をした後、ちょっと自身のオリジナリティを出したいなと思って、調べながらやってみたという内容です。
webスクレイピングの参考コード通りに実行したら、タイトルとかを抜きだす事ができた！がレベル1なら、今回はレベル2くらいの内容だと思います。
なので、勘違いの部分もあったりすると思うので、ご指摘あればコメントいただければと思います。

はじめに

環境

python 3.7.3
私は、visual studio codeで開発しました。

ライブラリのインポート

Pythonには"urlib2"というHTTPのライブラリがあるのですが、使い勝手がよくないので、"Requests"と"BeautifulSoup"のライブラリを使ってwebスクレイピングを行います。
RequestsでWebページを取得して、Beautiful SoupでそのHTMLを抽出します。

やっていこう！

スクレイピング内容

日経ビジネス電子版
https://business.nikkei.com/
から、新着記事の見出しとURLを取得してみようと思います。

Google Chromeでアクセスし、F12キーで、デベロッパーツール（検証モード）にアクセスします。

新着記事の部分が、HTMLのどの部分なのか知りたいので、Ctrl + Shift + Cで見出しにカーソルを動かします。

すると、classがcategoryという部分に、記事の連載名がある事が分かりました。

記事の見出しは、h3タグに入っています。
また、URLは、少し上のaタグの部分にある事が分かります。
この関係を図に示すと、以下のような構図になります。
あとで、プログラムと一緒に解説したいと思います。

コードの解説

code.py

import requests
from bs4 import BeautifulSoup
import re

urlName = "https://business.nikkei.com"
url = requests.get(urlName)
soup = BeautifulSoup(url.content, "html.parser")

requestsのライブラリでhttp接続をして、BeautifulSoupでhtmlの解析をします。

code.py

elems = soup.find_all("span")

まず、span要素を全てelemsに格納します。

code.py

for elem in elems: 
  try:
    string = elem.get("class").pop(0)
    if string in "category":
      print(elem.string)
      title = elem.find_next_sibling("h3")
      print(title.text.replace('\n',''))
      r = elem.find_previous('a')
      print(urlName + r.get('href'), '\n')
  except:
    pass

次に、span要素の中から、class名を取り出して、categoryかどうかを判別します。
もし、classがcategoryであった場合には、連載名のテキストを、.stringを用いて抜き出します。

そして、次は見出しの内容を取得します。
見出しは、h3タグにありました。
h3タグは、同じ深さで、すぐ下の部分にありました。
なので、find_next_sibling()を使って、要素以降の同じ深さにあるh3を検索します。

抜き出した文章には、画像もついている場合があり、改行が含まれている場合と含まれていない場合があったので、含まれている場合は削除しました。

最後に、URLを抜き出しにいきたいと思います。先ほどは、同じ深さでしたが、aタグは、ひとつ上の深さにいます。
したがって、find_previous()を使って、aタグを探しに行き、要素の指定の属性値を取得するgetメソッドを用いて、hrefのアドレスを取り出しました。

以下、実行結果の一部です。

池松由香のニューヨーク発直行便
米海軍の巨大病院船がNY入り　それでも足らない病床
https://business.nikkei.com/atcl/gen/19/00119/033100011/

市嶋洋平のシリコンバレーインサイ…
「需要2割経済」を生きる　ポストコロナ考え動く米外食業界
https://business.nikkei.com/atcl/gen/19/00137/033100002/

橋本宗明が医薬・医療の先を読む
中国平安保険と提携した塩野義・手代木社長の深謀遠慮
https://business.nikkei.com/atcl/gen/19/00110/033100012/ 

このように、取得する事ができました。

さいごに

まだまだ勉強中なので、勘違いしている部分があったり、もっと良い方法があるかとは思います。
少しずつ理解を深めながら実践していきたいと思います。

matplotlibで画像を表示させた時にちょっと気になった点がいくつかあったので，それをまとめておきました．
きっと気になっている人もいると思うので(特に使い始めの方)，ちょっとでも検索の手間が省けたら良いなぁと思ってます．
私はjupyter notebookをよく使うので，そこでの使用を前提としてます．

１．あれ？グレースケールのはずなのに色が付いてるな？

百聞は一見に如かず．とりあえず画像を表示させてみます．

# いつもの色々インポート
import cv2
import numpy as np
from matplotlib import pyplot as plt

im = np.array(cv2.imread('lena.png',0)) #画像を読み込んでnumpy形式にする

cv2.imread('lena.png',0) は"lena.png"を 0 すなわちグレースケールで読み込んでね，ということです．
ではmatplotlibで画像を表示させてみます．

# matplotlibで表示させます
plt.imshow(im)
plt.show()

するとこんな感じで表示されます．

なんだか色が付いちゃってます．
matplotlibはデフォルトでカラーマップが設定されているのでこんな感じになってしまうんですね．
これを解除するためにはplt.gray()を入れてあげるだけです．

plt.gray() # カラーマップを"gray"に設定
plt.imshow(im)
plt.show()

ちゃんと期待通りになりました．
実はここまではいろんなところに情報があるんですけど，私はこう思ってしまったことがあるんです．

でもなんかいい感じだったから元のカラーマップに戻したいな？

案外見つからなくて，「もとに戻したいと思ってる人いないのかな…」って思ってました．
デフォルトのカラーマップが何なのかわからなかったので，元に戻すのを諦めてましたが，公式にちゃんと書いてありました．

このviridisがデフォルトのカラーマップのようです．
早速設定してみましょう．

plt.viridis()
plt.imshow(im)
plt.show()

これで元に戻りました．めでたしめでたし．

２．あれ？コントラストがおかしいな？

これに気づけたのは奇跡だと思っているんですが，なんかちょっと色がおかしいなって思ったんです．
なので，テストでこんな画像を読み込ませてみました．

まっくろです．…が，matplotlibで表示させると

どうやら色の階調が正規化されているみたいです．これを解消するためには下記のようにvminとvmaxを指定してあげないといけないようです．

plt.imshow(im, vmin = 0, vmax = 255) 

目を凝らすとわずかに見えます．

きっと他にも私が気づいていないデフォルト設定があるんでしょうね．
以上です．間違いなどありましたら教えて下さい．

matplotlibで画像を表示させた時にちょっと気になった点がいくつかあったので，それをまとめておきました．
きっと気になっている人もいると思うので(特に使い始めの方)，ちょっとでも検索の手間が省けたら良いなぁと思ってます．
私はjupyter notebookをよく使うので，そこでの使用を前提としてます．

１．あれ？グレースケールのはずなのに色が付いてるな？

百聞は一見に如かず．とりあえず画像を表示させてみます．

# いつもの色々インポート
import cv2
import numpy as np
from matplotlib import pyplot as plt

im = np.array(cv2.imread('lena.png',0)) #画像を読み込んでnumpy形式にする

cv2.imread('lena.png',0) は"lena.png"を 0 すなわちグレースケールで読み込んでね，ということです．
ではmatplotlibで画像を表示させてみます．

# matplotlibで表示させます
plt.imshow(im)
plt.show()

するとこんな感じで表示されます．

なんだか色が付いちゃってます．
matplotlibはデフォルトでカラーマップが設定されているのでこんな感じになってしまうんですね．
これを解除するためにはplt.gray()を入れてあげるだけです．

plt.gray() # カラーマップを"gray"に設定
plt.imshow(im)
plt.show()

ちゃんと期待通りになりました．
実はここまではいろんなところに情報があるんですけど，私はこう思ってしまったことがあるんです．

でもなんかいい感じだったから元のカラーマップに戻したいな？

案外見つからなくて，「もとに戻したいと思ってる人いないのかな…」って思ってました．
デフォルトのカラーマップが何なのかわからなかったので，元に戻すのを諦めてましたが，公式にちゃんと書いてありました．

このviridisがデフォルトのカラーマップのようです．
早速設定してみましょう．

plt.viridis()
plt.imshow(im)
plt.show()

これで元に戻りました．めでたしめでたし．

２．あれ？コントラストがおかしいな？

これに気づけたのは奇跡だと思っているんですが，なんかちょっと色がおかしいなって思ったんです．
なので，テストでこんな画像を読み込ませてみました．

まっくろです．…が，matplotlibで表示させると

どうやら色の階調が正規化されているみたいです．これを解消するためには下記のようにvminとvmaxを指定してあげないといけないようです．

plt.imshow(im, vmin = 0, vmax = 255) 

目を凝らすとわずかに見えます．

きっと他にも私が気づいていないデフォルト設定があるんでしょうね．
以上です．間違いなどありましたら教えて下さい．

https://qiita.com/yohiro/items/04984927d0b455700cd1
の続き

matplotlibとnumpy

matplotlibは、グラフ描画ライブラリ。pyplotをよく使うっぽい。
numpyは、数値計算ライブラリ。
この二つを組み合わせて数学グラフの描画ができるようになる。

インポートの仕方

import matplotlib.pyplot as plt
import numpy as np

グラフ描画の仕方

$y = 2x + 1$であれば、

x = np.arange(0, 10, 0.1) # 0から10まで0.1刻みのリストxを作成
y = 2*x + 1 # x各要素の2x+1の写像としてyを作成

plt.plot(x, y) # (x, y)のグラフを描画
plt.show() # グラフ表示

と書くことで以下のグラフを描ける

シグモイド関数

以下のような曲線を描く関数のこと。
式は$y=\frac{1}{1 + e^x}$で表される。

ニューロンの「入力が閾値を超えたら発火する」という性質と親和性が高く、活性化関数としてよく使われる（と思われる）。
Pythonによる実装例は下記の通り。

import math
def sigmoid(a):
    return (1.0 / (1.0 + math.exp(-a)))

今回は、対数プロットすることにより各国の感染ピークの時期を予測することとする。

やったこと

・ちょっと理論
・各国の予測
・コード解説
・追加の理論
・各国の予測II

・ちょっと理論

前回、以下の式を導出した。

I =  I_0\exp(\gamma (\frac{\beta}{\gamma} \frac{S}{N} - 1)t)\\

この両辺を対数取ると、

\ln I =  \ln I_0+\gamma (\frac{\beta}{\gamma} \frac{S}{N} - 1)t\\

が得られる。すなわち、縦軸に$\ln I$横軸に時間$t$を取ったグラフを描くと、傾きが

\gamma (\frac{\beta}{\gamma} \frac{S}{N} - 1)\\

の直線が描ける。ここで、この傾きは$S$に依存しており時間と共に減少する未感染者数である。
そして、この傾きは前回導入した実効再生産数$R$を使うと以下のように書き換えられる。

\gamma (R - 1)\\
R = R_0 \frac{S}{N}\\
R_0 = \frac{\beta}{\gamma}\\

したがって、この傾きはR=1で０になり、すなわち再生産数＝０となり、これ以上感染者が増えないので感染ピークとなる。
今回は、この理論を根拠に各国のデータを$\log I\ vs\ t$でプロットし、それぞれの国の感染ピークの時期を見たいと思う。

また、あわよくば立ち上がりの傾きから、$S\fallingdotseq N$であり、$\gamma (R_0 -1)$が得られる。
これを初期値にして得られたグラフをフィッティングすると、$\gamma$,$\beta$,$R_0$が求められそうである。

・各国の予測

・Korea, South, Diamond Princess, Italyの状況

Italyが入っているのが不信に思われるかもしれませんが、結果は以下のとおりです。

この3つの国を並べてみると、割と様になっているのが分かります。
すなわち、適当にx,y軸を平行移動すると重なりそうです。つまり、傾きが同一な領域があることが分かります。
どのグラフも、ほぼ10日で10倍になる区間が見えます。
そうです、DPが最初に感染数曲線が飽和し、次に韓国、そしてこれからイタリアがまさしく飽和しそうな状況です。下のグラフは傾きをプロットしています。すなわち、０になると飽和したことを表します。
ということで、このグラフからも上記のことが分かります。
イタリアに関して言及すると、いろいろ言われてきましたがあと10日程度でピークに達すると思われます。

・ヨーロッパ各国の状況

イタリア以外の、Spain,Germany,Switzerland,United Kingdom, Netherlandsについて、同一のグラフで表示します。

重ねて書いてみると、イギリスが少し遅れていますが、その他はよく似ている曲線を描いています。
縦軸はそれぞれの事情に依存して（人口など）異なりますが、ほぼ同様な傾きの変化をしています。
そして、傾きの図を見ると2週間程度内にはまとめて飽和してきそうです。

・ちょっと残念な国

Japan, Iran, Bahrainの3つの国もある意味曲線が酷似しています。

立ち上がりの状況が異なるので、y軸が全く違うので医療体制の状況は異なると思います。
しかし、一旦立ち上がってしまったあとの傾きが割と似ており、しかもこの何日間でほぼ飽和しているにも関わらず、イランと日本は再度スーと増加してしまっています。バーレーンも傾きが負になっているのですが、再度正に戻ってしまっています。
これら3国の勾配は、安定してからの勾配として10倍／30日程度のいつでも飽和しそうな緩い勾配になっています。
これらの国はまだまだ感染が継続しており、緊張を持って社会全体で封じ込めて行く必要があります。

・噂のアメリカの状況

SwedenとUSを同一プロットしてみました。

一目瞭然ですが、立ち上がりは似ています。しかし、スウェーデンは途中から傾きが小さくなっています。しかし、下の傾きのグラフを見るとこのところほぼx軸に平行になっていて、0に向かっているようには見えません。つまり、近日中に飽和はしないように見えます。
ただし、58日から65日のデータがそろってきていて、傾きが0に向かって急激に近づいているようにも見えます。これを伸ばすとあと10日程度で感染ピークになる可能性もあることが分かります。
なお、両国の最大勾配は25倍／20日程度で先ほどの勾配より大きいのが分かります。

・コード解説

利用するLibとデータ読み込みは前回と同一です。
死亡数データも処理していますが、以下では使わないので削除しています。

import numpy as np
from scipy.integrate import odeint
from scipy.optimize import minimize
import matplotlib.pyplot as plt
import pandas as pd

#pandasでCSVデータ読む。
data = pd.read_csv('COVID-19/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv')
data_r = pd.read_csv('COVID-19/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_recovered_global.csv')

以下、データ加工関数で、処理した感染数データと勾配データを返します。

def data_city(city):
    #データを加工する
    t_cases = 0
    t_recover = 0
    t_deaths = 0
    for i in range(0, len(data_r), 1):
        if (data_r.iloc[i][1] == city): #for country/region
        #if (data_r.iloc[i][0] == city):  #for province:/state  
            print(str(data_r.iloc[i][0]) + " of " + data_r.iloc[i][1])
            for day in range(4, len(data.columns), 1):            
                confirmed_r[day - 4] += data_r.iloc[i][day]
            t_recover += data_r.iloc[i][day]      
    for i in range(0, len(data), 1):
        if (data.iloc[i][1] == city): #for country/region
        #if (data.iloc[i][0] == city):  #for province:/state  
            print(str(data.iloc[i][0]) + " of " + data.iloc[i][1])
            for day in range(4, len(data.columns), 1):
                confirmed[day - 4] += data.iloc[i][day] -  confirmed_r[day - 4]
                diff_confirmed[day - 4] += confirmed[day-4] /  confirmed_r[day - 4]

    tl_confirmed = 0
    dlog_confirmed = [0] * (len(data.columns) - 4)
    dlog_confirmed[0]=np.log(confirmed[0])
    dlog_confirmed[1]=np.log(confirmed[1])-np.log(confirmed[0])
    for i in range(2, len(confirmed)-2, 1):
        if confirmed[i] > 0:
            dlog_confirmed[i]=(np.log(confirmed[i+1])-np.log(confirmed[i-1]))/2
        else:
            continue
    tl_confirmed = confirmed[len(confirmed)-1] + confirmed_r[len(confirmed)-1] + confirmed_d[len(confirmed)-1]
    t_cases = tl_confirmed

    return confirmed, dlog_confirmed       

複数の国を同一図上にプロットするために以下のようにしています。
city_listに描きたい国を入れます。
上記のデータ処理関数で感染数データとその片対数の傾きを取得します。
これを描画します。
そのためにfig, (ax1,ax2) = ．．．の定義をfor文外で定義しました。
最後にまとめて、表示します。
colorはあえて定義せず、自動で選んでもらっています。

city_list={ "Spain","Germany","Switzerland","United Kingdom", "Netherlands" }
fig, (ax1,ax2) = plt.subplots(2,1,figsize=(1.6180 * 4, 4*2))
for city in city_list:
    confirmed = [0] * (len(data.columns) - 4)
    confirmed_r = [0] * (len(data_r.columns) - 4)
    confirmed_d = [0] * (len(data_d.columns) - 4)
    diff_confirmed = [0] * (len(data.columns) - 4)
    days_from_22_Jan_20 = np.arange(0, len(data.columns) - 4, 1)
    days_from_22_Jan_20_ = np.arange(0, len(data.columns) - 4, 1)
    recovered_rate = [0] * (len(data_r.columns) - 4)
    deaths_rate = [0] * (len(data_d.columns) - 4)

    confirmed, dlog_confirmed = data_city(city)
    #matplotlib描画
    lns1=ax1.semilogy(days_from_22_Jan_20, confirmed, ".-", label = str(city))
    lns4=ax2.plot(days_from_22_Jan_20_, dlog_confirmed, ".-",label = str(city))

ax1.set_xlabel("days from 22, Jan, 2020")
ax1.set_ylabel("casas ")
ax2.set_ylabel("slopes ")
ax1.legend(loc=0)
ax2.legend(loc=0)
ax1.grid()
ax2.grid()

plt.pause(1)
plt.savefig('./fig/logmulti_{}_.png'.format(city)) 
plt.close()

・追加の理論

実は、治癒数曲線の推移は回復に遅れはあるが、感染数曲線と類似しているだろうという想像ができる。
ということで、以下その話を眉唾だがやってみた。
まず、SIRモデルの第三方程式に$I$を代入すると以下のとおりとなる。
※再生産数のRと回復数のRが被るのでここでは回復数を小文字にした

\frac{dr}{dt} =  \gamma I_0\exp(\gamma (R - 1)t) \\

従って、両辺の対数を取ると

\ln (\frac{dr}{dt}) =  \ln \gamma + \ln I_0 + \gamma (R - 1)t \\

となる。すなわち、回復率の時間微分（傾き）は上式と傾きが同じ曲線になることが期待でき、y切片が$\ln \gamma$だけ異なるものが得られる。
しかし、これだと誤差も大きく扱い難そうなので、ちょっと変更します。

すなわち、（解析的には正しくないが）Rがコンスタントとして積分すると以下の式を得る。

r =  \frac{I_0}{(R-1)}\exp(\gamma (R - 1)t) + r_0\\

したがって、両辺対数取ると

\ln r =  \ln I_0 -\ln (R-1) + \gamma (R - 1)t +\ln r_0

と書ける。ということで上記の$\ln I$の式と比較するとほぼ平行でｙ切片が$-ln (R-1)+\ln r_0$異なる曲線が得られることが期待できる。

・各国の予測II

・武漢の状況

以下のように綺麗な曲線が得られました。予想どおり感染数曲線と治癒数曲線は平行に始まり感染数ピークを過ぎたところで交差します。
それに合わせて、二つの曲線の差分（$\log (A/B)=\log A-\log B$）は直線的に$10^0$を通過します。

・イタリアの状況

両曲線が平行になっており、上では10日程度で感染数ピークになりそうという結論を出しましたが、そのあと、なかなか交差しないように見えます。
このあと、きっと近づいてきて交差するものと考察されますが、治癒数データが怪しいかもしれません。

・ヨーロッパを代表してスペインの状況

着実に両曲線が近づきつつあり、感染数ピークのあと交差し終息が見えてくる状況です。

・日本の状況

日本はだんだん二つの曲線が近づいて来て、差は0に向かっていましたが、このところの感染数増加で差が広がってしまいました。ということで上記と同様、感染拡大は継続しているので抑え込みが必要という結論になります。

まとめ

・各国の感染数データを対数プロットして分類し、それぞれの感染数ピークの時期を予測した
・感染数と治癒数データを同時にプロットすることにより、終息時期を予測した

・日々変わる状況を追尾したいと思う
・追加部分の理論が怪しいのでもう少し整理し、パラメータ決定につなげようと思う

今回は、対数プロットすることにより各国の感染ピークの時期を予測することとする。

やったこと

・ちょっと理論
・各国の予測
・コード解説
・追加の理論
・各国の予測II

・ちょっと理論

前回、以下の式を導出した。

I =  I_0\exp(\gamma (\frac{\beta}{\gamma} \frac{S}{N} - 1)t)\\

この両辺を対数取ると、

\ln I =  \ln I_0+\gamma (\frac{\beta}{\gamma} \frac{S}{N} - 1)t\\

が得られる。すなわち、縦軸に$\ln I$横軸に時間$t$を取ったグラフを描くと、傾きが

\gamma (\frac{\beta}{\gamma} \frac{S}{N} - 1)\\

の直線が描ける。ここで、この傾きは$S$に依存しており時間と共に減少する未感染者数である。
そして、この傾きは前回導入した実効再生産数$R$を使うと以下のように書き換えられる。

\gamma (R - 1)\\
R = R_0 \frac{S}{N}\\
R_0 = \frac{\beta}{\gamma}\\

したがって、この傾きはR=1で０になり、すなわち再生産数＝０となり、これ以上感染者が増えないので感染ピークとなる。
今回は、この理論を根拠に各国のデータを$\log I\ vs\ t$でプロットし、それぞれの国の感染ピークの時期を見たいと思う。

また、あわよくば立ち上がりの傾きから、$S\fallingdotseq N$であり、$\gamma (R_0 -1)$が得られる。
これを初期値にして得られたグラフをフィッティングすると、$\gamma$,$\beta$,$R_0$が求められそうである。

・各国の予測

・Korea, South, Diamond Princess, Italyの状況

Italyが入っているのが不信に思われるかもしれませんが、結果は以下のとおりです。

この3つの国を並べてみると、割と様になっているのが分かります。
すなわち、適当にx,y軸を平行移動すると重なりそうです。つまり、傾きが同一な領域があることが分かります。
どのグラフも、ほぼ10日で10倍になる区間が見えます。
そうです、DPが最初に感染数曲線が飽和し、次に韓国、そしてこれからイタリアがまさしく飽和しそうな状況です。下のグラフは傾きをプロットしています。すなわち、０になると飽和したことを表します。
ということで、このグラフからも上記のことが分かります。
イタリアに関して言及すると、いろいろ言われてきましたがあと10日程度でピークに達すると思われます。

・ヨーロッパ各国の状況

イタリア以外の、Spain,Germany,Switzerland,United Kingdom, Netherlandsについて、同一のグラフで表示します。

重ねて書いてみると、イギリスが少し遅れていますが、その他はよく似ている曲線を描いています。
縦軸はそれぞれの事情に依存して（人口など）異なりますが、ほぼ同様な傾きの変化をしています。
そして、傾きの図を見ると2週間程度内にはまとめて飽和してきそうです。

・ちょっと残念な国

Japan, Iran, Bahrainの3つの国もある意味曲線が酷似しています。

立ち上がりの状況が異なるので、y軸が全く違うので医療体制の状況は異なると思います。
しかし、一旦立ち上がってしまったあとの傾きが割と似ており、しかもこの何日間でほぼ飽和しているにも関わらず、イランと日本は再度スーと増加してしまっています。バーレーンも傾きが負になっているのですが、再度正に戻ってしまっています。
これら3国の勾配は、安定してからの勾配として10倍／30日程度のいつでも飽和しそうな緩い勾配になっています。
これらの国はまだまだ感染が継続しており、緊張を持って社会全体で封じ込めて行く必要があります。

・噂のアメリカの状況

SwedenとUSを同一プロットしてみました。

一目瞭然ですが、立ち上がりは似ています。しかし、スウェーデンは途中から傾きが小さくなっています。しかし、下の傾きのグラフを見るとこのところほぼx軸に平行になっていて、0に向かっているようには見えません。つまり、近日中に飽和はしないように見えます。
ただしアメリカは、58日から65日のデータがそろってきていて、傾きが0に向かって急激に近づいているようにも見えます。これを伸ばすとあと10日程度で感染ピークになる可能性もあることが分かります。
なお、両国の最大勾配は25倍／20日程度で先ほどの勾配より大きいのが分かります。

・コード解説

利用するLibとデータ読み込みは前回と同一です。
死亡数データも処理していますが、以下では使わないので削除しています。

import numpy as np
from scipy.integrate import odeint
from scipy.optimize import minimize
import matplotlib.pyplot as plt
import pandas as pd

#pandasでCSVデータ読む。
data = pd.read_csv('COVID-19/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv')
data_r = pd.read_csv('COVID-19/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_recovered_global.csv')

以下、データ加工関数で、処理した感染数データと勾配データを返します。

def data_city(city):
    #データを加工する
    t_cases = 0
    t_recover = 0
    t_deaths = 0
    for i in range(0, len(data_r), 1):
        if (data_r.iloc[i][1] == city): #for country/region
        #if (data_r.iloc[i][0] == city):  #for province:/state  
            print(str(data_r.iloc[i][0]) + " of " + data_r.iloc[i][1])
            for day in range(4, len(data.columns), 1):            
                confirmed_r[day - 4] += data_r.iloc[i][day]
            t_recover += data_r.iloc[i][day]      
    for i in range(0, len(data), 1):
        if (data.iloc[i][1] == city): #for country/region
        #if (data.iloc[i][0] == city):  #for province:/state  
            print(str(data.iloc[i][0]) + " of " + data.iloc[i][1])
            for day in range(4, len(data.columns), 1):
                confirmed[day - 4] += data.iloc[i][day] -  confirmed_r[day - 4]
                diff_confirmed[day - 4] += confirmed[day-4] /  confirmed_r[day - 4]

    tl_confirmed = 0
    dlog_confirmed = [0] * (len(data.columns) - 4)
    dlog_confirmed[0]=np.log(confirmed[0])
    dlog_confirmed[1]=np.log(confirmed[1])-np.log(confirmed[0])
    for i in range(2, len(confirmed)-2, 1):
        if confirmed[i] > 0:
            dlog_confirmed[i]=(np.log(confirmed[i+1])-np.log(confirmed[i-1]))/2
        else:
            continue
    tl_confirmed = confirmed[len(confirmed)-1] + confirmed_r[len(confirmed)-1] + confirmed_d[len(confirmed)-1]
    t_cases = tl_confirmed

    return confirmed, dlog_confirmed       

複数の国を同一図上にプロットするために以下のようにしています。
city_listに描きたい国を入れます。
上記のデータ処理関数で感染数データとその片対数の傾きを取得します。
これを描画します。
そのためにfig, (ax1,ax2) = ．．．の定義をfor文外で定義しました。
最後にまとめて、表示します。
colorはあえて定義せず、自動で選んでもらっています。

city_list={ "Spain","Germany","Switzerland","United Kingdom", "Netherlands" }
fig, (ax1,ax2) = plt.subplots(2,1,figsize=(1.6180 * 4, 4*2))
for city in city_list:
    confirmed = [0] * (len(data.columns) - 4)
    confirmed_r = [0] * (len(data_r.columns) - 4)
    confirmed_d = [0] * (len(data_d.columns) - 4)
    diff_confirmed = [0] * (len(data.columns) - 4)
    days_from_22_Jan_20 = np.arange(0, len(data.columns) - 4, 1)
    days_from_22_Jan_20_ = np.arange(0, len(data.columns) - 4, 1)
    recovered_rate = [0] * (len(data_r.columns) - 4)
    deaths_rate = [0] * (len(data_d.columns) - 4)

    confirmed, dlog_confirmed = data_city(city)
    #matplotlib描画
    lns1=ax1.semilogy(days_from_22_Jan_20, confirmed, ".-", label = str(city))
    lns4=ax2.plot(days_from_22_Jan_20_, dlog_confirmed, ".-",label = str(city))

ax1.set_xlabel("days from 22, Jan, 2020")
ax1.set_ylabel("casas ")
ax2.set_ylabel("slopes ")
ax1.legend(loc=0)
ax2.legend(loc=0)
ax1.grid()
ax2.grid()

plt.pause(1)
plt.savefig('./fig/logmulti_{}_.png'.format(city)) 
plt.close()

・追加の理論

実は、治癒数曲線の推移は回復に遅れはあるが、感染数曲線と類似しているだろうという想像ができる。
ということで、以下その話を眉唾だがやってみた。
まず、SIRモデルの第三方程式に$I$を代入すると以下のとおりとなる。
※再生産数のRと回復数のRが被るのでここでは回復数を小文字にした

\frac{dr}{dt} =  \gamma I_0\exp(\gamma (R - 1)t) \\

従って、両辺の対数を取ると

\ln (\frac{dr}{dt}) =  \ln \gamma + \ln I_0 + \gamma (R - 1)t \\

となる。すなわち、回復率の時間微分（傾き）は上式と傾きが同じ曲線になることが期待でき、y切片が$\ln \gamma$だけ異なるものが得られる。
しかし、これだと誤差も大きく扱い難そうなので、ちょっと変更します。

すなわち、（解析的には正しくないが）Rがコンスタントとして積分すると以下の式を得る。

r =  \frac{I_0}{(R-1)}\exp(\gamma (R - 1)t) + r_0\\

したがって、両辺対数取ると

\ln r =  \ln I_0 -\ln (R-1) + \gamma (R - 1)t +\ln r_0

と書ける。ということで上記の$\ln I$の式と比較するとほぼ平行でｙ切片が$-ln (R-1)+\ln r_0$異なる曲線が得られることが期待できる。

・各国の予測II

・武漢の状況

以下のように綺麗な曲線が得られました。予想どおり感染数曲線と治癒数曲線は平行に始まり感染数ピークを過ぎたところで交差します。
それに合わせて、二つの曲線の差分（$\log (A/B)=\log A-\log B$）は直線的に$10^0$を通過します。

・イタリアの状況

両曲線が平行になっており、上では10日程度で感染数ピークになりそうという結論を出しましたが、そのあと、なかなか交差しないように見えます。
このあと、きっと近づいてきて交差するものと考察されますが、治癒数データが怪しいかもしれません。

・ヨーロッパを代表してスペインの状況

着実に両曲線が近づきつつあり、感染数ピークのあと交差し終息が見えてくる状況です。

・日本の状況

日本はだんだん二つの曲線が近づいて来て、差は0に向かっていましたが、このところの感染数増加で差が広がってしまいました。ということで上記と同様、感染拡大は継続しているので抑え込みが必要という結論になります。

まとめ

・各国の感染数データを対数プロットして分類し、それぞれの感染数ピークの時期を予測した
・感染数と治癒数データを同時にプロットすることにより、終息時期を予測した

・日々変わる状況を追尾したいと思う
・追加部分の理論が怪しいのでもう少し整理し、パラメータ決定につなげようと思う