雑記

参考文献: client_python/metrics.py at master · prometheus/client_python · GitHub

あなたは弁当販売を始めることにしました。A か B の type を指定してリクエストしてもらい、それに応じてA弁当(バターチキンカレー)かB弁当(ローストビーフ)を渡したいと思います。以下のようにすれば、localhost:8080/bentou?type=A にアクセスするとバターチキンカレーが返却されます。

import asyncio
from aiohttp import web

class BentouHandler:
    async def handle(self, request):
        params = {k: v for k, v in request.query.items() if v}
        if 'type' not in params:
            raise ValueError('type がありません')
        type_ = params['type']
        bentou = None
        if type_ == 'A':
            await asyncio.sleep(2)
            bentou = {'主食': 'ナン', 'おかず': 'バターチキンカレー'}
        elif type_ == 'B':
            await asyncio.sleep(1)
            bentou = {'主食': 'ごはん', 'おかず': 'ローストビーフ'}
        else:
            raise ValueError('type が不正です')
        return web.json_response(bentou, status=200)

handler = BentouHandler()
app = web.Application()
app.add_routes([web.get(f'/bentou', handler.handle)])

if __name__ == '__main__':
    web.run_app(app, port=8080)

ここであなたは、これまでに弁当が何個売れたか、実際に弁当を提供するのにかかった時間のヒストグラムはどうなっているかを知りたいと思いました。以下のようにすると localhost:8080/metrics からこれまでの総リクエスト数とレスポンス時間の累積分布が確認できます。

import asyncio
from aiohttp import web
from prometheus_client import Counter, Histogram, generate_latest

class BentouHandler:
    async def handle(self, request):
        params = {k: v for k, v in request.query.items() if v}
        if 'type' not in params:
            raise ValueError('type がありません')
        type_ = params['type']
        bentou = None
        if type_ == 'A':
            await asyncio.sleep(2)
            bentou = {'主食': 'ナン', 'おかず': 'バターチキンカレー'}
        elif type_ == 'B':
            await asyncio.sleep(1)
            bentou = {'主食': 'ごはん', 'おかず': 'ローストビーフ'}
        else:
            raise ValueError('type が不正です')
        return web.json_response(bentou, status=200)

class BentouMetrics:
    def __init__(self):
        # 必要なカウンタを用意
        self.get_req_counter = Counter('get_request_count', 'リクエスト数')
        self.resp_time_hist = Histogram('resp_time_hist', 'レスポンスタイム', buckets=[0.5, 1.5, 2.5])

    @web.middleware
    async def wrapper(self, request, handler):  # ハンドラをラップして各種メトリクスを計測する
        # 目的のエンドポイント以外へのリクエストは計測対象外
        if request.path not in ['/bentou']:
            return await handler(request)

        # 各種メトリクスを計測する
        self.get_req_counter.inc()  # リクエスト数をインクリメント
        with self.resp_time_hist.time():  # レスポンス時間を計測
            response = await handler(request)
        return response

    async def exposer(self, request):  # 現在のメトリクスを返却する
        return web.Response(body=generate_latest(), content_type='text/plain')

    @classmethod
    def setup(cls, app):
        metrics = cls()
        app.middlewares.append(metrics.wrapper)
        app.router.add_get('/metrics', metrics.exposer)

handler = BentouHandler()
app = web.Application()
BentouMetrics.setup(app)
app.add_routes([web.get(f'/bentou', handler.handle)])
...
# HELP get_request_count_total リクエスト数
# TYPE get_request_count_total counter
get_request_count_total 2.0
# HELP resp_time_hist レスポンスタイム
# TYPE resp_time_hist histogram
resp_time_hist_bucket{le="0.5"} 0.0
resp_time_hist_bucket{le="1.5"} 1.0
resp_time_hist_bucket{le="2.5"} 2.0
resp_time_hist_bucket{le="+Inf"} 2.0
...

しかし上記の方法だと正常に弁当を返せなかったときのレスポンス時間もヒストグラムに含まれてしまいます。正常に弁当を返せなかったときは、回数は記録するがレスポンス時間は計測しないようにするには以下のようにすると思います。

import time

