LightsailのWordPressの開発環境を立てる

このBlogはWordPressの設定を最低限度にしか変更せずに運用しています。(記事執筆時時点)
ただ、それでも少々修正したい部分は発生しており、その一部はどうも管理画面では修正できず、PHPファイルを修正する必要があるようです。
(具体的に挙げと、フッターの、 「プライバシーポリシー Proudly powered by WordPress」 の部分など)

僕自身が普段はPHPを書いておらず、WordPressに不慣れなのに本番環境を直接いじって修正するのはリスクが高すぎるので、
開発環境が欲しいなと思ったので調べてみました。

このBlogはLightsailで動いているので、同じようにLightsailでWordpressインスタンスを立てて、
今入れているテーマやプラグインを追加していくのも方法の一つです。

ただ、そのようなことををしなくてもこのサーバーをそのままコピーできるらしいことがわかりました。

手順は以下の通りです。
1. Lightsail の管理画面に入り、このBlogが動いているインスタンスを選択。
2. スナップショットを選択し、手動スナップショットをとる。(自動スナップショットを取っている人はこの手順は不要)
3. 取得したスナップショットから、「新規インスタンスを作成」

これで、記事等も全部入った状態でインスタンスのコピーが立ち上がります。

あとは、
http://{IPアドレス}/
でアクセスしたら、開発環境にアクセスできて万歳、となると思っていたのですが、そうはならず、
何度試しても、本番環境(このBlog)にリダイレクトされてしまいました。

原因は以前設定したリダイレクト設定です。
参考: httpのアクセスをhttpsにリダイレクトする

開発サーバーにSSHログインし、bitnami.confに参考記事で追記した3行を消し、apacheを再起動すると、開発環境にアクセスできるようになります。

お金がかかるので、必要な修正の検討が終わったら作った開発環境とスナップショットは消しておきましょう。

Amazon Linux 2 (EC2)にphpMyAdminを導入

MySQLやその互換のDBを管理するphpMyAdminという便利なツールがあります。
僕が私用で使っている Aurora Serverless の管理もphpMyAdminで行おうとしたら、予想よりも苦戦したのでそのメモです。
参考: Amazon Aurora Serverlessを使ってみる

昔、Amazon Linux (無印)のインスタンスを建てたときは、ApachとPHPを入れて、phpMyAdminのソースファイルを配置しただけで
当時使っていた同じインスタンス内のMySQLの管理をすぐできるようになりました。
しかし当時の手順を元に、Amazon Linux 2で同様の操作を行ったら全然うまくいきませんでした。

改めて手順を探してみると、「チュートリアル: Amazon Linux 2 に LAMP ウェブサーバーをインストールする
というドキュメントが用意されており、これを参考にすることでうまくいきました。
(このチュートリアルでは、同サーバーにMariaDBを入れますが、僕はRDSのAurora Serverlessを使うので、MariaDB関係の手順は省略します。)

EC2インスタンス作成

最初にインスタンスを作ります。
DBをEC2インスタンスにインストールしている場合はそのインスタンスを使えば良いのでこの手順は不要です。

Amazon Linux 2 の AMIを選択し、通常通りインスタンスを作成します。
sshログインと、Webアクセスができるように、22番と80番のポートが開いたセキュリティグループを設定しておきます。

インスタンスができたら、sshで入って、yumをアップデートします。


$ sudo yum update -y

PHPのインストール

最初のハマりポイントがここです。
普通に、 
$ yum install -y php
とすると、バージョン 5.4.16 の非常に古いPHPが入ってしまいます。

チュートリアルに沿って、 amazon-linux-extras と言うのを使って 7系のPHPをインストールします。


$ sudo amazon-linux-extras install -y php7.2

# 入ったバージョンの確認
$ php -v
# 以下出力
# PHP 7.2.34 (cli) (built: Oct 21 2020 18:03:20) ( NTS )
# Copyright (c) 1997-2018 The PHP Group
# Zend Engine v3.2.0, Copyright (c) 1998-2018 Zend Technologies

Apacheのインストール

PHPの次は Apache のインストールと起動設定です。
Amazon Linux 2になって、無印の時とコマンドが変わっているので注意が必要です。


# Apache のインストール
$ sudo yum install -y httpd
# 起動
$ sudo systemctl start httpd
# サーバー起動時に自動的に起動するようにする
$ sudo systemctl enable httpd

