NumPyを使って線形モデルのパラメータを最小二乗法で推定する

統計学 Python

線形モデルと最小二乗法について調べることがありましたので、まとめておきます。

今回は青本の2章を参考に、繰り返し測定、1元配置、多項式の3つのモデルで例示します。また、NumPyを使った行列計算でナイーブに最小二乗法を実装します。

自然科学の統計学 (基礎統計学)

自然科学の統計学 (基礎統計学)

線形モデル

ある観測値について未知母数の線形式＋誤差を仮定するモデルを、線形モデルと言います。誤差は正規分布 N(0,\sigma^2) を仮定します。

繰り返し測定

飲料中のある成分の含有量を繰り返し測定した場合に以下の観測値が得られたとします。

No. 1 2 3 4 5 6 7 8 9
含有量(%) 37.55 37.74 37.69 37.71 37.73 37.81 37.53 37.58 37.74

こうした繰り返し実験では、以下のようにモデル化できます。

$$ y_i=\mu+\epsilon_i (i=1, \cdots , n) $$

1元配置

例えば、4種類の脂肪で揚げたドーナツの脂肪吸収量(g)を、各6回計測し、以下の値が得られたとします。

脂肪1 脂肪2 脂肪3 脂肪4
64 78 75 55
72 91 93 66
68 97 78 49
77 82 71 64
56 85 63 70
95 77 76 68

このようにa通りの条件で各r回の独立した実験を行った場合などは、1元配置のモデルとして以下のように表現できます。なお 1, \cdots , i , \cdots , a 、 1, \cdots , j , \cdots , r です。

$$ y_{ij}=\mu+\alpha_i+\epsilon_{ij} $$

多項式

液体中のある成分の含有量によって曇点の変化を観測した実験を想定します。散布図にプロットしてみると、曲線を描いていることが確認できます。

含有量(%) 0 1 2 3 4 5 6 7 8 0 2 4 6 8 10 0 3 6 9
曇点(℃) 22.1 24.5 26.0 26.8 28.2 28.9 30.0 30.4 31.4 21.9 26.1 28.5 30.3 31.5 33.1 22.8 27.3 29.8 31.8

f:id:ohke:20180616095332p:plain

2次多項式を仮定すると、以下の式でモデル化できます。

$$ y_i=\beta_0+\beta_1x_i+\beta_2x_i^2+\epsilon_i $$

行列による一般化

観測値 \vec{y}=(y_1, \cdots , y_n)^T 、未知母数 \vec{\theta}=(\theta_1, \cdots , \theta_p)^T 、誤差 \vec{\epsilon}=(\epsilon_1, \cdots , \epsilon_n)^T 、係数行列 {\bf X} (n行p列のデザイン行列)を用いると、線形モデルは以下で一般化して表現できます。

$$ \vec{y} = {\bf X}\vec{\theta}+\vec{\epsilon} \Leftrightarrow \begin{pmatrix} y_{1} \\ \vdots \\ y_{n} \end{pmatrix}= \begin{pmatrix} x_{1,1} & \ldots & x_{1,p} \\ \vdots & \cdots & \vdots \\ x_{n,1} & \ldots & x_{n,p} \end{pmatrix} \begin{pmatrix} \theta_{1} \\ \vdots \\ \theta_{p} \end{pmatrix} + \begin{pmatrix} \epsilon_{1} \\ \vdots \\ \epsilon_{n} \end{pmatrix} $$

最小二乗法

観測値 \vec{y} とその期待値の偏差二乗和(以下式)を最小にする解 \hat{\vec{\theta}} を求めます。

$$ S(\vec{\theta})=||\vec{y} - {\bf X} \vec{\theta}||^2=(\vec{y} - {\bf X} \vec{\theta})^{T}(\vec{y} - {\bf X} \vec{\theta}) $$

  • 繰り返し測定のモデルでは、 S(\vec{\theta})=\sum(y_i-\mu)^2
  • 1元配置のモデルでは、 S(\vec{\theta})=\sum\sum(y_{i,j}-\mu-\alpha_i)^2
  • 多項式のモデルでは、 S(\vec{\theta})=\sum(y_i-\beta_0-\beta_1x_i-\beta_2x_i^2)^2

