賢くなりたいトイプードルの日記

データサイエンス系の話をメインにしていきます

Rails6とSidekiqのアプリをDockerで作ってHerokuでリリース

docker heroku rails6 Ruby On Rails sidekiq

Ruby3.0、Rails6でSidekiqを使ったアプリをDocker上に作る際に引っかかったので、備忘録を書いておきます。

Herokuにリリースします。

参考にしたURL：https://blog.u6k.me/2018/04/29/sample-sidekiq.html

$ mkdir myapp
$ cd myapp

GemfileとGemfile.lock、Dockerfile、docker-compose.ymlを作成して、Gemfileを編集してgem 'sidekiq'とgem 'sidekiq-cron'を追加

$ touch Gemfile Gemfile.lock Dockerfile docker-compose.yml
$ cat Gemfile
source 'https://rubygems.org'
git_source(:github) { |repo| "https://github.com/#{repo}.git" }

ruby '3.0.1'
gem 'rails', '~> 6'

railsコマンドをDockerコンテナ化する。Dockerfileを作って以下のように編集

FROM ruby:3.0.1

RUN apt-get update && \
    apt-get install -y \
        nodejs postgresql-client && \
    curl -sS https://dl.yarnpkg.com/debian/pubkey.gpg | apt-key add - && \
    echo "deb https://dl.yarnpkg.com/debian/ stable main" | tee /etc/apt/sources.list.d/yarn.list && \
    apt-get update && apt-get install -y yarn && \
    apt-get clean

VOLUME /myapp
WORKDIR /myapp

COPY Gemfile .
COPY Gemfile.lock .

RUN bundle install

EXPOSE 3000

CMD ["rails", "server", "-b", "0.0.0.0"]

docker-compose.ymlを作って以下のように編集

version: '3'

services:
  web:
    build: .
    environment:
      - "RAILS_ENV=development"
      - "REDIS_URL=redis://redis:6379"
    volumes:
      - ".:/myapp"
    ports:
      - "3000:3000"
    depends_on:
      - db
      - redis
  worker:
    build: .
    environment:
      - "RAILS_ENV=development"
      - "REDIS_URL=redis://redis:6379"
    volumes:
      - ".:/myapp"
    depends_on:
      - "redis"
    command: sidekiq
  redis:
    image: redis:4.0
    volumes:
      - "redis:/data"
    command: redis-server --appendonly yes
  db:
    image: postgres
    ports:
      - '5432:5432'
    volumes:
      - postgresql-data:/var/lib/postgresql/data
    environment:
      - POSTGRES_PASSWORD=password

volumes:
  redis:
    driver: local
  postgresql-data:
    driver: local

Docker上でrails new

$ docker-compose run web rails new . --force --no-deps --database=postgresql --skip-bundle

Gemfileを編集してsidekiqを追加

# sidekiq
gem 'sidekiq'
gem 'sidekiq-cron'

イメージをビルドしてwebpackerをインストール

$ docker-compose build
$ docker-compose run web yarn install --check-files
$ docker-compose run web bundle exec rails webpacker:install

トップページをwelcomeページから変えて、config/routes.rbでrequireやmountをしておきます。

$ docker-compose exec web rails g controller home top
require 'sidekiq/web'
require 'sidekiq/cron/web'

Rails.application.routes.draw do
  mount Sidekiq::Web, at: '/sidekiq'
  get "/" => "home#top"
end

config/database.ymlを編集

default: &default
  adapter: postgresql
  encoding: unicode
  host: db
  username: postgres
  password: password
  pool: 5

development:
  <<: *default
  database: myapp_development

test:
  <<: *default
  database: myapp_test

production:
  <<: *default
  database: myapp_production
  username: myapp
  password: <%= ENV['MYAPP_DATABASE_PASSWORD'] %>

コンテナを起動

$ docker-compose up -d

dbを作ってマイグレート

$ docker-compose run web rake db:create
$ docker-compose run web rails db:migrate

これでlocalhost:3000を開けば表示されます。

Hello ワーカーを作成します