以前は、起動は
$ service httpd start
で、自動起動は、
$ chkconfig httpd on
でしたね。

この時点でWebサーバーは立ち上がっているので、ブラウザを起動し、
http://{サーバーのIPアドレス}
にアクセスできることを確認します。

PHPの動作確認

ApacheでPHPが動くことを確認します。

昔は、Apacheの設定ファイル( /etc/httpd/conf/httpd.conf ) に
<FilesMatch \.php$>
SetHandler applocation/x-httpd-php
</FilesMatch>
など書いて設定しないと動かなかった覚えがあるのですが、それをしなくても動きます。

さて、動作確認に進みましょう。
ドキュメントルート(つまり、 /var/www/html)に移動し、phpinfo.phpというテキストファイルを生成します。


$ cd /var/www/html
$ vim phpinfo.php

# phpinfo.php に以下の内容を書き込む。
<?php phpinfo(); ?>

ブラウザから
http://{サーバーのIPアドレス}/phpinfo.php
にアクセスし、PHPのバージョンやシステム情報などのテーブルが表示されることが確認できたらPHPは動作しています。

phpMyAdminのインストール

ここから、phpMyAdmin本体のインストールです。

まず、phpの依存モジュールをインストールします。
(これはyumでいいらしいです。不思議です。)


# 必要な依存ファイルをインストール
$ sudo yum install php-mbstring -y
# Apache を再起動
$ sudo systemctl restart httpd
# php-fpm を再起動
$ sudo systemctl restart php-fpm

次に、phpMyAdminの本体ファイルをダウンロードして配置します。
ec2-userはドキュメントルート配下に書き込み権限を持っていないと思うので、
権限設定を適切にするか、面倒であればルートになって作業しましょう。


# ドキュメントルートに移動
$ cd /var/www/html
# phpMyAdminのファイルを取得する。
$ wget https://www.phpmyadmin.net/downloads/phpMyAdmin-latest-all-languages.tar.gz
# ファイル展開
$ mkdir phpMyAdmin && tar -xvzf phpMyAdmin-latest-all-languages.tar.gz -C phpMyAdmin --strip-components 1
# 不要ファイル削除
$ rm phpMyAdmin-latest-all-languages.tar.gz

以前は、 phpMyAdminのダウンロードページに訪問して、
最新バージョンを調べて、
wget https://files.phpmyadmin.net/phpMyAdmin/5.0.4/phpMyAdmin-5.0.4-all-languages.tar.gz
みたいに、バージョン指定して落としていたのですが、 latest で取れたんですね。
これは知りませんでした。

あとはブラウザで、
http://{IPアドレス}/phpMyAdmin/
にアクセスし、ログイン画面が表示されればインストールできています。

同サーバー内のDBの管理をするのであればこのまま使えます。

RDSへの接続を設定する

さて、冒頭の通り、僕が管理したいサーバーはAurora Serverlessなので、その設定を行います。
方法は設定ファイルにエンドポイントを指定するだけです。
まず、設定ファイルを作成します。
デフォルトでは設定ファイルは存在しておらず、config.sample.inc.phpと言うテンプレートを、
config.inc.php という有効なファイル名にコピーして使います。


# デフォルトの設定ファイルをコピーして設定ファイルを生成する
cd /var/www/html/phpMyAdmin
cp config.sample.inc.php config.inc.php

/var/www/html/phpmyadmin/config.inc.php の以下の行を書き換える
# 元の記述
$cfg['Servers'][$i]['host'] = 'localhost';
# 修正後の記述
$cfg['Servers'][$i]['host'] = '{RDSのエンドポイント}';

ログイン確認

以上の操作が全て終わったら、実際にphpMyAdminにログインして確認します。

ユーザーとパスワードは RDSのものを使います。
Aurora Serverless なので、停止状態の場合は起動に少し時間がかかりますが、数分程度待てば無事にログインできます。

gensimのモデルを保存するときのフォーマットを調べてみた

※ この記事を書いたときの僕の環境のgensimのバージョンは 3.8.0 です。

gensimでword2vecなり、トピックモデル(LDA)なり、そのための辞書なりを学習したとき、
saveメソッドで学習したモデルを保存し、同様にloadメソッドで読み込んで使うことができます。
このとき、ファイルに保存するので当然ファイル名を決めないといけないのですが、何形式で保存しているのかよくわからず拡張子に悩んでいました。