\hat{\vec{\theta}} は最良線形不偏推定量(BLUE)となるのが、最小二乗法の原理です。

正規方程式

最小二乗法の解 \hat{\vec{\theta}} を求める方程式は以下となります。これを正規方程式といいます。

$$ {\bf X}^T{\bf X}\hat{\vec{\theta}}={\bf X}^T\vec{y} \Leftrightarrow \hat{\vec{\theta}}=({\bf X}^T{\bf X})^{-1}{\bf X}^T\vec{y} $$

ただし、上の通り {\bf X } のランクが未知母数の次元数よりも小さくなる場合は、不定解となります。後ほど見ていきますが、1元配置の例がこれに当たります。

Pythonによる実装

NumPyによる実装を示します。

import numpy as np

繰り返し測定

  • デザイン行列 {\bf X} は1行6列で、要素は全て1となります
  • \hat{\mu}=37.68 となりますが、これは \bar{y} と一致します
# 9個のデータ
y = np.array([37.55, 37.74, 37.69, 37.71, 37.73, 37.81, 37.53, 37.58, 37.74])

# 1行9列のデザイン行列
X = np.array([np.ones(y.shape)]).T

theta_hat = np.linalg.inv(X.T @ X) @ X.T @ y
# array([37.67555556])

1元配置

1元配置の場合は上の式をそのまま使えません。というのも、全ての方程式において \mu+\alpha_i の形で現れており、どう操作しても \mu と \alpha_i を分離できず、逆行列 {\bf X}^T を計算できないのです(ランク落ち)。

そこで純粋な逆行列の代わりに、一般逆行列 {\bf X}^{-} を採用した \hat{\vec{\theta}}=({\bf X}^T{\bf X})^{-}{\bf X}^T\vec{y} を解きます。
一般逆行列については、こちら↓の投稿が参考になります。

https://www.slideshare.net/wosugi/ss-79624897

zellij.hatenablog.com

  • 一般逆行列はnumpy.linalg.pinvで計算できます
  • デザイン行列 {\bf X} は24行(r*a)5列(a+1)
    • \mu に対応する1列目は全て1
    • \alpha_1 に対応する2列目は1〜6行目が1、 \alpha_2 に対応する3列目は7〜12行目が1, .... として、それ以外は0
  • 結果は \hat{\vec{\theta}}=(\mu, \alpha_1, \alpha_2, \alpha_3, \alpha_4)=(59, 13, 26, 17, 3) 
# 1元配置
y = np.array([
    64, 72, 68, 77, 56, 95,
    78, 91, 97, 82, 85, 77,
    75, 93, 78, 71, 63, 76,
    55, 66, 49, 64, 70, 68
]).T

# 24行5列
X = np.array([
    np.ones(24),
    np.concatenate([np.ones(6), np.zeros(18)]),
    np.concatenate([np.zeros(6), np.ones(6), np.zeros(12)]),
    np.concatenate([np.zeros(12), np.ones(6), np.zeros(6)]),
    np.concatenate([np.zeros(18), np.ones(6)])
]).T

# 行列式 np.linalg.det(X.T @ X) の値は0なので逆行列が存在しない
#theta_hat = np.linalg.inv(X.T @ X) @ X.T @ y

# 逆行列の代わりに一般逆行列を使って推定する
theta_hat = np.linalg.pinv(X.T @ X) @ X.T @ y
# array([59., 13., 26., 17.,  3.])

多項式

  • デザイン行列 {\bf X} は19行3列
    • \beta_0 に対応する1列目は全て1
    • \beta_1 に対応する2列目は x , \beta_2 に対応する3列目は x^2 とする
  • 結果は \hat{\vec{\theta}}=(\beta_0, \beta_1, \beta_2)=(22.561, 1.668, -0.068) 
# 19個のデータ
x = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 0, 2, 4, 6, 8, 10, 0, 3, 6, 9])
y = np.array([22.1, 24.5, 26.0, 26.8, 28.2, 28.9, 30.0, 30.4, 31.4, 21.9, 26.1, 28.5, 30.3, 31.5, 33.1, 22.8, 27.3, 29.8, 31.8])