class BentouMetrics:
    def __init__(self):
        # 必要なカウンタを用意
        self.get_req_counter = Counter('get_request_count', '正常に返せたリクエスト数')
        self.get_err_req_counter = Counter('get_error_request_count', '正常に返せなかったリクエスト数')
        self.resp_time_hist = Histogram('resp_time_hist', '正常レスポンスタイム', buckets=[0.5, 1.5, 2.5])

    @web.middleware
    async def wrapper(self, request, handler):  # ハンドラをラップして各種メトリクスを計測する
        # 目的のエンドポイント以外へのリクエストは計測対象外
        if request.path not in ['/bentou']:
            return await handler(request)

        # 各種メトリクスを計測する
        time_0 = time.perf_counter()
        try:
            response = await handler(request)
            self.get_req_counter.inc()  # リクエスト数をインクリメント
            self.resp_time_hist.observe(time.perf_counter() - time_0)  # レスポンス時間を記録
        except Exception as e:
            self.get_err_req_counter.inc()  # 正常に返せなかったリクエスト数をインクリメント
            raise e
        return response

雑記: 分布収束関連の話

お気付きの点がありましたらご指摘いただけますと幸いです。

  1. Amazon | Asymptotic Statistics (Cambridge Series in Statistical and Probabilistic Mathematics, Series Number 3) | van der Vaart, A. W. | Applied
    • 本記事はこの 11~12 ページの内容に尾ひれを付けたものである。筆者の誤りは筆者に帰属する。
  1. 雑記 - クッキーの日記( 参考文献 [1] の全体図があるだけ )
  2. 雑記: t分布の話 - クッキーの日記( 母分布が正規分布の場合の t 統計量の話 )

まとめ

  • スラツキーの定理は便利である。
  • 分布収束の収束先が連続分布であるであるならば、分布関数は一様収束する。
 t 統計量の漸近分布(+ α )

 t 統計量とは以下のような、確率変数の関数になっているような確率変数でしたよね。分母の  S_n^2 はいわゆる不偏分散です。母分布の分散  \sigma^2 ではないことに留意してください。

定義〈  t 統計量(※) 〉 Y_1, Y_2, \cdots, Y_n をそれぞれ独立に同一の分布(平均  \mu )にしたがう確率変数とする。このとき、以下の T t 統計量という。
\displaystyle T = T (Y_1, Y_2, \cdots, Y_n) = \frac{\sqrt{n} \bigl( \overline{Y}_n - \mu \bigr)}{\sqrt{S_n^2}} \qquad \left( S_n^2 \equiv \frac{1}{n-1} \sum_{i=1}^n (Y_i - \overline{Y}_n)^2 \right)
※ 例えば現代数理統計学(竹村)では明示的に「正規分布からの無作為標本(中略)を t 統計量という」とあるので正規分布からのサンプルのときに t 統計量というのだという感じにみえます。他方、逆に Wikipedia では t 統計量を単に何らかのパラメータの推定量とその推定量標準偏差で構成されるものと定義しています。「独立同一分布にしたがうサンプルから分布の平均を推定する」という文脈に限っていません。上にかいた定義は [参考文献 1] の定義がこれだろうかと思ってかいたものですが、よく読むと「独立同一分布にしたがうサンプルから分布の平均を推定する」という文脈のときのみ t 統計量というとはどこにもかいていない気がするので文脈を制限しすぎかもしれません。何にせよ、t 統計量の話をするときは相手と定義を確認してください。
特に Y_1, Y_2, Y_3, \cdots, Y_n がしたがう「同一の分布」が正規分布 \mathcal{N}(\mu, \sigma^2) であるとき、T は自由度 n-1 t 分布にしたがうのでした [関連記事 イ]。 t 統計量は母分布の分散を使用しないので、母分散がわからないときの母平均の値の存在範囲に関する仮説の検定などに用いられますね。つまり、「もし母平均が \mu = \mu_0 であるならば、T はある  t 分布にしたがうので、T の値がその分布の真ん中付近からあまりに外れていたら \mu = \mu_0 であるか疑わしい」という論法ですね。

では、「同一の分布」が正規分布ではない場合はどうだろう。

正規分布ではない場合ですか? そういわれましても、[関連記事 イ] で  t 統計量がしたがう分布を導出したときは全体的に母分布が正規分布であることを利用していましたし、分布について何もわからなかったら同じ要領でアプローチすることはできませんね……。

うん。だから、n がとても大きいときを考えてみよう。「n が大きいときは正規分布に近づく」というような定理があったよね?

n が大きいときは正規分布に近づく」……中心極限定理ですね。