以前書いた、gensimでword2vec という記事のサンプルコードでは、ファイル形式がわからないのでとりあえず .model としています。
公式サイトの、Usage examples もそうなってるのですよ。

LDAの方は、
gensim.models.ldamodel.LdaModel.save
の Note を読むと、どうやら、pickleを使ってるように見えます。

ソースコードを確認してみましょう。
https://github.com/RaRe-Technologies/gensim/blob/master/gensim/models/ldamodel.py

def save(~略~):
のところを見ていくと、
salf.state ってのを使って、saveしていますね。

def __init__(~略~):
のところを確認すると、salf.state には、 LdaState というクラスのインスタンスが格納されており、
LdaStateは、
class LdaState(utils.SaveLoad):
とある通り、utils.SaveLoadを継承しています。
どうやらこれが保存と読み込みの本体のようです。

そのソースコードがこちらです。
https://github.com/RaRe-Technologies/gensim/blob/master/gensim/utils.py

このファイル内の、
class SaveLoad:
のところを見ていくと、普通にpickleを使ってファイルに書き出したり読み込んだ理されていることがわかりました。

LDAの次はWord2Vecです。
こちらは話が単純です。

ソースコードを見てみます。
https://github.com/RaRe-Technologies/gensim/blob/master/gensim/models/word2vec.py

class Word2Vec(utils.SaveLoad):

というふうに、モデル自体が、utils.SaveLoadを継承して作られており、先ほどのLDAと同様に保存と読み込みにはpickleが使われています。

Dictionary も同様です。
https://github.com/RaRe-Technologies/gensim/blob/master/gensim/corpora/dictionary.py

pickleは、このブログでも以前記事にしたことがあるようにPythonのオブジェクトを手軽にファイルに書き出せる形式です。
参考: pickleを使ってpythonのオブジェクトをファイルに保存する
ただ、これを使うとなると、いちいちwith openとかいろいろ書かないといけなくてややこしいので、
saveやloadなどわかりやすいメソッド名でラッピングしてもらえているのはありがたいですね。

さて、冒頭に挙げた問題の拡張子名ですが、.pickle あたりを使えば良さそうです。

gensimのDictionaryオブジェクトに含まれれる単語を出現頻度で絞り込む

最近久々にgensimのトピックモデルを使う機会がありました。
そのとき、出現する単語を頻度で絞り込みたったので方法を調べました。

トピックモデルの方法自体は、既に記事を書いてますのでこちらをご参照ください。
参考: gensimでトピックモデル(LDA)をやってみる

さて、gensimのLDAは、学習するコーパスを(単語, 出現回数) というタプルの配列に変換して読み込ませる必要があり、
その形への変換に、gensim.corpora.dictionary.Dictionaryを使います。
この辞書は、何も指定しないと、1回以上出現した単語を全部学習してしまいます。
それを、Scikit-learnのCountVectorizerで、min_dfを指定したときみたいに、n回以上出現した単語のみ、と足切りしたいというのが今回の記事の目的です。

Dictionaryの語彙学習時に指定できる引数の中に、CountVectorizerのmin_dfに相当するものがなかったので、てっきり指定できないのかと思っていたのですが、
じつは、学習した後に、語彙を絞り込む関数である、filter_extremesが用意されていることがわかりました。

使いかたを説明するために、まず適当な単語の羅列でコーパスを作って、辞書を学習しておきます。


import numpy as np
from gensim.corpora.dictionary import Dictionary

# 単語リスト作成
words = [
    "White",
    "Black",
    "Grey",
    "Red",
    "Orange",
    "Yellow",
    "Green",
    "Purple",
    "Blue",
    "Cyan",
    "Magenta",
]

# 再現性のためシードを固定する
np.random.seed(2)
# 単語を適当に選んで文章データを生成
documents = [
    np.random.choice(words, size=np.random.randint(3, 7)).tolist() for _ in range(10)
]