# 19行3列のデザイン行列
X = np.array([np.ones(y.shape), x, x**2]).T

theta_hat = np.linalg.inv(X.T @ X) @ X.T @ y
# array([22.56123063,  1.66802044, -0.06795836])

まとめ

今回は線形モデルとそのパラメータを求める最小二乗法について、3つの例で見ていきました。またNumPyで正規方程式を解く実装を行いました。

PythonでRedisを参照・更新する

Python Redis

仕事でPythonアプリケーションからアクセスするRedisの導入を検討した際に、redis-pyでRedisを参照・更新する方法について調べましたので、備忘録にしておきます。

redis-pyのドキュメントはこちらです。
http://redis-py.readthedocs.io/

DockerでRedisコンテナを起動

実験用の環境をDockerで立ち上げ、localhost:6379で公開します。

$ docker run --name redis -d -p 6379:6379 redis redis-server --appendonly yes
$ docker ps
CONTAINER ID        IMAGE                          COMMAND                  CREATED              STATUS              PORTS                    NAMES
0481cc795f8f        redis                          "docker-entrypoint.s…"   About a minute ago   Up About a minute   0.0.0.0:6379->6379/tcp   redis

インストール

いつもどおりpipでインストールします。

$ pip install redis
$ pip list | grep redis
redis              2.10.6

接続

StrictRedisクラスを生成することでRedisとのコネクションを張ります。
なお、Redisクラスというのもあり、これでも接続できますが、後方互換性維持のために残されているインタフェースなので、新しく作るアプリケーションではStrictRedisを使ったほうが良さそうです。

import redis

# ホスト、ポート、DB番号を指定
conn = redis.StrictRedis(host='localhost', port=6379, db=0)

ConnectionPoolクラスでコネクションプールも作れます。
こちら↓の方の投稿では、接続まで倍くらい早いようです。

PythonでRedisを扱う(redis-pyの基本) - [Dd]enzow(ill)? with DB and Python

redis_pool = redis.ConnectionPool(host='localhost', port=6379, db=0, max_connections=4)
conn = redis.StrictRedis(connection_pool=redis_pool)

文字列値のset/get

まずはkey-valueのvalueが単一の文字列の場合を見ていきます。

基本は、StrictRedisクラスのset(key, value)で値の更新、get(key)で値の参照となります。

# 値のset
conn.set('key01', 'value01')  # 成功するとTrueが返る

# 値のget
value = conn.get('key01')  # バイナリ値で取得

print(type(value))
# <class 'bytes'>

print(str(value, encoding='utf-8'))  # str()で変換
# 'value01'

存在しないキーでgetするとNoneが返されます。

print(conn.get('keyXX'))
# None

また、存在しているキーでsetすると上書きされます。

conn.set('key01', 'new value')
conn.get('key01')
# b'new value'

なお、キーの削除はdelete()

conn.delete('key01')  # 返り値は削除したキー数
conn.delete('key02', 'key03')  # 複数同時削除も可能

有効期限

キーに有効期限を設定することもできます。
引数exに秒数を指定すると、その秒数経過後にgetしても値は取り出せなくなります。ミリ秒で指定する場合は引数pxを使います。

conn.set('key02', 'value02', ex=10)

print(conn.get('key02'))
# b'value02'

# 10秒後に再度get
print(conn.get('key02'))
# None

set後に有効期限を更新する場合、expire(key, seconds)で、有効期限を秒数で設定することができます。また、expireat(key, when)ならdatetime型の値を渡すことで"いつ無効にするのか"といったこともできます。

import datetime

# "key03"を10秒後に削除
conn.expire('key03', 10)

# "key04"を2018/6/3 17時に削除
conn.set('key04', 'value04')
conn.set('key04', datetime.datetime(2018, 6, 3, 17))

リスト

valueに使えるのは単一の値だけではありません。複数の値を関連付けられるvalueの型として、リスト、ハッシュ、セットというのが提供されています。

リスト型は、順序を持つ連結リストです。左端(先頭)または右端(末尾)から値を挿入・削除できます。

  • lpush / rpush : で値を追加
  • lrange : インデックスで値を参照 (参照された値は削除されない)
  • lpop / rpop : 単一の値を参照 (参照された値は削除される)