$ docker-compose exec web rails g sidekiq:worker Hello

app/worker/hello_worker.rbにputs "hello"を追加

class HelloWorker
  include Sidekiq::Worker

  def perform(*args)
    puts "hello"
  end
end

コンテナを再起動してrailsコンソールを開きます

$ docker-compose down -v
$ docker-compose up -d
$ docker-compose exec web rails c

ワーカーをキューに登録します。

> HelloWorker.perform_async

正常に登録された場合、文字列が表示されます。

キューに登録したワーカーはすぐに実行されたことをworkerコンテナのログを確認します

$ docker-compose logs worker
worker_1  | hello

という行があれば実行されています。

sidekiqのcron化をしていきます。

コンテナを再起動してRailsコンソールを開き、

$ docker-compose down -v
$ docker-compose up -d --build
$ docker-compose exec web rails c

sidekiqでクーロンを登録します

Sidekiq::Cron::Job.create name: "Hello Job", cron: "* * * * *", class: "HelloWorker"

ログを確認すると、約1分ごとにジョブが実行されていることがわかります。

docker-compose logs worker
orker_1  | 2021-06-14T21:07:33.930Z pid=1 tid=641 INFO: Booted Rails 6.1.3.2 application in development environment
worker_1  | 2021-06-14T21:07:33.931Z pid=1 tid=641 INFO: Running in ruby 3.0.1p64 (2021-04-05 revision 0fb782ee38) [x86_64-linux]
worker_1  | 2021-06-14T21:07:33.932Z pid=1 tid=641 INFO: See LICENSE and the LGPL-3.0 for licensing details.
worker_1  | 2021-06-14T21:07:33.932Z pid=1 tid=641 INFO: Upgrade to Sidekiq Pro for more features and support: https://sidekiq.org
worker_1  | 2021-06-14T21:07:33.933Z pid=1 tid=641 INFO: Booting Sidekiq 6.2.1 with redis options {:url=>"redis://redis:6379"}
worker_1  | 2021-06-14T21:10:10.749Z pid=1 tid=eox class=HelloWorker jid=4b94056a98e35c24d72cf61c INFO: start
worker_1  | hello
worker_1  | 2021-06-14T21:10:11.100Z pid=1 tid=eox class=HelloWorker jid=4b94056a98e35c24d72cf61c elapsed=0.348 INFO: done
worker_1  | 2021-06-14T21:11:03.202Z pid=1 tid=euh class=HelloWorker jid=1e975d07184b6d8b5cd52b45 INFO: start
worker_1  | hello
worker_1  | 2021-06-14T21:11:03.206Z pid=1 tid=euh class=HelloWorker jid=1e975d07184b6d8b5cd52b45 elapsed=0.003 INFO: done
worker_1  | 2021-06-14T21:12:08.061Z pid=1 tid=ev1 class=HelloWorker jid=188bc7a35b13bef271afb3c0 INFO: start
worker_1  | hello
worker_1  | 2021-06-14T21:12:08.065Z pid=1 tid=ev1 class=HelloWorker jid=188bc7a35b13bef271afb3c0 elapsed=0.004 INFO: done
worker_1  | 2021-06-14T21:13:22.170Z pid=1 tid=evl class=HelloWorker jid=6d5f82a01ae29804c9663749 INFO: start
worker_1  | hello
worker_1  | 2021-06-14T21:13:22.173Z pid=1 tid=evl class=HelloWorker jid=6d5f82a01ae29804c9663749 elapsed=0.003 INFO: done

Herokuでリリース

# ログイン
$ heroku login

# herokuのコンテナレジストリにログイン
$ heroku container:login

# アプリ作成
$ heroku create アプリ名

$ heroku config:set HOST=0.0.0.0

$ rm tmp/pids/server.pid      

# イメージを作成してコンテナレジストリにpush
$ heroku container:push web

# postgresqlアドオンの無料プランを追加
$ heroku addons:create heroku-postgresql:hobby-dev

# イメージをherokuへデプロイ
$ heroku container:release web

# DBセットアップ
$ heroku run rails db:migrate