print(documents)
"""
[['Blue', 'Green', 'Grey'],
 ['Blue', 'Purple', 'Grey', 'Black', 'Yellow', 'Magenta'],
 ['Orange', 'Yellow', 'Purple'],
 ['Green', 'Orange', 'Magenta', 'Red', 'Purple', 'Green'],
 ['Black', 'Magenta', 'Red', 'Yellow', 'Blue'],
 ['Green', 'Red', 'Cyan'],
 ['Grey', 'White', 'Orange', 'Grey', 'Orange', 'Magenta'],
 ['Black', 'Purple', 'Blue', 'Grey', 'Magenta', 'Cyan'],
 ['Blue', 'Purple', 'Black', 'Green'],
 ['Magenta', 'Yellow', 'Cyan']]
"""

# 辞書の作成
dictionary = Dictionary(documents)

# 学習した単語リスト
for word, id_ in dictionary.token2id.items():
    print(f"id: {id_}, 単語: {word}, 出現ドキュメント数: {dictionary.dfs[id_]}, 出現回数: {dictionary.cfs[id_]}")

"""
id: 0, 単語:  Blue,    出現ドキュメント数: 5, 出現回数: 5
id: 1, 単語:  Green,   出現ドキュメント数: 4, 出現回数: 5
id: 2, 単語:  Grey,    出現ドキュメント数: 4, 出現回数: 5
id: 3, 単語:  Black,   出現ドキュメント数: 4, 出現回数: 4
id: 4, 単語:  Magenta, 出現ドキュメント数: 6, 出現回数: 6
id: 5, 単語:  Purple,  出現ドキュメント数: 5, 出現回数: 5
id: 6, 単語:  Yellow,  出現ドキュメント数: 4, 出現回数: 4
id: 7, 単語:  Orange,  出現ドキュメント数: 3, 出現回数: 4
id: 8, 単語:  Red,     出現ドキュメント数: 3, 出現回数: 3
id: 9, 単語:  Cyan,    出現ドキュメント数: 3, 出現回数: 3
id: 10, 単語: White,   出現ドキュメント数: 1, 出現回数: 1
"""

これを4個以上の文章に登場した単語だけに絞りこみたいとすると、
filter_extremes(no_below=4)
を実行すれば良いよいうに思えます。
それでやってみたのがこちら。


dictionary.filter_extremes(no_below=4)
for word, id_ in dictionary.token2id.items():
    print(f"単語: {word}, 出現ドキュメント数: {dictionary.dfs[id_]}")
"""
単語: Blue, 出現ドキュメント数: 5
単語: Green, 出現ドキュメント数: 4
単語: Grey, 出現ドキュメント数: 4
単語: Black, 出現ドキュメント数: 4
単語: Purple, 出現ドキュメント数: 5
単語: Yellow, 出現ドキュメント数: 4
"""

Orange/Red/Cyan/White が消えましたね。Orangeは出現回数自体は4でしたが、ドキュメント数が3だったので消えています。
ここで注意なのが、出現ドキュメント数が6だった、Magentaも消えていることです。

これは、filter_extremesのデフォルトの引数が、(no_below=5, no_above=0.5, keep_n=100000, keep_tokens=None) と、
no_above=0.5 も指定されていることに起因します。
つまり、全体の0.5=50%よりも多く出現している単語も一緒に消してしまうわけです。
逆に、no_above だけ指定しても、no_belowは5扱いなので、4文書以下にしか登場しない単語は足切りされます。

この挙動が困る場合は、忘れないように、no_belowとno_aboveを両方指定する必要があります。


# もう一度辞書の作成
dictionary = Dictionary(documents)
dictionary.filter_extremes(no_below=4, no_above=1)
for word, id_ in dictionary.token2id.items():
    print(f"単語: {word}, 出現ドキュメント数: {dictionary.dfs[id_]}")
"""
単語: Blue, 出現ドキュメント数: 5
単語: Green, 出現ドキュメント数: 4
単語: Grey, 出現ドキュメント数: 4
単語: Black, 出現ドキュメント数: 4
単語: Magenta, 出現ドキュメント数: 6
単語: Purple, 出現ドキュメント数: 5
単語: Yellow, 出現ドキュメント数: 4
"""

出現回数で足切りするのではなく、残す単語数を指定したい場合は、keep_n を使えます。
(これにもデフォルト引数が入ってるので気をつけてください。元の単語数が100000を超えていたら、意図せず動作します)