# 末尾へ"A", "B", "C"の順で追加
conn.rpush('key01', 'A')  # 1
conn.rpush('key01', 'B')  # 2
conn.rpush('key01', 'C')  # 3

# 先頭へ"Z"を追加
conn.lpush('key01', 'Z')  # 4

# 先頭から2つを参照
conn.lrange('key01', 0, 1)  # [b'Z', b'A']
# 先頭から4つを参照
conn.lrange('key01', 0, 3) # [b'Z', b'A', b'B', b'C']

# 先頭から1つを取り出す (取り出された"Z"は削除されます)
conn.lpop('key01')  # b'Z'
conn.lrange('key01', 0, 2)  # [b'A', b'B', b'C']

# 末尾から1つを取り出す (取り出された"C"は削除されます)
conn.rpop('key01')  # b'C'
conn.lrange('key01', 0, 1)  # [b'A', b'B']

ハッシュ

値が(Pythonで言うところの)ディクショナリ型となっている、ハッシュ型という型もあります。

幅(width)と高さ(height)を持つ長方形を格納する場合を考えてみましょう。

  • hset : ハッシュ型の追加
  • hget : ハッシュ型の値の取得
  • hgetall : キーに紐づく全ての値をdict型で
  • hmset : ハッシュ型の値をdict型で一括追加
# ハッシュ型の追加
conn.hset('rect1', 'width', '10.0')
conn.hset('rect1', 'height', '15.0')

# ハッシュ型の参照
conn.hget('rect1', 'width')  # b'10.0'
# キーに紐づくすべての値をdict型で返す
conn.hgetall('rect1')  # {b'width': b'10.0', b'height': b'15.0'}

# dict型で追加
conn.hmset('rect2', {'width':'7.5', 'height':'12.5'})
conn.hgetall('rect2')  # {b'width': b'7.5', b'height': b'12.5'}

セット

集合内に重複した値を持たないのがセット型です。Pythonのset型の値と対応します。

  • sadd : セット型の値の追加
  • smembers : セット型の値の参照で、set型で返される
  • sinterやsunionなどの集合演算が使える
# セット型の追加
conn.sadd('key01', 'A')

# セット型の参照
conn.smembers('key01')  # {b'A'}

conn.sadd('key01', 'B')
conn.sadd('key01', 'C')
conn.sadd('key01', 'A')  # 重複キーでsaddした場合、返り値は0

# 重複する"A"は1つになる
conn.smembers('key01')  # {b'A', b'C', b'B'}

conn.sadd('key02', 'B', 'D', 'B', 'A')
conn.smembers('key02')  # {b'D', b'A', b'B'}

# 積集合
conn.sinter('key01', 'key02')  # {b'A', b'B'}

# 和集合
conn.sunion('key01', 'key02')  # {b'D', b'A', b'C', b'B'}

ソート済みセット

上述のセットは順序を持ちませんが、ソート済みセットはある浮動小数(スコアと呼ばれる)によってソートされたセットです。

例えば、4人の生徒の数学のテストの点数をソート済みセットで表現してみます。

  • zadd : ソート済みセットの追加
  • zrange : ソート済みセットの取得
    • 引数withscoresでスコア付きのタプルリストが返る
    • 引数descをTrueにすると降順ソート(デフォルトは昇順)
# ソート済みセットの追加
conn.zadd('math', 75.0, 'Tanaka')

# ソート済みセットの取得
conn.zrange('math', 0, 1)  # [b'Tanaka']
# スコア付き
conn.zrange('math', 0, 0, withscores=True)  # [(b'Tanaka', 75.0)]

conn.zadd('math', 95.0, 'Suzuki', 60.0, 'Sato', 75.0, 'Komeda')
conn.zcard('math')  # 4

# 昇順
conn.zrange('math', 0, 3, withscores=True)  # [(b'Sato', 60.0), (b'Komeda', 75.0), (b'Tanaka', 75.0), (b'Suzuki', 95.0)]