# アクセス
$ heroku open

Pythonで微分方程式とラプラス変換

ラプラス変換 微分方程式 数学

よびのり先生のYouTube動画で、物理工学の教授二人と話しているのがあり、その中でこう語られていました。

・・・・なんかすごそう！

ということで、ラプラス変換について理解を深めていきたいと思います。

ラプラス変換は微分方程式を簡単にするとはどういうことなのか。まずは微分方程式のややこしさを理解していきます。

微分方程式

微分方程式とは、$\frac{dy}{dx} = y$を満たす関数を見つけるための方程式。この解は$ y = ce^x 　(c: 任意数) $です。

では問題を解いていきます。

例1

$$ y' = 2x+3, y(0) = 0(初期条件) $$

解:

$ y = x^2 + 3x + c　(一般解) $
$ y(0) = 0 $より
$ y = x^2 + 3x （特殊解）$

例2 :

$$ y^{\prime\prime} = e^x + 1, y(0) = 0, y'(0) = 1 $$

解：

$ y' = e^x + x + c_1 $
$ y = e^x + \frac{1}{2} x^2 + c_1 x + c_2 $
$ y(0) = 0 $より、$ 1 + c_2 = 0 $　$\therefore c_2 = -1$
$ y'(0) = 1 $より、$ 1 + c_1 = 1 $　$\therefore c_1 = 0$
以上より
$ y = e^x + \frac{1}{2} x^2 - 1 $

ここまでの微分方程式はかなり単純な微分方程式でした。

ただし、実際の現象をモデル化して微分方程式（＝数理モデル）を得て、それを解いて得た解を見てまた微分方程式を修正して、というのを繰り返して得られた微分方程式は、かなり複雑で解くのがかなり難しいものばかり。

まずは簡単な人工モデルを見ていきます。

例：人工モデル 「個体数の増加率はその時の個体数に比例する」（マルサスの法則）

つまり、

$$ y' = ay, y(0) = y_0 $$

これを解くと、

$ y = ce^{at} $
$ y(0) = y_0 $より、$ c = y_0 $
$ \therefore y = y_0 e^{at} $

この関数は指数関数的増加をするが、実際には資源に限りがあるので、微分方程式を以下のように修正する必要がある。

$$ y' = ay - by^2 $$

これを解くと、最初は指数関数的増加をしますが、あるところで一定になっていきます。

人口モデルはシンプルでしたが、たとえば人工衛星についての微分方程式はかなり複雑です。

こちらによれば、以下の常微分方程式（未知関数が本質的にただ一つの変数を持つ微分方程式）を解くと、木星の周りを回る人工衛星の位置と速度がわかります。

てすと

$$ \frac{d2x}{dt2} = - \frac{GM(x + tV_j)}{((x + tV_j)2 + y2)^{\frac{3}{2}}} $$

$$ \frac{d2y}{dt2} = - \frac{GMy}{((x + tV_j)2 + y2)^{\frac{3}{2}}} $$

$$ \frac{dx}{dt} = v_x $$

$$ \frac{dy}{dt} = v_y $$

$G$…万有引力定数、$ M $…木星の質量、$V_j$…木星の平均公転速度

ラプラス変換

ラプラス変換とは、座標変換の一種で、「微分方程式が簡単になる」という大きなメリットがあります。

微分方程式を解くのは大抵の場合難しいですが、微分方程式をラプラス変換して$s$世界に関する解を求めた後、それをラプラス逆変換して$t$世界に戻せば、微分方程式を解くことができるというわけです。

では、ラプラス変換とは何か。ラプラス変換は「関数は指数関数$e^{st}$のまとまりで表せる」という大胆な思想です。ちなみにフーリエ変換は「関数は三角関数のまとまりで表せる」という思想により、$t$（時間）空間から$\omega$（周波数）空間に移すもので、テイラー展開は「関数はその関数の色々な次数関数のまとまりで表せる」という思想でした。

微分方程式のひとつである運動方程式を考えてみます。