5単語に絞り込むコードはこうなります。
単語は出現頻度が高い順に選ばれます。
no_below や no_aboveも同時に作用するので、これらの設定次第では、keep_nで指定したよりも少ない単語しか残らないことがあります。


# もう一度辞書の作成
dictionary = Dictionary(documents)
dictionary.filter_extremes(no_below=1, no_above=1, keep_n=5)
for word, id_ in dictionary.token2id.items():
    print(f"単語: {word}, 出現ドキュメント数: {dictionary.dfs[id_]}")
"""
単語: Blue, 出現ドキュメント数: 5
単語: Green, 出現ドキュメント数: 4
単語: Grey, 出現ドキュメント数: 4
単語: Magenta, 出現ドキュメント数: 6
単語: Purple, 出現ドキュメント数: 5
"""

あとは、あまり使わなさそうですが、 keep_tokens に単語を指定することで、no_belowや、no_aboveに関係なく、
その単語を残すことができます。


# もう一度辞書の作成
dictionary = Dictionary(documents)
dictionary.filter_extremes(no_below=5, no_above=1, keep_tokens=["White"])
for word, id_ in dictionary.token2id.items():
    print(f"単語: {word}, 出現ドキュメント数: {dictionary.dfs[id_]}")
"""
単語: Blue, 出現ドキュメント数: 5
単語: Magenta, 出現ドキュメント数: 6
単語: Purple, 出現ドキュメント数: 5
単語: White, 出現ドキュメント数: 1
"""

小ネタですが、Dictionaryオブジェクトは、各単語が出現したドキュメント数をdfs, 出現した回数をcfsという変数に保有しています。
filter_extremes を実行すると、dfsの方は単語が絞り込まれた上でidも振り直されるのですが、
cfsは単語が絞り込まれるだけで、idが振り直されません。
(なぜこんな仕様になっているのかは謎です。将来的に修正されるような気がします。)
直前のサンプルコードを動かした時点で、 dfsとcfs の中身を見たものがこちらです。
単語数が4個に減っているのは共通ですが、cfsの方はidが5とか10とか、元のままであることがわかります。


print(dictionary.dfs)
# {0: 5, 2: 5, 1: 6, 3: 1}
print(dictionary.cfs)
# {0: 5, 5: 5, 4: 6, 10: 1}

SciPyで多次元のカーネル密度推定

以前も紹介した、SciPyのカーネル密度推定のメソッド、gaussian_kdeの話です。
参考: SciPyによるカーネル密度推定

最近、多次元(と言っても2次元のデータですが)に対して、カーネル密度推定を行いたいことがあり、
どうせ1次元の場合と同じように使えるのだろうと適当に書いたら思うような動きにならず苦戦しました。

何をやろうとしたかというと、
[[x0, y0], [x1, y1], [x2, y2], …, [xn, yn]]
みたいなデータをそのままgaussian_kdeに渡してしまっていました。