# 降順
conn.zrange('math', 0, 3, desc=True, withscores=True)  # [(b'Suzuki', 95.0), (b'Tanaka', 75.0), (b'Komeda', 75.0), (b'Sato', 60.0)]

# 降順でトップ3つ
conn.zrange('math', 0, 2, desc=True, withscores=True)  # [(b'Suzuki', 95.0), (b'Tanaka', 75.0), (b'Komeda', 75.0)]

おわりに

redis-pyを使って、PythonからRedisをCRUDする方法についてまとめました。
ここで紹介した以外にも、ビット配列型やインクリメンタルに参照するスキャンなどもredis-pyで扱えますので、ドキュメントを参考にしてみてください。

Google AnalyticsのデータをBigQueryで集計・分析するときのテクニック集

BigQuery Google Analytics

先週の投稿で、Google AnalyticsのサンプルデータをBigQueryでクエリできるようにしました。

BigQueryを有効化してGoogle Analyticsのサンプルデータにクエリできるようにする - け日記

今回はBigQueryを使ってGoogle Analytics (GA)のデータを集計・分析するとに役立ったテクニックを紹介します。
ここではbigquery-public-dataデータセットのgoogle_analytics_sampleテーブルを使います。

スキーマ定義

GAのデータをBigQueryにインポートすると、以下のスキーマでテーブルが生成されます。項目の意味を確認したい時に、都度都度参照します。

support.google.com

かなりの項目数ですので、適宜タブ補完も活用しましょう。

Standard SQLを使う

BigQueryはLegacy SQLとStandard SQLの2つが利用でき、デフォルトではLegacy SQLです。

Standard SQLはSQL 2011と互換性があり、この後説明するWITH句やSELECT内のサブクエリなど、Standard SQLでしかできないこともあります。そのため、これから作るクエリについてはStandard SQLで書くことをおすすめします。以降のテクニックもStandard SQLを前提としています。
2つの違いについては、こちらの投稿でまとめられています。

tech.starttoday-tech.com

画面から Show Options > SQL Dialect で"Use Legacy SQL"のチェックを外すことで無効化されます。

また、SQLの1行目で#standardSQLと記述すると、Standard SQLで解釈・実行されます。SQLをコピーして共有したりすることもあるかと思いますので、おすすめはこちらの方法です。

#standardSQL

SELECT
  SUM(totals.transactions) AS transactions,
  SUM(totals.transactionRevenue) / 1000000 AS revenue
FROM `bigquery-public-data.google_analytics_sample.ga_sessions_20170801`

パーティションを指定する

GAのデータは日付ごとにパーティションされており、ga_sessions_YYYYMMDDという名前になっています。
そのため、最初に日付を範囲指定して1テーブルにして取り出す必要があります(実質的にはUNION)。

例えば、2017/7/3〜2017/7/30の4週間分を取り出すには、以下のようにクエリします。FROM句のテーブル末尾*でワイルドカード指定できるようにし、_TABLE_SUFFIX BETWEEN 'YYYYMMDD' AND 'YYYYMMDD'とすることで指定の日付でフィルタリングしています。

#standardSQL

# 2017/7/3〜2017/7/30のGAのデータ
SELECT *
FROM `bigquery-public-data.google_analytics_sample.ga_sessions_*`
WHERE _TABLE_SUFFIX BETWEEN '20170703' AND '20170730'

WITH句でクエリをまとめる

Standard SQLの大きな利点の1つが、WITH句を使ってクエリを再利用できることです。
Legacy SQLでは、FROM句やJOIN句で別々に記述する、あるいは、ビューを作る、といったことをしなければならず、非常に読みづらいクエリになりがちでした。
GAのデータはデバイス毎やランディングページ毎、リファラ毎で別々に集計することがよくあるため、WITH句を使って1箇所で定義してそれに名前をつけてあげることで、クエリの見通しがとても良くなります(よく使う定型のクエリについては、ビューを検討した方が良いかもしれません)。

例えば、上でパーティション指定して取得した行データから、PCとスマートデバイス(SD)のセッションを抽出し、最後に日毎・デバイス毎のセッション数を集計するクエリは以下のようになります。
row_dataという名前で複数のクエリから使っていること、別のWITHクエリ(pc_sessionとsd_session)内からでも定義済みのrow_dataでクエリできることがポイントです。