定理〈 中心極限定理 Y_1, Y_2, \cdots, Y_n をそれぞれ独立に同一の分布(平均  \mu,分散 \sigma^2 < \infty )にしたがう確率変数とする。このとき、 \sqrt{n} \bigl( \overline{Y}_n - \mu \bigr)\mathcal{N}(0, \sigma^2) に分布収束する。
なるほど、この  \sqrt{n} \bigl( \overline{Y}_n - \mu \bigr) t 統計量の分子にあたる部分です。中心極限定理によれば、母分布の分散が有限でさえあれば、その部分は正規分布に収束していくというわけですね。しかし、分母の  \sqrt{S_n^2} も確率変数ですから、 t 統計量全体としてどのような分布にしたがうかは依然としてわかりません。

そうだね。ここで、いまは  \mu=0 であるということしよう。それで、後々の都合上、以下の2つの定理をつかって、 \sqrt{S_n^2} が何にどんな収束をするかを考えてみてほしい。

定理〈 大数の弱法則(※) 〉 Y_1, Y_2, \cdots, Y_n をそれぞれ独立に同一の分布(平均  \mu < \infty,分散 \sigma^2 < \infty )にしたがう確率変数とする。このとき、 \overline{Y}_n\mu に確率収束する。
\sigma^2 < \infty がなくても成り立ちますがこの条件込みで弱法則としてチェビシェフの不等式で証明することが多いと思います(例.清水本)。ただ後から気付いたのですが、[参考文献 1] の15ページではこの条件がないものを大数の弱法則とよんでいました。証明はそのページをみてください。
定理〈 連続写像定理(確率収束版) 〉 g: \mathbb{R} \to \mathbb{R}P(Y \in C) = 1 であるような C 上の任意の点で連続な関数であるとする。このとき、以下が成り立つ。
Y_n \xrightarrow{\mathrm{P}} Y \quad \Rightarrow \quad g(Y_n) \xrightarrow{\mathrm{P}} g(Y)

 \mu=0 とするのですか? そして誘導が強い……まずは [関連記事 イ] と同様に  S_n^2 を展開しましょう。

 \displaystyle \begin{split} S_n^2 &= \frac{n}{n -1} \left( \frac{1}{n} \sum_{i=1}^{n} Y_{i}^2 - \bar{Y_n}^2 \right) \equiv \frac{n}{n -1} \left( A_n - \bar{Y_n}^2 \right) \end{split}
いま、大数の弱法則より \overline{Y}_n \xrightarrow{\mathrm{P}} \mu = 0 であり、また、A_n \xrightarrow{\mathrm{P}}〈母分布の2乗平均〉  = \sigma^2 です。よって、連続写像定理も用いて、A_n - \bar{Y_n}^2 \xrightarrow{\mathrm{P}} \sigma^2 です。これを確率収束の定義に差し戻してかくと以下の (*) ですが、 (*) が成り立つならば (**) も成り立つとわかります。
\displaystyle \forall \varepsilon, \; \, \lim_{n \to \infty} P \bigl( \{ \omega \in \Omega \, | \, \| A_n - \bar{Y_n}^2 - \sigma^2 \| > \varepsilon \} \bigr) = 0 \tag{*}
\displaystyle \forall \varepsilon, \; \, \lim_{n \to \infty} P \bigl( \{ \omega \in \Omega \, | \, \| \frac{n}{n -1} \left( A_n - \bar{Y_n}^2 \right) - \sigma^2 \| > \varepsilon \} \bigr) = 0  \tag{**}
つまり、 S_n^2\sigma^2 に確率収束します。であれば連続写像定理より、 \sqrt{S_n^2}\sqrt{\sigma^2} に確率収束しますね。

ここまでをまとめると以下です。

  • 独立に同一の分布(平均  \mu = 0,分散 \sigma^2 < \infty )にしたがう Y_1, Y_2, \cdots, Y_n に対する  t 統計量について以下のことがいえる。
    • 分子の  \sqrt{n} \bigl( \overline{Y}_n \bigr)\mathcal{N}(0, \sigma^2) に分布収束する。
    • 分母の  \sqrt{S_n^2}\sqrt{\sigma^2} に確率収束する。
……「分子はこう」「分母はこう」とわかりましたが、だから何なのでしょうか。いま分子と分母が独立であるのかもわかりませんし(母分布が正規分布の場合は独立であることが示せましたが [関連記事 イ] )。やはり全体としての確率変数がどうふるまうかわからないと思うのですが。

ここまできたら以下の補題を適用することができるよ。

定理〈 スラツキーの定理(除算版) 〉 X_nX に分布収束し、Y_nc(ゼロでない定数)に分布収束するとする。このとき、 Y_n^{-1} X_nc^{-1} X に分布収束する。