ドキュメント: scipy.stats.gaussian_kde
をよく読むと、
Datapoints to estimate from. In case of univariate data this is a 1-D array, otherwise a 2-D array with shape (# of dims, # of data).
と書いてあります。

僕が渡そうとしたデータは shapeが[データ件数, 2(=次元数)]になっていたのですが、実際は
[2(=次元数), データ件数]の型で渡さないといけなかったわけです。

丁寧なことに、Examplesに取り上げられている例も1Dではなく2Dの例で、np.vstackとか使って書かれているので、以前の記事を書いた時にもっとしっかり読んでおけばよかったです。

使い方がわかったので、2次元のデータに対してやってみます。
サンプルデータは今回はscikit-learnのmake_moonsを使いました。


from sklearn.datasets import make_moons
import matplotlib.pyplot as plt

# データ生成。 今回はラベルは不要なのでデータだけ取得する。
data, _ = make_moons(
    n_samples=200,
    noise=0.1,
    random_state=0
)

fig = plt.figure(facecolor="w")
ax = fig.add_subplot(111, title="サンプルデータ")
ax.scatter(data[:, 0], data[:, 1])
plt.show()

サンプルデータはこんな感じです。

生成したデータを、gaussian_kdeにそのまま渡すとうまくいきません。
出来上がるモデルは、200次元のデータ2個を学習したものになっているからです。


from scipy.stats import gaussian_kde
# これはエラーは出ないが誤り
kde = gaussian_kde(data)

# 密度推定した結果を得ようとするとエラーになる。
print(kde.pdf([0, 0]))
# ValueError: points have dimension 1, dataset has dimension 200

きちんと、2次元のデータ200個を学習させるには、データを転置させて渡します。


from scipy.stats import gaussian_kde
# (データを、(次元数, データ件数) の型で渡す)
kde = gaussian_kde(X.T)

kde.evaluate を使って、推定した結果の値を得る時も、(次元数, データ件数)の形で、データを渡す必要があります。
(一点のみなら長さが次元数分の配列を渡せば良いです。)


# 点(1, 1)での値
print(kde.evaluate([1, 1]))
# [0.03921808]
# 4点(-0.5, 0), (0, 1), (0.5, 1), (1, 0) での値
print(kde.evaluate([[-0.5, 0, 0.5, 1], [0, 1, 1, 0]]))
# [0.07139624 0.2690079  0.2134083  0.16500181]

等高線を引いて図示する場合は次のように行います。(公式ドキュメントでは等高線ではなく、imshowで可視化していますね。)
こちらの記事も参考にしてください。
参考: matplotlibで等高線


# 等高線を引く領域のx座標とy座標のリストを用意する
x = np.linspace(-1.5, 2.5, 41)
y = np.linspace(-0.8, 1.3, 22)
# メッシュに変換
xx, yy = np.meshgrid(x, y)
# kdeが受け取れる形に整形
meshdata = np.vstack([xx.ravel(), yy.ravel()])
# 高さのデータ計算
z = kde.evaluate(meshdata)

# 可視化
fig = plt.figure(facecolor="w")
ax = fig.add_subplot(111, title="カーネル密度推定")
ax.scatter(data[:, 0], data[:, 1], c="b")
ax.contourf(xx, yy, z.reshape(len(y), len(x)), cmap="Blues", alpha=0.5)
plt.show()

結果がこちらです。

最後に、xもしくはy座標を固定して、断面をみる方法を紹介しておきます。
これはシンプルに、固定したい方は定数値でもう一方のデータと同じ長さの配列を作って、
固定しない方のデータを動かしてプロットするだけです。
x=-1.0, 0.5, 2.0 の3つの直線で切ってみた断面を可視化すると次のようなコードになります。


fig = plt.figure(facecolor="w")
ax = fig.add_subplot(111)
for x_ in [-1.0, 0.5, 2.0]:
    ax.plot(y, kde.pdf([[x_]*len(y), y]), label=f"x={x_}")
ax.set_ylabel("z")
ax.set_xlabel("y")
ax.legend()
plt.show()

結果がこちら。

だいぶ使い方の感覚が掴めてきました。

matplotlibで複数のグラフを含む図に全体のタイトルをつける

matplotlibで複数のグラフを含む図を作る時、全体にタイトルをつけたいことがあります。
個々のグラフはadd_subplotsする時にtitle引数で指定するか、もしくは、set_titleすれば指定できるのですが、
figureについては、plt.figure()する時にtitle引数を受け取ってくれないですし、
set_titleのようなメソッドも持っていません。

そのため、ちょっと迷っていたのですが、fig.suptitleでタイトルをつけることができることがわかりました。

ドキュメントはここです。: suptitle

適当な4グラフでやってみます。


import numpy as np
import matplotlib.pyplot as plt

fig = plt.figure(figsize=(10, 8), facecolor="w")
# 全体にタイトルをつける
fig.suptitle("全体のタイトル")

# 以下、適当なグラフを作る
ax = fig.add_subplot(2, 2, 1, title="棒グラフ")
ax.bar(range(5), np.random.randint(100, size=5))
ax = fig.add_subplot(2, 2, 2, title="折れ線グラフ")
ax.plot(range(5), np.random.randint(100, size=5))
ax = fig.add_subplot(2, 2, 3, title="散布図")
ax.scatter(np.random.randn(10), np.random.randn(10))
ax = fig.add_subplot(2, 2, 4, title="円グラフ")
ax.pie(np.random.randint(1, 10, size=5))

plt.show()

出力がこちら。

手軽ですね。

もしタイトルの位置を微調整したい場合は、
引数 x, y で位置を変えることができます。
デフォルトの値は x=0.5 と y=0.98 です。
それずれタイトルを左右中央に配置することと、図形の上部に置くことを意味しています。
xを0にすれば左揃えになりますし、 yを0付近にすれば図の下部にタイトルを置くことができます。

Numpyだけで重回帰分析

興味本位でNumPyの多項式回帰(polyfit)のソースコードを読んでいたのですが、
その中でNumPyにも重回帰分析のメソッドが用意されていて使われているのを見つけました。
てっきり重回帰分析は、scikit-learnかstatsmodelsを使うか、もしくはNumpyでやるならスクラッチ実装しないといけないと思い込んでいたので意外でした。
使ってみるとかなり手軽に使えたのでこの記事で紹介します。

ちなみに多項式回帰については既に記事を書いているのでご参照ください。
参考記事: NumPyで多項式回帰

紹介する関数はこちらです。
numpy.linalg.lstsq
(正直この名前はドキュメントでは探しにくい。おそらく、least-squaresの略語です。)

とりあえず、使ってみましょう。
ダミーデータとして、
$$
y = 3x_0 -2 x_1 + 5x_2 + \varepsilon
$$
のデータを作っておきます。$\varepsilon$はノイズです。


import numpy as np
X = np.random.randn(100, 3)
y = X@np.array([[3], [-2], [5]]) + np.random.randn(100).reshape(100, 1)

さて、早速lstsqを使ってみます。
使い方は簡単で、先ほどのXとyを渡してあげて、あと、rcondという引数を指定するだけです。
rcondは指定しないとwarningが出ますが、指定しなくても動きます。
小さい特異値を切り捨てる割合を指定する方法で、Noneか-1か指定しておけば良さそうです。


print(np.linalg.lstsq(X, y, rcond=None))
"""
(array([[ 2.97422787],
       [-2.01082975],
       [ 4.81883873]]),
       array([112.27261964]),
       3,
       array([10.54157065,  9.62381814,  8.55906245]))
"""

さて、ご覧の通り、結構いろいろな値がタプルで戻ってきました。
最初のArrayの
array([[ 2.97422787],
[-2.01082975],
[ 4.81883873]]),
の部分が推定した係数です。正解の 3, -2, 5 に近い値になっているのがわかります。
次の、[112.27261964]の値は残差の平方和です。
そして、 3 は Xのrank、次の array([10.54157065, 9.62381814, 8.55906245]) は Xの特異値です。

残渣平方和が 112.27261964 になるのは計算してみておきましょう。


coef, rss, rank, s = np.linalg.lstsq(X, y, rcond=None)
# 予測値を計算
p = X@coef
# 残差の平方和を計算
print(((y-p)**2).sum())
# 112.2726196447206

バッチリですね。

ここまでの流れで、気付いた方もいらっしゃると思いますが、lstsqで重回帰分析すると、定数項が出てきません。
定数項を含めて重回帰分析するには、Xに値が全部1になる列を追加して、それを渡す必要があります。
Numpyだけでもできますが、 statsmodels の add_constant あたりを使ってもいいでしょう。

例えば、
$$
y = 3x_0 -2 x_1 + 5x_2 + 4 + \varepsilon
$$
をダミーデータを作って、回帰分析するとこまでやると次のようになります。


import numpy as np
import statsmodels.api as sm

# ダミーデータ生成
X = np.random.randn(100, 3)
y = X@np.array([[3], [-2], [5]]) + 4 + np.random.randn(100).reshape(100, 1)

# Xに定数項を追加したデータを生成
X_add_const = sm.add_constant(X)

# 回帰分析
coef, rss, rank, s = np.linalg.lstsq(X_add_const, y, rcond=None)

# 推定した係数を表示
print(coef)
"""
[[ 4.01847308]
 [ 3.02763718]
 [-1.98363017]
 [ 4.99985177]]
"""

最初の 4.018…が定数項で残りが係数です。

matplotlibで複数のグラフを含むアニメーションを作成する

再びmatplotlibの動画の話です。
先日、一枚のfigureの中に複数のグラフを持っている図をアニメーションにしたいことがありました。
これもすぐにできると思ったのですが、方法が理解できるまで結構手間取ったので記録として残しておきます。

ドキュメントの ArtistAnimation のページには明記されてないし、
Examples の Animationの一覧 をみても、
サンプルは全部1枚の図に1つのグラフのもので、2グラフを動かしている例はないんですよね。(※記事執筆時点の情報です)

さて、方法ですが結論としては、次のようなコードで実現できました。
matplotlibでgif動画生成 の記事のコードがグラフ一つの例なので、見比べてみてください。


import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import ArtistAnimation
fig = plt.figure(facecolor="w")
ax1 = fig.add_subplot(2, 1, 1)
ax2 = fig.add_subplot(2, 1, 2)
# 0 <=x < 2pi の範囲の点列を作成。
x = np.linspace(0, 2*np.pi, 101)[: -1]
# 各コマの画像を格納する配列
image_list = []
# 1周だと短すぎたので5回繰り返す
for _ in range(5):
    for i in range(100):
        # ずらしながらsinカーブを描写し、配列に格納
        y1 = np.sin(np.roll(x, -i))
        y2 = np.cos(np.roll(x, -i))
        # 一つ目のグラフを描写する
        image1 = ax1.plot(x, y1, c="b")
        # 二つ目のグラフを描写する
        image2 = ax2.plot(x, y2, c="g")
        # 同時に描写したいグラフを連結したものを配列に追加する。
        image_list.append(image1+image2)
# アニメーションを作成
ani = ArtistAnimation(fig, image_list, interval=10)
# mp4ファイルに保存
ani.save('animation.mp4', writer='ffmpeg')
# gifファイルに保存する場合
# ani.save('animation.gif', writer='pillow')

ポイントになるのはこの部分です。
作成したグラフの配列を連結して、それをコマの配列に追加していっています。


image_list.append(image1+image2)

image1(とimage2)はそれぞれ、グラフの配列(要素は一個だけ)です。


print(image1)
# [<matplotlib.lines.Line2D object at 0x7fbe7dc4eed0>]
print(image2)
# [<matplotlib.lines.Line2D object at 0x7fbe7db8cf90>]

これを足すことで配列が連結されています。


print(image1+image2)
# [<matplotlib.lines.Line2D object at 0x7fbe7dc4eed0>, <matplotlib.lines.Line2D object at 0x7fbe7db8cf90>]

これ一コマとして、それらの配列を作り、ArtistAnimation メソッドに渡すことで動画にできます。
出来上がったのがこちらです。

最初、順番に更新したらいいかと思ってこういうのを試しましたが、これは2つのグラフが同時には表示されず、
1枚ずつ表示されるので非常にチカチカした見た目になります。


# 失敗例1
image_list.append(image1)
image_list.append(image2)

次に試みたのがこれ。これはこの後の ArtistAnimation でエラーになりました。
(配列の階層が深くなりすぎたようです)


# 失敗例2
image_list.append([image1, image2])

FFmpegでgif動画をmp4動画に変換する

以前、matplotlibでgif動画を作る方法を紹介しましたが、実際に業務で使っていると出来上がったgifファイルのサイズの大きさに困ることがよくありました。
参考: matplotlibでgif動画生成

gifファイルのままサイズを圧縮する方法をいろいろ調べていたのですが、どうやらmp4形式に変換するのが一番圧縮率が良さそうということがわかりました。
(mp4よりgifの方が軽いと勘違いしていたのですが、完全にただの思い込みだったようです。)

ということで、Macでgifファイルをmp4に変換する方法を紹介します。

利用するのは、FFmpegというツールです。
公式サイト: FFmpeg

HomeBrewでインストールできるという情報が各所にあったのですが、今時点の環境では失敗しました。


$ brew install ffmpeg 
#
# いろんなメッセージ (略)
#
Error: The following formula
  [#<Dependency: "python@3.9" []>, #<Options: []>]
cannot be installed as binary package and must be built from source.
Install the Command Line Tools:
  xcode-select --install

最後に出てくる、 xcode-select --install も試しても動かないのでお手上げです。

その一方で、 condaでもインストールできるらしいことがわかったので、僕はcondaでインストールしました。
(バージョンが古いという噂もあります。)


$ conda install ffmpeg

これでインストールしてしまえば、次のコマンドでmp4形式に変換できます。


$ ffmpeg -i [元のgifファイル名].gif -pix_fmt yuv420p [作成するmp4ファイル名].mp4 

先日の記事の sinカーブを動かす動画であれば、元々1.4MBもあったのが、33KBまで軽くなりました。