#standardSQL

WITH row_data AS (
  SELECT *
  FROM `bigquery-public-data.google_analytics_sample.ga_sessions_*`
  WHERE _TABLE_SUFFIX BETWEEN '20170703' AND '20170730'
),
pc_session AS (
  SELECT date, fullVisitorId, visitId
  FROM row_data
  WHERE device.deviceCategory = 'desktop'
  GROUP BY date, fullVisitorId, visitId 
),
sd_session AS (
  SELECT date, fullVisitorId, visitId
  FROM row_data
  WHERE device.deviceCategory <> 'desktop'
  GROUP BY date, fullVisitorId, visitId 
)
SELECT date, 'PC', COUNT(*) FROM pc_session GROUP BY date
UNION ALL
SELECT date, 'SD', COUNT(*) FROM sd_session GROUP BY date
ORDER BY 1,2

RECORD型をUNNEST関数で展開する

GAのデータを集計する上でやっかいなのが、hitsやcustomDimensionsなどの、繰り返し(repeated)のRECORD型です。

例えば、ランディングページでの滞在時間を集計したいとします。
ランディングページでの滞在時間は、hits.hitNumber=2のときのhits.time(つまりセッション内の2ページ目までの滞在時間)を取得する必要があります。

以下のようにhits.Numberやhits.timeで直接アクセスすることはできません。

#standardSQL

WITH row_data AS (
  SELECT *
  FROM `bigquery-public-data.google_analytics_sample.ga_sessions_*`
  WHERE _TABLE_SUFFIX BETWEEN '20170703' AND '20170730'
)
SELECT date, AVG(hits.time) / 1000
FROM row_data
WHERE row_data.hits.hitNumber = 2
GROUP BY date
ORDER BY 1

こういうときは、以下のようにUNNEST関数でRECORD型を展開します。こうするとRECORD型中の複数の要素が展開され、1要素1行になります (ちょうどINNER JOINで自己結合したようになります)。

#standardSQL

WITH row_data AS (
  SELECT *
  FROM `bigquery-public-data.google_analytics_sample.ga_sessions_*`
  WHERE _TABLE_SUFFIX BETWEEN '20170703' AND '20170730'
)
SELECT date, AVG(hits.time) / 1000
FROM row_data, UNNEST(row_data.hits) AS hits 
WHERE hits.hitNumber = 2
GROUP BY date
ORDER BY 1

試しに、以下のクエリでセッションごとのhitsからhitNumberとtimeを抽出すると、1行1セッションだったデータが、複数行に分かれていることがわかります。

#standardSQL

WITH row_data AS (
  SELECT *
  FROM `bigquery-public-data.google_analytics_sample.ga_sessions_*`
  WHERE _TABLE_SUFFIX BETWEEN '20170703' AND '20170730'
)
SELECT date, fullVisitorId, visitId, hits.hitNumber, hits.time
FROM row_data, UNNEST(row_data.hits) AS hits

SELECT句のサブクエリで集計する

Standard SQLではSELECT句内でサブクエリができます。

例えば、セッション内である商品(hits.product.productSKU='GGOEGBPB021199')が閲覧されたかどうかを出力するクエリは以下のようになります。
SELECT内でサブクエリを発行していますが、サブクエリ内でもUNNESTが使えますので、そこでproductを展開し、productSKUが一致するかどうかを判定しています (こんなふうにCASE WHENで横持ちにして集計するのもあるあるですね)。

#standardSQL

WITH row_data AS (
  SELECT *
  FROM `bigquery-public-data.google_analytics_sample.ga_sessions_*`
  WHERE _TABLE_SUFFIX BETWEEN '20170703' AND '20170730'
)
SELECT date, fullVisitorId, visitId, 
  MAX((SELECT COUNT(DISTINCT productSKU) FROM UNNEST(hits.product)))
FROM row_data, UNNEST(row_data.hits) AS hits 
GROUP BY date, fullVisitorId, visitId
ORDER BY 1, 2, 3

おわりに

自分がGAのデータをBigQueryで集計・分析する際に役立ったテクニックをまとめました。