なんと、そのような補題が……確率収束するならば分布収束しますから、スラツキーの定理を現在の状況に適用すれば、結局、 t 統計量は \mathcal{N}(0, 1) に分布収束することになりますね。平均がゼロで分散が有限な同一分布から独立なサンプルをどんどん得ていくとき、その  t 統計量は標準正規分布に分布収束するといえるのですか……。

スラツキーの定理は便利だよね。スラツキーの定理をつかうと以下のような例も扱えるよ。

例.
  • あなたは毎日何らかの確率変数(たち)を観測しているとする。 \theta, \, \sigma^2 は毎日観測する確率変数(たち)がしたがう分布の何らかのパラメータである。あなたは特に  \theta の真の値に興味がある。が、 \theta を直接観測することはできず、観測できる変数から推測するしかない。
  • ここで、n 日目までの観測結果から以下のような統計量  T_n, \, S_n^2 を構成できるとわかった。
    •  \sqrt{n} \bigl( T_n - \theta \bigr)\mathcal{N}(0, \sigma^2) に分布収束する。
    •  S_n^2\sigma^2 に確率収束する。
  • このとき、スラツキーの定理より、 Z_n \equiv \sqrt{n} \bigl( T_n - \theta \bigr) / \sqrt{S_n^2}\mathcal{N}(0, 1) に分布収束する。
  • 上の Z_n が標準正規分布にしたがうとすれば、あなたは \theta の信頼区間を求めることができるだろう。
要するに、知りたいパラメータ  \theta を中心に分布する T_n を構成できたとして、肝心のばらつきがわからないと結局信頼区間がわからないよね。でも、ばらつきに近づく S_n^2 も構成できたなら、信頼区間がわかるんだ。これはスラツキーの定理を適用しないといえないよ。そうでないと部長がさっきいったように、「分子はこう、分母はこう、だから何?」という話だからね。

はあ……うーん、あの、基礎的な統計学で信頼区間を求めたり、統計的仮説検定をしたりするときは、「この統計量はこの分布にしたがう」がきちんと求まっていましたよね? 対していまわかっているのは、「この統計量はこの分布に『分布収束する』」です。もちろん分布収束するのであれば観測をたくさん重ねていけばその分布に近づきはするでしょうが……しかし、確か分布収束の定義は「分布関数の『各点』が収束すればよい」といったものでしたよね? それだと、観測をある程度たくさん重ねたとき、統計量の分布が収束先の分布に近い形をしている保証はないのではないですか? 各点収束だったら、「こちらのエリアではもう収束先の分布に近いですが、あちらのエリアはまだまだ近くないですね」ということもあり得るでしょう? そのような状況かもしれないのに信頼区間を議論するのは乱暴なのでは?

実際そういう心配はあるかもね。ただ、ここまでの例については一応、分布関数は一様収束しているけどね。

あれ、そうなんですか?

収束先が連続分布であるときの分布収束

連続分布なら以下が成り立つので示してみよう( [関連記事 1] の12ページ)。

補題〈 収束先が連続分布であるときの分布収束 〉 k 次元確率変数 X_n は分布関数が連続であるような X に分布収束するとする。このとき、 \sup_x \bigl| P(X_n \leqq x) - P(X \leqq x) \bigr| \to 0 が成り立つ。

以下のように示せますね……確かに連続分布であれば一様収束します(2022-11-14 追記: 以下の手書きノートでは1次元の場合も2次元の場合も「N'=* をとって、k=* とすれば、」とありますが、N' は k に依存するので「k=* として、N'=* をとれば、」が正しいです)


連続分布であれば、領域を 100 分割して、個々のエリアに対応する確率が 1/100 というようにできますね(連続分布でない場合はできる保証がありません)。そして、個々の領域での差  \bigl| F_n(x) - F(x) \bigr| はどちらかの端っこでの差 + 1/100 で上から抑えられますね。各分割領域の端っこは有限個ですから、無限個の点の上での話( \sup_x)を、有限個(いまは99個)の点の上での話にできるわけです。99個の点を先に固定すれば、各点では収束しますから、すべての点での差の上界が十分小さくなるように n を大きくすることが可能です。任意の正数に応じて領域の分割数をじゅうぶん細かくして、何個かの点を固定して、n をじゅうぶん進めればよいですね。2次元以上の場合も同様の議論が可能です。上図で水色に塗っている領域は 2 × 1/5 で上から抑えられますが、d 次元であれば d × 1/5 で上から抑えられます。

終わり