$$
\begin{equation}
\left\{ \,
\begin{aligned}
& v = \frac{dx}{dt} \\
& m\frac{dv}{dt} = F　(運動方程式) \\
& \frac{dv}{dt} = -\frac{k}{m}x　(バネの周期的な動きを表す微分方程式) \\
\end{aligned}
\right.
\end{equation}
$$

これより、

$$
\begin{equation}
\frac{d^2x}{dt^2} = -\frac{k}{m}x \\
\end{equation}
$$

となりますが、２回微分すると元に戻って$-\frac{k}{m}$という係数が前につく関数が$x$ということになります。２回微分すると元に戻る関数を考えると、たとえば

$$
\begin{equation}
\begin{aligned}
& sin(t) \\
& cos(t) \\
& Ae^{st}　(A=定数, s=\pm{i \sqrt{\frac{k}{m}}})
\end{aligned}
\end{equation}
$$

などがあります。なので、$x$は指数関数とか三角関数だろうと予想できますが、今のは計算ではなく、代入法です。頭が良い人ならわかるかもしれませんが、計算でわかるようなものではなさそうです。

そこで関数を指数関数のまとまりで表すことで、その関数の微分は指数関数$e^{st}$の$s$がおりてきた一次関数に、2階微分は指数関数$e^{st}$の$s$が２回おりてきて二次関数になり、一次関数や二次関数なら誰でも解くことができるようになります。つまり、「ラプラス変換により微分方程式を$s$の$n$次関数で表すことができる」ということです。

またオイラーの公式$ e^{i\theta} = cos\theta + i sin\theta $より、指数関数は複素数を使うと、三角関数を含むようになります。

関数を指数関数のまとまりで表すことで表すことができるということは、

$$ f(t) = \int F(s) e^{st} ds 　(原関数)$$

と表せますが、この$F(s)$を計算する方法のことをラプラス変換といいます。

$F(s)$だけを取り出すには$e^{st}$を消す必要があるので、ラプラス変換は以下のような数式で表せます。

$$ F(s) = \mathcal{L}[f(t)] = \int^{\infty}_0 f(t)e^{-st} dt 　(像関数)$$

$f(t)$が$F(s)$になってます。つまり、$t$の関数が$s$の関数に変換されています。

たとえば水力発電を例に取ってみると、川の流れの速さを$x$、水車の回る速度が$y$、モーターで発生する電圧が$z$、得られる電力が$w$とする時、普通$x$と$y$の間、$y$と$z$の間、$z$と$w$の間は微分方程式によって繋がれていますが、これが$x$から$w$を得る際に大変であり、ラプラス変換を使えば簡単に得られるということです。

$f(t) = 1$をラプラス変換すると、$F(s) = \frac{1}{s}$
$f(t) = t$をラプラス変換すると、$F(s) = \frac{1}{s^2}$
$f(t) = \frac{t^2}{2}$をラプラス変換すると、$F(s) = \frac{1}{s^3}$

となります。つまり、$f(t)$を$t$で微分すると、$ F(s)$を$s$倍することになっています。なので、

$f(t) = sin(t)$をラプラス変換すると、$F(s) = \frac{1}{s^2+1}$であれば
$f(t) = cos$のラプラス変換は、$F(s) = \frac{s}{s^2+1}$となります。

時間領域で解析しようと思っていたけれど微分方程式が難しいので、微分法則が簡単になる$F(s)$だけ見れば良い、ということができます。

$f(t) = x(t)$をラプラス変換すると、$F(s) = X(s)$とすると、
$f(t) = \dot{x}(t)$をラプラス変換すると、$F(s) = sX(s)$になるため、$x(t)$の部分積分を頑張らなくても、$x(t)$のラプラス変換を求めていればそれを$s$倍すれば良いだけ。つまり、

$$ \mathcal{L}[\dot{x}] = s\mathcal{L}[x] $$

が成り立ちます。

例題1

$$ \ddot{x} - 3\dot{x} + 2x = 3t^2 $$

この式をラプラス変換すると、

$ s^2X - 3sX + 2X = \frac{6}{s^3} $
$ (s^2 - 3s + 2)X = \frac{6}{s^3} $
$ X = \frac{1}{(s^2 - 3s + 2)} * \frac{6}{s^3} $

と求められます。そしてもし

$$ \ddot{x} - 3\dot{x} + 2x = f(t) $$

とするなら、そのラプラス変換による解は

$ X(s) = \frac{1}{(s^2 - 3s + 2)} * F(s) $

と表せます。$\frac{1}{(s^2 - 3s + 2)}$の部分は伝達関数といい、$G(s)$として

$$ X(s) = G(s)* F(s) $$

と表すこともあります。このように入出力の関係が見やすくなることがわかります。

そしてこれを時間領域$t$に逆変換すれば微分方程式を簡単に解くことができます。しかしラプラス変換のすごさはそういうことではなく、$ G(s) $がわかれば、システムの安定性がだいたいわかるということにあります。

つまり伝達関数$ G(s) $は制御システムに対する入力と出力の関係を表すものであり、入出力システムの挙動や安定性を評価することができる指標として考えることもできるということです。

このような入出力の関係を表す関数には入力信号の種類に関係なく、入力の設定部やセンサによる検出部、必要な制御入力を作る調節部などの制御要素に固有の関数となる事が要求されます。

ではシステムの安定性とは何かといえば、それは「入力なしの時に出力が収束するかどうか」と言えます。

Pythonでマルサスの法則を見てみる

マルサスの法則をPythonの記号計算ライブラリsympyを使って見てみます。

まずはライブラリのインポートと設定、変数や関数の定義をしていきます。

import sympy as sym
from sympy.plotting import plot
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()
sym.init_printing(use_unicode=True)
%matplotlib inline

x,y = sym.symbols("x y")
f = sym.Function('f')

マルサスの法則は$ y' = ay, y(0) = y_0 $でした。これを実装していきます。

a = 0.0078 # 日本の人工増加率

eq = sym.Eq(f(x).diff(x,1) - a*f(x))
eq

$ - 0.0078 f{\left(x \right)} + \frac{d}{d x} f{\left(x \right)} = 0 $

# 一般解
ans = sym.dsolve(eq)
ans

$f{\left(x \right)} = C_{1} e^{0.0078 x}$

# 特殊解
f_0 = 123731000 # 日本の人口
ans = sym.dsolve(eq, ics={f(0):f_0})
ans

可視化してみる

plot(ans.rhs, (x, 9000, 12000))

線形代数の練習問題をPythonでやってみた

python Python データサイエンス データサイエンスのための数学 掃き出し法 数学 線形代数 練習問題 逆行列

「データサイエンスのための数学」の練習問題を解いていて、手計算がとにかく苦手でミスしまくりなのでPythonで確認しながらやってみました。

この本には練習問題の解説がついてないので、独自の解説っぽいものも書いていきます。

※この記事は随時追記していきます。

練習問題 1.2 基本変形で逆行列を求める

問： 行列の基本変形によって以下の行列の逆行列をもとめよ

$$ \begin {pmatrix}
1 & 2 \\
2 & -1 \\
\end {pmatrix} $$

解法： $ \textbf{A}\textbf{A}^{-1} = \textbf{E} $より、$\textbf{A}$を$\textbf{E}$に変形する基本変形と同じ基本変形を$\textbf{E}$にすれば、それは$\textbf{A}^{-1}$になる。（$\textbf{E}$は単位行列）

pythonで掃き出し法をやっていきます。一番上の行から、順に対角成分を1にして同じ列を0にする、という処理を行うのが効率が良いとよびのり先生が言っておりました。なのでその方法でやっていきます。

a = np.array([[1,2],[2,-1]], dtype=float)
e = np.identity(2)
combine = [a,e]

for mat in combine:
  print(mat)
[[ 1.  2.]
 [ 2. -1.]]
[[1. 0.]
 [0. 1.]]

とりあえずnp.linalg.invでaの逆行列を見てみます

np.linalg.inv(a)
array([[ 0.2,  0.4],
       [ 0.4, -0.2]])

これが正解です。では掃き出し法をやっていきます。

# (2,1)成分を0にする
for mat in combine:
  mat[1] = mat[1] + mat[0]*(-2)
  print(mat)
[[ 1.  2.]
 [ 0. -5.]]
[[ 1.  0.]
 [-2.  1.]]
# (2,2)成分を1にする
for mat in combine:
  mat[1] = mat[1] * (-1/5)
  print(mat)
[[ 1.  2.]
 [-0.  1.]]
[[ 1.   0. ]
 [ 0.4 -0.2]]
# (1,2)成分を0にする
for mat in combine:
  mat[0] = mat[0] - mat[1]*2
  print(mat)
[[ 1.  0.]
 [-0.  1.]]
[[ 0.2  0.4]
 [ 0.4 -0.2]]

ということで逆行列を掃き出し法で求めることができました。

ちなみに、

5*e
array([[ 1.,  2.],
       [ 2., -1.]])

となるので、

$ \frac{1}{5}
\begin {pmatrix}
1 & 2 \\
2 & -1 \\
\end {pmatrix} $

とスカラーとの積で綺麗にできます。

練習問題1.3 逆行列で連立方程式を解く

問： 逆行列を求めることで以下の連立方程式を解け

$$ (1) \begin{equation}
\left\{ \,
\begin{aligned}
& 2x + 3y = 350 \\
& 5x + 6y = 800
\end{aligned}
\right.
\end{equation}$$

解法： この連立方程式を$ \textbf{A}x = b $の形にすると、

$$ \begin {pmatrix}
2 & 3 \\
5 & 6 \\
\end {pmatrix}
\begin {pmatrix}
x \\
y \\
\end {pmatrix}
=
\begin {pmatrix}
350 \\
800 \\
\end {pmatrix} $$

また、

$ \textbf{A}^{-1} \textbf{A} \textbf{x} = \textbf{A}^{-1} \textbf{b} $
$ \textbf{E} \textbf{x} = \textbf{A}^{-1} \textbf{b} $
$ \textbf{x} = \textbf{A}^{-1} \textbf{b} $

なので、$ \textbf{A} $の逆行列と$ \textbf{b} $の積が$ \textbf{x} $になることがわかる。

ではPythonで逆行列を求めていく。

a = np.array([[2,3],[5,6]], dtype=float)
e = np.identity(2)
combine = [a,e]

print(np.linalg.inv(a))
# out:
# array([[-2.        ,  1.        ],
#        [ 1.66666667, -0.66666667]])

for mat in combine:
  print(mat)
# out:
# [[2. 3.]
#  [5. 6.]]
# [[1. 0.]
#  [0. 1.]]

for mat in combine:
  mat[0] = mat[0] * (1/2)  # (1,1)成分を1にする
  mat[1] = mat[1] - mat[0]*5  # (2,1)成分を0にする
  mat[1] = mat[1]*(-2/3)  # (2,2)成分を1にする
  mat[0] = mat[0] + mat[1]*(-1.5)  # (1,2)成分を0にする
  print(mat)
# out:
# [[ 1.  0.]
#  [-0.  1.]]
# [[-2.          1.        ]
#  [ 1.66666667 -0.66666667]]

これで逆行列を求めることができた。

$ \textbf{x} = \textbf{A}^{-1} \textbf{b} $より、

x = e.dot(np.array([350, 800]))
print(x)
# out:
# array([100.,  50.])

ということで、$x = 100, y = 50$が解である。

では次の問題

$$ (2) \begin{equation}
\left\{ \,
\begin{aligned}
& x + 2y - z = 20 \\
& 2x - y + 2z = 110 \\
& -x - y + z = -10
\end{aligned}
\right.
\end{equation}$$

a = np.array([[1,2,-1],[2,-1,2],[-1,-1,1]], dtype=float)
e = np.identity(3)
combine = [a,e]

for mat in combine:
  mat[1] = mat[1] + mat[2]*2
  mat[2] = mat[2] + mat[0]

  mat[1] = mat[1] * (-1/3)
  mat[0] = mat[0] + mat[1]*(-2)
  mat[2] = mat[2] + mat[1]*(-1)

  mat[2] = mat[2]*(3/4)
  mat[0] = mat[0] + mat[2]*(-5/3)
  mat[1] = mat[1] + mat[2]*(4/3)

  print(mat)
# out: 
# [[ 1.00000000e+00  0.00000000e+00 -2.22044605e-16]
#  [ 0.00000000e+00  1.00000000e+00  0.00000000e+00]
#  [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
# [[-0.25  0.25 -0.75]
#  [ 1.    0.    1.  ]
#  [ 0.75  0.25  1.25]]

print(4 * e)
# out: 
# array([[-1.,  1., -3.],
#        [ 4.,  0.,  4.],
#        [ 3.,  1.,  5.]])

print(e.dot(np.array([20,110,-10])))
# out:
# array([30., 10., 30.])

よって逆行列は

$ \frac{1}{4}
\begin {pmatrix}
-1 & 1 & -3 \\
4 & 0 & 4 \\
3 & 1 & 5
\end {pmatrix} $

であり、この連立方程式の解は

$ x=30, y=10,z=30 $となる

練習問題1.11 左右から行列を掛けて標準化

問： $p$個の観測項目についての$n$人のデータである$n\times p$行列$\textbf{X}$に対し、$\textbf{A}\textbf{X}\textbf{B}$のように左右から掛けることで行列$\textbf{X}$を標準化する行列$\textbf{A}, \textbf{B}$を求めよ。ただし、各観測項目の平均値と分散はそれぞれ$ (\bar{x_1}, \sigma_1^2), ... ,(\bar{x_n}, \sigma_n^2) $とする。

解法： 標準化は$ \frac{x_{ij} - \mu_j}{\sigma_j} $で表されるように、それぞれのデータから各項目ごとの平均値を引いて、それを各項目ごとの標準偏差で割れば良い。

よって、$\textbf{A},\textbf{B}$を以下のような行列にする。

解：$\textbf{A} = \textbf{E}_n - \frac{1}{n} \textbf{1}_n $, $\textbf{B} = diag(\frac{1}{\sigma_i})$

この$\textbf{A}$のような行列を、中心化作用素行列といいます。$\frac{1}{n} \textbf{1}_n$は全ての要素が$\frac{1}{n}$である$n$次元ベクトルで、もしこれだけを行列に左から掛けると以下のように、それぞれの列の平均値が要素であるm次元ベクトルになることがわかります。

$$ \begin{eqnarray}
\frac{1}{n} \textbf{1}_n^T \textbf{X} &=& \frac{1}{n} (1, ... 1)
\begin{pmatrix}
x_{11} & \cdots & x_{1m} \\
\vdots & & \vdots \\
x_{n1} & \cdots & x_{nm} \\
\end{pmatrix} \\
&=& ( \sum_{i=1}^n \frac{1}{n} x_{i1}, ... ,\sum_{i=1}^n \frac{1}{n} x_{im} ) \\
&=& (\frac{1}{n} \sum_{i=1}^n x_{i1}, ... ,\frac{1}{n} \sum_{i=1}^n x_{im} )
\end{eqnarray} $$

$\textbf{A} \textbf{X}$を考えると、これは

$ \begin{eqnarray}
\textbf{A} \textbf{X} &=& \textbf{E}_n \textbf{X} - \frac{1}{n} \textbf{1}_n \textbf{X} \\
&=& \textbf{X} - \frac{1}{n} \textbf{1}_n \textbf{X}
\end{eqnarray}
$

となり、これは行列$\textbf{X}$の全ての要素から、その要素の列全体の平均値を引いた行列になります。

また$\textbf{B} = diag(\frac{1}{\sigma_i})$について考えると、これは、

$ \begin{pmatrix}
\frac{1}{\sigma_{1}} & \cdots & 0 \\
\vdots & \ddots & \vdots \\
0 & \cdots & \frac{1}{\sigma_{p}} \\
\end{pmatrix} $

ですが、これを右から掛けることにより、行列$\textbf{X}$の全ての要素を、要素の列全体の標準偏差で割るということになります。

Pythonで実際にやってみます。データはおなじみirisを使っていきます。

import seaborn as sns
import pandas as pd
iris = sns.load_dataset('iris')

iris.describe()

sepal_lengthsepal_widthpetal_lengthpetal_width
count150.000000150.000000150.000000150.000000
mean5.8433333.0573333.7580001.199333
std0.8280660.4358661.7652980.762238
min4.3000002.0000001.0000000.100000
25%5.1000002.8000001.6000000.300000
50%5.8000003.0000004.3500001.300000
75%6.4000003.3000005.1000001.800000
max7.9000004.4000006.9000002.500000

これを標準化するのは以下のようにscikit-learnを使えば簡単にできます。

from sklearn.preprocessing import StandardScaler
sc = StandardScaler()
data_std = sc.fit_transform(iris.iloc[:,:-1])

pd_data_std = pd.DataFrame(data_std)
pd_data_std.describe()

0123
count1.500000e+021.500000e+021.500000e+021.500000e+02
mean-4.736952e-16-7.815970e-16-4.263256e-16-4.736952e-16
std1.003350e+001.003350e+001.003350e+001.003350e+00
min-1.870024e+00-2.433947e+00-1.567576e+00-1.447076e+00
25%-9.006812e-01-5.923730e-01-1.226552e+00-1.183812e+00
50%-5.250608e-02-1.319795e-013.364776e-011.325097e-01
75%6.745011e-015.586108e-017.627583e-017.906707e-01
max2.492019e+003.090775e+001.785832e+001.712096e+00


ではこれを、左右から行列を掛けることで再現していきます。

X = iris.drop(columns="species")

n = 150
p = 4
A = np.identity(n) - (1/n)*np.ones(n)

pd.DataFrame(A.dot(X.values)).describe()

0123
count1.500000e+021.500000e+02150.0000001.500000e+02
mean4.263256e-16-1.136868e-150.000000-4.263256e-16
std8.280661e-014.358663e-011.7652987.622377e-01
min-1.543333e+00-1.057333e+00-2.758000-1.099333e+00
25%-7.433333e-01-2.573333e-01-2.158000-8.993333e-01
50%-4.333333e-02-5.733333e-020.5920001.006667e-01
75%5.566667e-012.426667e-011.3420006.006667e-01
max2.056667e+001.342667e+003.1420001.300667e+00

$\textbf{A}\textbf{X}$で、平均値が0になりました。標準偏差は変化していません。

sigmas = [1/(X[c].std()) for c in X.columns]
B = np.diag(sigmas)

df_axb = pd.DataFrame(A.dot(X.values).dot(B))
df_axb

0123
0-0.8976741.015602-1.335752-1.311052
1-1.139200-0.131539-1.335752-1.311052
2-1.3807270.327318-1.392399-1.311052
3-1.5014900.097889-1.279104-1.311052
4-1.0184371.245030-1.335752-1.311052
...............
1451.034539-0.1315390.8168591.443994
1460.551486-1.2786800.7035640.919223
1470.793012-0.1315390.8168591.050416
1480.4307220.7861740.9301541.443994
1490.068433-0.1315390.7602110.788031

統計量を調べてみます。

df_axb.describe()

0123
count1.500000e+021.500000e+021.500000e+021.500000e+02
mean5.684342e-16-2.676378e-15-9.473903e-17-5.684342e-16
std1.000000e+001.000000e+001.000000e+001.000000e+00
min-1.863780e+00-2.425820e+00-1.562342e+00-1.442245e+00
25%-8.976739e-01-5.903951e-01-1.222456e+00-1.179859e+00
50%-5.233076e-02-1.315388e-013.353541e-011.320673e-01
75%6.722490e-015.567457e-017.602115e-017.880307e-01
max2.483699e+003.080455e+001.779869e+001.706379e+00

平均値が0、標準偏差が1になっているので、標準化できています。