「やりなおし統計」一覧

[mathjax] Record Linkage、Entity Recognitionなど、いわゆる「名寄せ」問題には、割とさまざまな解決策が 提案されている。その1つに確率論的な数学的背景を持つFellegi-Sunterモデルを実装したSplinkがある。 数学的背景について咀嚼して説明してくれているサイトはあまりなく、Splinkの公式ドキュメントが唯一 に近いソースだろうと思う。 今回は5万人の人物データを使って重複を除去するSplinkの主要な機能のサンプルを追ってみた。 数式を読んだところで実りは薄い、ということもあり、なるべく数式を読まずに読み進めたが、 モデルのパラメタ、予測結果を解釈する際にどうしても数式を読まないとわからない箇所があった。 特にベイズ統計について理解がかなり浅くて、正直なところ理解しきれてはいない。 Splinkの公式ドキュメントを読み進めていった際に気付いたメモを書いてみることにした。 [arst_toc tag=\"h4\"] 環境構築・準備 環境構築の実装部分は主な学習要素ではないのでClaude Codeでサクッと作ってみた。 まずはPlan modeで計画を立ててもらう。これぐらいだとほぼ言い値で使える。 DuckDB、Splinkを使う最小構成の環境を作ります。環境はDocker composeで構築します。計画を立ててください。 ⏺ DuckDB、Splinkを使うDocker Compose環境の計画を完了しました。 計画概要 1. 設計文書: .tmp/design.mdで要件、アーキテクチャ、技術仕様を定義 2. タスク管理: .tmp/task.mdで詳細な実装手順を3フェーズに分けて整理 主要コンポーネント - Jupyter Lab環境: DuckDB、Splink、必要なライブラリがプリインストール - データボリューム: サンプルデータ、notebook、出力用のディレクトリ - 最小構成: 開発・学習用途に特化した軽量環境 次のステップ 実装を開始する場合は、Phase 1から順番に進めます： 1. Docker設定ファイル作成 2. ディレクトリ構造構築 3. サンプル実装 4. テストと検証 計画の詳細は.tmp/design.mdと.tmp/task.mdをご確認ください。実装を開始しますか？ 実際にファイルを作ってもらうと以下のような構成となる。 得られる中間成果物である docker-compose.yml、Dockerfile は本記事では省略する。 . ├── data ├── docker-compose.yml ├── Dockerfile ├── notebooks ├── output └── requirements.txt 普通に http://localhost:8888 で JupyterLab が開く。 使用するサンプルデータ 5万人の人物データを使って名寄せを行うサンプル。おそらくSplinkの用途として最初に思いつくやつ。 Splinkにデータをロードする前に必要なデータクリーニング手順について説明がある。 公式によると、まずは行に一意のIDを割り当てる必要がある。 データセット内で一意となるIDであって、重複除去した後のエンティティを識別するIDのことではない。 [clink implicit=\"false\" url=\"https://moj-analytical-services.github.io/splink/demos/tutorials/01_Prerequisites.html\" imgurl=\"https://user-images.githubusercontent.com/7570107/85285114-3969ac00-b488-11ea-88ff-5fca1b34af1f.png\" title=\"Data Prerequisites\" excerpt=\"Splink では、リンクする前にデータをクリーンアップし、行に一意の ID を割り当てる必要があります。このセクションでは、Splink にデータをロードする前に必要な追加のデータクリーニング手順について説明します。\"] 使用するサンプルデータは以下の通り。 from splink import splink_datasets df = splink_datasets.historical_50k df.head() データの分布を可視化 splink.exploratoryのprofile_columnsを使って分布を可視化してみる。 from splink import DuckDBAPI from splink.exploratory import profile_columns db_api = DuckDBAPI() profile_columns(df, db_api, column_expressions=[\"first_name\", \"substr(surname,1,2)\"]) 同じ姓・名の人が大量にいることがわかる。 ブロッキングとブロッキングルールの評価 テーブル内のレコードが他のレコードと「同一かどうか」を調べるためには、 基本的には、他のすべてのレコードとの何らかの比較操作を行うこととなる。 全てのレコードについて全てのカラム同士を比較したいのなら、 対象のテーブルをCROSS JOINした結果、各カラム同士を比較することとなる。 SELECT ... FROM input_tables as l CROSS JOIN input_tables as r あるカラムが条件に合わなければ、もうその先は見ても意味がない、 というケースは多い。例えば、まず first_name 、surname が同じでなければ、 その先の比較を行わない、というのはあり得る。 SELECT ... FROM input_tables as l INNER JOIN input_tables as r ON l.first_name = r.first_name AND l.surname = r.surname このような考え方をブロッキング、ON句の条件をブロッキングルールと言う。 ただ、これだと性と名が完全一致していないレコードが残らない。 そこで、ブロッキングルールを複数定義し、いずれかが真であれば残すことができる。 ここでポイントなのが、ブロッキングルールを複数定義したとき、 それぞれのブロッキングルールで重複して選ばれるレコードが発生した場合、 Splinkが自動的に排除してくれる。 このため、ブロッキングルールを重ねがけすると、最終的に残るレコード数は一致する。 ただ、順番により、同じルールで残るレコード数は変化する。 逆に言うと、ブロッキングルールを足すことで、重複除去後のOR条件が増えていく。 積算グラフにして、ブロッキングルールとその順番の効果を見ることができる。 from splink import DuckDBAPI, block_on from splink.blocking_analysis import ( cumulative_comparisons_to_be_scored_from_blocking_rules_chart, ) blocking_rules = [ block_on(\"substr(first_name,1,3)\", \"substr(surname,1,4)\"), block_on(\"surname\", \"dob\"), block_on(\"first_name\", \"dob\"), block_on(\"postcode_fake\", \"first_name\"), block_on(\"postcode_fake\", \"surname\"), block_on(\"dob\", \"birth_place\"), block_on(\"substr(postcode_fake,1,3)\", \"dob\"), block_on(\"substr(postcode_fake,1,3)\", \"first_name\"), block_on(\"substr(postcode_fake,1,3)\", \"surname\"), block_on(\"substr(first_name,1,2)\", \"substr(surname,1,2)\", \"substr(dob,1,4)\"), ] db_api = DuckDBAPI() cumulative_comparisons_to_be_scored_from_blocking_rules_chart( table_or_tables=df, blocking_rules=blocking_rules, db_api=db_api, link_type=\"dedupe_only\", ) 積算グラフは以下の通り。積み上がっている数値は「比較の数」。 要は、論理和で条件を足していって、次第に緩和されている様子がわかる。 DuckDBでは比較の数を2,000万件以内、Athena,Sparkでは1億件以内を目安にせよとのこと。 比較の定義 Splinkは Fellegi-Sunter model モデル (というかフレームワーク) に基づいている。 https://moj-analytical-services.github.io/splink/topic_guides/theory/fellegi_sunter.html 各カラムの同士をカラムの特性に応じた距離を使って比較し、重みを計算していく。 各カラムの比較に使うためのメソッドが予め用意されているので、特性に応じて選んでいく。 以下では、first_name, sur_name に ForenameSurnameComparison が使われている。 dobにDateOfBirthComparison、birth_place、ocupationにExactMatchが使われている。 import splink.comparison_library as cl from splink import Linker, SettingsCreator settings = SettingsCreator( link_type=\"dedupe_only\", blocking_rules_to_generate_predictions=blocking_rules, comparisons=[ cl.ForenameSurnameComparison( \"first_name\", \"surname\", forename_surname_concat_col_name=\"first_name_surname_concat\", ), cl.DateOfBirthComparison( \"dob\", input_is_string=True ), cl.PostcodeComparison(\"postcode_fake\"), cl.ExactMatch(\"birth_place\").configure(term_frequency_adjustments=True), cl.ExactMatch(\"occupation\").configure(term_frequency_adjustments=True), ], retain_intermediate_calculation_columns=True, ) # Needed to apply term frequencies to first+surname comparison df[\"first_name_surname_concat\"] = df[\"first_name\"] + \" \" + df[\"surname\"] linker = Linker(df, settings, db_api=db_api) ComparisonとComparison Level ここでSplinkツール内の比較の概念の説明。以下の通り概念に名前がついている。 Data Linking Model ├─-- Comparison: Date of birth │ ├─-- ComparisonLevel: Exact match │ ├─-- ComparisonLevel: One character difference │ ├─-- ComparisonLevel: All other ├─-- Comparison: First name │ ├─-- ComparisonLevel: Exact match on first_name │ ├─-- ComparisonLevel: first_names have JaroWinklerSimilarity > 0.95 │ ├─-- ComparisonLevel: first_names have JaroWinklerSimilarity > 0.8 │ ├─-- ComparisonLevel: All other モデルのパラメタ推定 モデルの実行に必要なパラメタは以下の3つ。Splinkを用いてパラメタを得る。 ちなみに u は \"\'U\'nmatch\"、m は \"\'M\'atch\"。背後の数式の説明で現れる。 No パラメタ 説明 1 無作為に選んだレコードが一致する確率 入力データからランダムに取得した2つのレコードが一致する確率 (通常は非常に小さい数値) 2 u値(u確率) 実際には一致しないレコードの中で各 ComparisonLevel に該当するレコードの割合。具体的には、レコード同士が同じエンティティを表すにも関わらず値が異なる確率。例えば、同じ人なのにレコードによって生年月日が違う確率。これは端的には「データ品質」を表す。名前であればタイプミス、別名、ニックネーム、ミドルネーム、結婚後の姓など。 3 m値(m確率) 実際に一致するレコードの中で各 ComparisonLevel に該当するレコードの割合。具体的には、レコード同士が異なるエンティティを表すにも関わらず値が同じである確率。例えば別人なのにレコードによって性・名が同じ確率 (同姓同名)。性別は男か女かしかないので別人でも50%の確率で一致してしまう。 無作為に選んだレコードが一致する確率 入力データからランダムに抽出した2つのレコードが一致する確率を求める。 値は0.000136。すべての可能なレコードのペア比較のうち7,362.31組に1組が一致すると予想される。 合計1,279,041,753組の比較が可能なため、一致するペアは合計で約173,728.33組になると予想される、 とのこと。 linker.training.estimate_probability_two_random_records_match( [ block_on(\"first_name\", \"surname\", \"dob\"), block_on(\"substr(first_name,1,2)\", \"surname\", \"substr(postcode_fake,1,2)\"), block_on(\"dob\", \"postcode_fake\"), ], recall=0.6, ) > Probability two random records match is estimated to be 0.000136. > This means that amongst all possible pairwise record comparisons, > one in 7,362.31 are expected to match. > With 1,279,041,753 total possible comparisons, > we expect a total of around 173,728.33 matching pairs u確率の推定 実際には一致しないレコードの中でComparisonの評価結果がPositiveである確率。 基本、無作為に抽出したレコードは一致しないため、「無作為に抽出したレコード」を 「実際には一致しないレコード」として扱える、という点がミソ。 probability_two_random_records_match によって得られた値を使ってu確率を求める。 estimate_u_using_random_sampling によって、ラベルなし、つまり教師なしでu確率を得られる。 レコードのペアをランダムでサンプルして上で定義したComparisonを評価する。 ランダムサンプルなので大量の不一致が発生するが、各Comparisonにおける不一致の分布を得ている。 これは、例えば性別について、50%が一致、50%が不一致である、という分布を得ている。 一方、例えば生年月日について、一致する確率は 1%、1 文字の違いがある確率は 3%、 その他はすべて 96% の確率で発生する、という分布を得ている。 linker.training.estimate_u_using_random_sampling(max_pairs=5e6) > ----- Estimating u probabilities using random sampling ----- > > Estimated u probabilities using random sampling > > Your model is not yet fully trained. Missing estimates for: > - first_name_surname (no m values are trained). > - dob (no m values are trained). > - postcode_fake (no m values are trained). > - birth_place (no m values are trained). > - occupation (no m values are trained). m確率の推定 「実際に一致するレコード」の中で、Comparisonの評価がNegativeになる確率。 そもそも、このモデルを使って名寄せ、つまり「一致するレコード」を見つけたいのだから、 モデルを作るために「実際に一致するレコード」を計算しなければならないのは矛盾では..となる。 無作為抽出結果から求められるu確率とは異なり、m確率を求めるのは難しい。 もしラベル付けされた「一致するレコード」、つまり教師データセットがあるのであれば、 そのデータセットを使ってm確率を求められる。 例えば、日本人全員にマイナンバーが振られて、全てのレコードにマイナンバーが振られている、 というアナザーワールドがあるのであれば、マイナンバーを使ってm確率を推定する。(どういう状況??) ラベル付けされたデータがないのであれば、EMアルゴリズムでm確率を求めることになっている。 EMアルゴリズムは反復的な手法で、メモリや収束速度の点でペア数を減らす必要があり、 例ではブロッキングルールを設定している。 以下のケースでは、first_nameとsurnameをブロッキングルールとしている。 つまり、first_name, surnameが完全に一致するレコードについてペア比較を行う。 この仮定を設定したため、first_name, surname (first_name_surname) のパラメタを推定できない。 training_blocking_rule = block_on(\"first_name\", \"surname\") training_session_names = ( linker.training.estimate_parameters_using_expectation_maximisation( training_blocking_rule, estimate_without_term_frequencies=True ) ) > ----- Starting EM training session ----- > > Estimating the m probabilities of the model by blocking on: > (l.\"first_name\" = r.\"first_name\") AND (l.\"surname\" = r.\"surname\") > > Parameter estimates will be made for the following comparison(s): > - dob > - postcode_fake > - birth_place > - occupation > > Parameter estimates cannot be made for the following comparison(s) since they are used in the blocking rules: > - first_name_surname > > Iteration 1: Largest change in params was 0.248 in probability_two_random_records_match > Iteration 2: Largest change in params was 0.0929 in probability_two_random_records_match > Iteration 3: Largest change in params was -0.0237 in the m_probability of birth_place, level `Exact match on > birth_place` > Iteration 4: Largest change in params was 0.00961 in the m_probability of birth_place, level `All other >comparisons` > Iteration 5: Largest change in params was -0.00457 in the m_probability of birth_place, level `Exact match on birth_place` > Iteration 6: Largest change in params was -0.00256 in the m_probability of birth_place, level `Exact match on birth_place` > Iteration 7: Largest change in params was 0.00171 in the m_probability of dob, level `Abs date difference Iteration 8: Largest change in params was 0.00115 in the m_probability of dob, level `Abs date difference Iteration 9: Largest change in params was 0.000759 in the m_probability of dob, level `Abs date difference Iteration 10: Largest change in params was 0.000498 in the m_probability of dob, level `Abs date difference Iteration 11: Largest change in params was 0.000326 in the m_probability of dob, level `Abs date difference Iteration 12: Largest change in params was 0.000213 in the m_probability of dob, level `Abs date difference Iteration 13: Largest change in params was 0.000139 in the m_probability of dob, level `Abs date difference Iteration 14: Largest change in params was 9.04e-05 in the m_probability of dob, level `Abs date difference <= 10 year` 同様にdobをブロッキングルールに設定して実行すると、dob以外の列についてパラメタを推定できる。 training_blocking_rule = block_on(\"dob\") training_session_dob = ( linker.training.estimate_parameters_using_expectation_maximisation( training_blocking_rule, estimate_without_term_frequencies=True ) ) > ----- Starting EM training session ----- > > Estimating the m probabilities of the model by blocking on: > l.\"dob\" = r.\"dob\" > > Parameter estimates will be made for the following comparison(s): > - first_name_surname > - postcode_fake > - birth_place > - occupation > > Parameter estimates cannot be made for the following comparison(s) since they are used in the blocking rules: > - dob > > Iteration 1: Largest change in params was -0.474 in the m_probability of first_name_surname, level `Exact match on first_name_surname_concat` > Iteration 2: Largest change in params was 0.052 in the m_probability of first_name_surname, level `All other comparisons` > Iteration 3: Largest change in params was 0.0174 in the m_probability of first_name_surname, level `All other comparisons` > Iteration 4: Largest change in params was 0.00532 in the m_probability of first_name_surname, level `All other comparisons` > Iteration 5: Largest change in params was 0.00165 in the m_probability of first_name_surname, level `All other comparisons` > Iteration 6: Largest change in params was 0.00052 in the m_probability of first_name_surname, level `All other comparisons` > Iteration 7: Largest change in params was 0.000165 in the m_probability of first_name_surname, level `All other comparisons` > Iteration 8: Largest change in params was 5.29e-05 in the m_probability of first_name_surname, level `All other comparisons` > > EM converged after 8 iterations > > Your model is not yet fully trained. Missing estimates for: > - first_name_surname (some u values are not trained). モデルパラメタの可視化 m確率、u確率の可視化。 マッチウェイトの可視化。マッチウェイトは (log_2 (m / u))で計算される。 linker.visualisations.match_weights_chart() モデルの保存と読み込み 以下でモデルを保存できる。 settings = linker.misc.save_model_to_json( \"./saved_model_from_demo.json\", overwrite=True ) 以下で保存したモデルを読み込める。 import json settings = json.load( open(\'./saved_model_from_demo.json\', \'r\') ) リンクするのに十分な情報が含まれていないレコード 「John Smith」のみを含み、他のすべてのフィールドがnullであるレコードは、 他のレコードにリンクされている可能性もあるが、潜在的なリンクを明確にするには十分な情報がない。 以下により可視化できる。 linker.evaluation.unlinkables_chart() 横軸は「マッチウェイトの閾値」。縦軸は「リンクするのに十分な情報が含まれないレコード」の割合。 マッチウェイト閾値=6.11ぐらいのところを見ると、入力データセットのレコードの約1.3%が リンクできないことが示唆される。 訓練済みモデルを使って未知データのマッチウェイトを予測 上で構築した推定モデルを使用し、どのペア比較が一致するかを予測する。 内部的には以下を行うとのこと。 blocking_rules_to_generate_predictionsの少なくとも1つと一致するペア比較を生成 Comparisonで指定されたルールを使用して、入力データの類似性を評価 推定された一致重みを使用し、要求に応じて用語頻度調整を適用して、最終的な一致重みと一致確率スコアを生成 df_predictions = linker.inference.predict(threshold_match_probability=0.2) df_predictions.as_pandas_dataframe(limit=1) > Blocking time: 0.88 seconds > Predict time: 1.91 seconds > > -- WARNING -- > You have called predict(), but there are some parameter estimates which have neither been estimated or > specified in your settings dictionary. To produce predictions the following untrained trained parameters will > use default values. > Comparison: \'first_name_surname\': > u values not fully trained records_to_plot = df_e.to_dict(orient=\"records\") linker.visualisations.waterfall_chart(records_to_plot, filter_nulls=False) predictしたマッチウェイトの可視化、数式との照合 predictしたマッチウェイトは、ウォーターフォール図で可視化できる。 マッチウェイトは、モデル内の各特徴量によって一致の証拠がどの程度提供されるかを示す中心的な指標。 (lambda)は無作為抽出した2つのレコードが一致する確率。(K=m/u)はベイズ因子。 begin{align} M &= log_2 ( frac{lambda}{1-lambda} ) + log_2 K \\ &= log_2 ( frac{lambda}{1-lambda} ) + log_2 m - log_2 u end{align} 異なる列の比較が互いに独立しているという仮定を置いていて、 2つのレコードのベイズ係数が各列比較のベイズ係数の積として扱う。 begin{eqnarray} K_{feature} = K_{first_name_surname} + K_{dob} + K_{postcode_fake} + K_{birth_place} + K_{occupation} + cdots end{eqnarray} マッチウェイトは以下の和。 begin{eqnarray} M_{observe} = M_{prior} + M_{feature} end{eqnarray} ここで begin{align} M_{prior} &= log_2 (frac{lambda}{1-lambda}) \\ M_{feature} &= M_{first_name_surname} + M_{dob} + M_{postcode_fake} + M_{birth_place} + M_{occupation} + cdots end{align} 以下のように書き換える。 begin{align} M_{observe} &= log_2 (frac{lambda}{1-lambda}) + sum_i^{feature} log_2 (frac{m_i}{u_i}) \\ &= log_2 (frac{lambda}{1-lambda}) + log_2 (prod_i^{feature} (frac{m_i}{u_i}) ) end{align} ウォーターフォール図の一番左、赤いバーは(M_{prior} = log_2 (frac{lambda}{1-lambda}))。 特徴に関する追加の知識が考慮されていない場合のマッチウェイト。 横に並んでいる薄い緑のバーは (M_{first_name_surname} + M_{dob} + M_{postcode_fake} + M_{birth_place} + M_{occupation} + cdots)。 各特徴量のマッチウェイト。 一番右の濃い緑のバーは2つのレコードの合計マッチウェイト。 begin{align} M_{feature} &= M_{first_name_surname} + M_{dob} + M_{postcode_fake} + M_{birth_place} + M_{occupation} + cdots \\ &= 8.50w end{align} まとめ 長くなったのでいったん終了。この記事では教師なし確率的名寄せパッケージSplinkを使用してモデルを作ってみた。 次の記事では、作ったモデルを使用して実際に名寄せをしてみる。 途中、DuckDBが楽しいことに気づいたので、DuckDBだけで何個か記事にしてみようと思う。

[mathjax] 分散と標準偏差を計算しやすく変形できる。 いちいち偏差(x_i-bar{x})を計算しておかなくても、2乗和(x_i^2)と平均(bar{x})がわかっていればOK。 begin{eqnarray} s^2 &=& frac{1}{n} sum_{i=1}^n (x_i - bar{x})^2 \\ &=& frac{1}{n} sum_{i=1}^n ( x_i^2 -2 x_i bar{x} + bar{x}^2 ) \\ &=& frac{1}{n} ( sum_{i=1}^n x_i^2 -2 bar{x} sum_{i=1}^n x_i + bar{x}^2 sum_{i=1}^n 1) \\ &=& frac{1}{n} ( sum_{i=1}^n x_i^2 -2 n bar{x}^2 + nbar{x}^2 ) \\ &=& frac{1}{n} ( sum_{i=1}^n x_i^2 - nbar{x}^2 ) \\ &=& frac{1}{n} sum_{i=1}^n x_i^2 - bar{x}^2 \\ s &=& sqrt{frac{1}{n} sum_{i=1}^n x_i^2 - bar{x}^2 } end{eqnarray} 以下みたい使える。平均と標準偏差と2乗和の関係。 begin{eqnarray} sum_{i=1}^n (x_i - bar{x})^2 &=& sum_{i=1}^n x_i^2 - nbar{x}^2 \\ ns^2 &=& sum_{i=1}^n x_i^2 - nbar{x}^2 \\ sum_{i=1}^n x_i^2 &=& n(s^2 + bar{x} ) end{eqnarray}

[mathjax] 2項分布に従う母集団の母平均を推測するために有意水準を設定して95%信頼区間を求めてみた。 母平均のあたりがついていない状況だとやりにくい。 [clink url=\"https://ikuty.com/2019/01/11/sampling/\"] (hat{p})がどんな値であっても下限は(hat{p})の関数で抑えられると思ったので、 気になって(hat{p})を変数のまま残すとどうなるかやってみた。 begin{eqnarray} 1.96sqrt{frac{hat{p}(1-hat{p})}{n}} le 0.05 \\ frac{1.96}{0.05}sqrt{hat{p}(1-hat{p})} le sqrt{n} \\ 39.2^2 hat{p}(1-hat{p}) le n end{eqnarray} 左辺を(f(hat{p}))と置くと (f(hat{p}))は下に凸の2次関数であって、 (frac{d}{dhat{p}}f(hat{p})=0)の時に最大となる。というか(hat{p}=0.5)。 (hat{p}=0.5)であるとすると、これはアンケートを取るときのサンプル数を求める式と同じで、 非常に有名な以下の定数が出てくる。 begin{eqnarray} 1537 * 0.5 (1-0.5) le n \\ 384 le n end{eqnarray} (hat{p})がどんな値であっても、サンプル数を400とれば、 有意水準=5%の95%信頼区間を得られる。 だから、アンケートの(n)数はだいたい400で、となる。 さらに、有意水準を10%にとれば、(n)の下限は100で抑えられる。 なるはやのアンケートなら100、ちゃんとやるには400、というやつがこれ。

[mathjax] ちょっと不思議な計算をしてみる。 仮定に仮定を積み重ねた素人の統計。 成功か失敗かを応答する認証装置があったとする。 1回の試行における成功確率(p)は試行によらず一定でありベルヌーイ試行である。 (n)回の独立な試行を繰り返したとき、成功数(k)を確率変数とする離散確率変数に従う。 二項分布の確率密度関数は以下の通り。 begin{eqnarray} P(X=k)= {}_n C_k p^k (1-p)^{n-k} end{eqnarray} 期待値、分散は、 begin{eqnarray} E(X) &=& np \\ V(X) &=& np(1-p) end{eqnarray} (z)得点（偏差値,つまり平均からの誤差が標準偏差何個分か？）は、 begin{eqnarray} z &=& frac{X-E(X)}{sigma} \\ &=& frac{X-E(X)}{sqrt{V(X)}} \\ &=& frac{X-np}{sqrt{np(1-p)}} end{eqnarray} であり、(z)は標準正規分布に従う。 これを標本比率(hat{p}=frac{X}{n})を使うように式変形する。 begin{eqnarray} z &=& frac{frac{1}{n}}{frac{1}{n}} frac{X-np}{sqrt{np(1-p)}} \\ &=& frac{frac{X}{n}-p}{sqrt{frac{p(1-p)}{n}}} \\ &=& frac{hat{p}-p}{sqrt{frac{p(1-p)}{n}}} end{eqnarray} (n)が十分に大きいとき、(z)は標準正規分布(N(0,1))に従う。 従って、(Z)の95%信頼区間は以下である。 begin{eqnarray} -1.96 le Z le 1.96 end{eqnarray} なので、 begin{eqnarray} -1.96 le frac{hat{p}-p}{sqrt{frac{p(1-p)}{n}}} le 1.96 end{eqnarray} (hat{p})は(p)の一致推定量であるから、(n)が大なるとき(hat{p}=p)とすることができる。 begin{eqnarray} -1.96 le frac{hat{p}-p}{sqrt{frac{hat{p}(1-hat{p})}{n}}} le 1.96 \\ end{eqnarray} (p)について解くと(p)の95%信頼区間が求まる。 begin{eqnarray} hat{p}-1.96 sqrt{frac{hat{p}(1-hat{p})}{n}} le p le hat{p}+1.96 sqrt{frac{hat{p}(1-hat{p})}{n}} end{eqnarray} 上記のにおいて、標準誤差(1.96sqrt{frac{hat{p}(1-hat{p})}{n}})が小さければ小さいほど、 95%信頼区間の幅が狭くなる。この幅が5%以内であることを言うためには以下である必要がある。 (有意水準=5%) begin{eqnarray} 1.96sqrt{frac{hat{p}(1-hat{p})}{n}} le 0.05 end{eqnarray} 観測された(hat{p})が(0.9)であったとして(n)について解くと、 begin{eqnarray} 1.96sqrt{frac{0.9(1-0.9)}{n}} le 0.05 \\ frac{1.96}{0.05} sqrt{0.09} le sqrt{n} \\ 11.76 le sqrt{n} \\ 138.2 le n end{eqnarray} 139回試行すれば、100回中95回は(p)は以下の95%信頼区間に収まる。 つまり95%信頼区間は以下となる。 begin{eqnarray} hat{p}-1.96 sqrt{frac{hat{p}(1-hat{p})}{n}} &le& p le hat{p}+1.96 sqrt{frac{hat{p}(1-hat{p})}{n}} \\ 0.9-1.96 frac{sqrt{0.09}}{sqrt{139}} &le& p le 0.9 + 1.96 frac{sqrt{0.09}}{sqrt{139}} \\ 0.9-1.96 frac{0.3}{11.78} &le& p le 0.9+1.96 frac{0.3}{11.78} \\ 0.85 &le& p le 0.95 end{eqnarray} (n)を下げたい場合は有意水準を下げれば良い。 統計的に有意水準=10%まで許容されることがある。 有意水準が10%であるとすると、(n)は35以上であれば良いことになる。 begin{eqnarray} 1.96sqrt{frac{0.9(1-0.9)}{n}} le 0.1 \\ frac{1.96}{0.1} sqrt{0.09} le sqrt{n} \\ 5.88 le sqrt{n} \\ 34.6 le n end{eqnarray} 信頼区間と有意水準の式において(p)を標本から取ってきたけど、 アンケートにおいてYes/Noを答える場合、(p)は標本における最大値（つまり0.5)を 設定して(n)を求める。 つまり、(p)として利用するのは標本比率ではないのかな？と。 このあたり、(hat{p})を変数として残すとどういうことがわかった。 [clink url=\"https://ikuty.com/2019/01/13/sampling_with2/\"]

[mathjax] zozoスーツではないけれども、標本が正規分布に従うというのは、 真実の母平均に対して正規分布に従う計測誤差を含む分布を観測しているのと同じ。 母平均(mu)が未知である現象を計測誤差がわかっている計測手段で計測する話。 (n)回計測を行って得られた標本(X_1,X_2,cdots,X_n)は、母平均を中心として誤差分振れているはず。 つまり、(X_1=mu+e_1,X_2=mu+e_2,cdots,X_n=mu+e_N)。 誤差(e_i)は平均0、分散(sigma^2)の正規分布(N(0,sigma^2))に従うと考えると、 標本(X_i)は(mu)だけオフセットした正規分布(N(mu,sigma^2))に従うと考えられる。 標本平均は(bar{X}=frac{1}{n}(X_1+X_2+cdots+X_n))だから、 (n)に関係なく(E(bar{X})=mu)であって、(V(bar{X})=frac{sigma^2}{n})から(lim_{nrightarrow infty}V(bar{X})=0)。 中心極限定理により大なる(n)のとき(bar{X})の分布は正規分布(N(mu,frac{sigma^2}{n}))で近似できる。 標本(X_i)の標準偏差が(sigma)である一方、標本平均の標準偏差は(frac{sigma}{sqrt{N}})だから、 標本の分布より、標本平均の分布の方が裾が狭い。 正規分布(N(mu,frac{sigma^2}{n}))を標準化しておくと、標準正規分布の累積度数表を使って 平均(mu)、標準偏差(sigma)を評価できるようになる。z得点は以下の変換。 begin{eqnarray} Z=frac{bar{X}-mu}{frac{sigma}{sqrt{n}}} end{eqnarray} 分布(Z)は平均0、標準偏差1の標準正規分布になる。 見方としては、残差が標準偏差何個分か？の分布。全部足して1になる。 (bar{X},mu,sigma,n)として具体的な値を入れると数値(Z)が決まる。 ちなみに確率密度関数と累積度数は以下の通り。 begin{eqnarray} f(x) &=& frac{1}{sqrt{2pi}} exp left( -frac{x^2}{2} right) \\ int_{-infty}^{infty} f(x) dx &=& 1 end{eqnarray} (x=0)から(x=z)の面積(int_0^{z} frac{1}{sqrt{2pi}} left( -frac{x^2}{2} right) )を(Phi(z))とおき、 (Phi(z)=a)となる点を上側(a)パーセント点という名前が付いている。 (Phi(z))の積分は解析的に計算できないけれど、有用だし決まった数値なので、 ここみたいに表ができているからルックアップすれば良い様子。 (Z)得点が1.96であったとすると、標準正規分布表から(Phi(z=1.96)=0.475)であることがわかる。 これは上側確率が0.475という意味なので、両側確率は2をかけて0.975ということになる。 逆に言うと、(mu)だけが不明で、既知の母分散と標本平均から(mu)を推測することに、 この話を使うことができる。つまり、(-1.96 le z le +1.96)という式を立てると、 (mu)の信頼区間を作ることができる。つまり、(n)個の標本を取る操作を100回繰り返すと97.5回は 信頼区間が母平均を含まない区間になっている。 例 確率変数(X)が平均2、分散10の正規分布(N(2,10))に従うとする。 95%信頼区間は(-1.96 lt z lt 1.96)から、 (-1.96 sqrt{10} + 2 lt X lt 1.96 sqrt{10} + 2)。 (-4.2 lt X lt 8.2)。 100回試すと97.5回は母平均がこの区間にある。 (X)が負になる確率は、(Z=frac{X-2}{sqrt{10}})から、(sqrt{10}Z+2lt 0)、(Z lt -frac{2}{sqrt{10}})、(Z lt - 0.633)。 (P(X lt 0)=P(Z lt -0.633)=1-P(z lt 0.633))。

[mathjax] 正規母集団からの推定をやる前に、正規分布の再生性の理解が必要だったのでまとめてみる。 独立な確率変数(X_1)、(X_2)がそれぞれ確率分布(g(x))、(h(x))に従うとする。 各確率変数の和(X_1+X_2)が従う確率分布を(k(z))とする。 確率(P(X_1+X_2=z))を考えると、(X_1+X_2=z)となるのは、 (X_1=x, X_2=z-x)としたとき、両者を足して(z)になる全ての組み合わせ。 (X_1)は(g(x))、(X_2)は(h(z-x))に従うので、両者が同時に起こる確率は(g(x)h(z-x))。 これをまとめて書くと、 begin{eqnarray} k(z) = sum_x g(x)h(z-x) end{eqnarray} この形が「たたみこみ(convolution)」。　 (k = g * x)と書く。 確率変数(X_1)、(X_2)が独立で、それぞれ平均(mu_1)、(mu_2)、分散(sigma_1^2)、(sigma_2^2)の正規分布に従うなら、 以下が成り立つ。 begin{eqnarray} N(mu_1,sigma_1^2) * N(mu_2,sigma_2^2) = N (mu_1+mu_2, sigma_1^2 + sigma_2^2) end{eqnarray} これ、モーメント母関数を使って証明できる様子。 ある分布のモーメント母関数があったとして、モーメント母関数を(n)回微分して変数を(0)と置くと、 分布の期待値、分散、歪度、突度など統計量を求められるやつ。 [clink url=\"https://ikuty.com/2018/09/22/moment_generating_funuction/\"] 正規分布の確率密度関数とモーメント母関数は以下の通り。 begin{eqnarray} f(x) &=& frac{1}{sqrt{2pisigma}} expleft( - frac{(x-mu)^2}{2sigma^2} right) \\ M(t) &=& exp left( mu t + frac{sigma^2 t^2}{2} right) end{eqnarray} もちろん、(N(mu_1,sigma_1^2))、(N(mu_2,sigma_2^2))のモーメント母関数は, begin{eqnarray} M_1(t) &=& exp left( mu_1 t + frac{sigma_1^2 t^2}{2} right) \\ M_2(t) &=& exp left( mu_2 t + frac{sigma_2^2 t^2}{2} right) \\ end{eqnarray} かけると、以下の通り(N(mu_1+mu_2,sigma_1^2+sigma_2^2))のモーメント母関数となる。 begin{eqnarray} M_1(t) M_2(t) &=& expleft( mu_1 t +frac{sigma_1^2 t^2}{2} right) expleft( mu_2 t + frac{sigma_2^2 t^2}{2} right) \\ &=& expleft( (mu_1+mu_2) t +frac{(sigma_1^2 + sigma_2^2) t^2}{2} right) end{eqnarray} たたみこみの操作は、独立な確率変数(X_1,X_2)について(X_1+X_2)の確率分布を求める操作だから、 この結果は独立な確率変数(X_1,X_2)が(N(mu_1,sigma_1^2))、(N(mu_2,sigma_2^2))に従うとき、 (X_1+X_2)が(N(mu_1+mu_2,sigma_1^2+sigma_2^2))に従うことを意味する。 ある確率分布のたたみ込みの結果が同じ確率分布になることを再生性(reproductive)というらしい。 正規分布の再生性を使った演算 正規分布には再生性があるので、以下みたいな演算ができる。 (X_1,X_2,cdots,X_n)が独立で、それぞれ正規分布(N(mu_1,sigma_1^2),N(mu_2,sigma_2^2),cdots,N(mu_N,sigma_N^2) )に 従うとき、(X_1+X_2+cdots+X_n)は(N(mu_1+mu2+cdots,mu_N,sigma_1^2+sigma_2^2+cdots+sigma_N^2))に従う。 (X_1,X_2,cdots,X_n)が全て同じ(N(mu,sigma^2))に従うなら、 (X_1+X_2+cdots+X_n)は、(N(nmu, nsigma^2))に従う。 (bar{X}=frac{X_1+X_2+cdots+X_n}{n})は(N(mu,frac{sigma^2}{n}))に従う。

[mathjax] 不偏分散は(frac{1}{n} sum_{i=1}^n (X_i-bar{X})^2)ではなく、(frac{1}{n-1} sum_{i=1}^n (X_i-bar{X})^2)。 分母から1を引く必要がある。なんでか調べてみたので書いてみる。 標本平均は(n)の大小によらず母平均の近傍にあって、母平均に確率収束する。 標本平均は(n)の大小に関係なく、その期待値と母平均が等しい（不偏）。 begin{eqnarray} E(bar{X}) &=& frac{1}{n}nmu = mu \\ lim_{n rightarrow infty} V(bar{X}) &=& 0 end{eqnarray} 100個のデータがあって、その中から5個取ったときの平均と、50個取ったときの平均に 母平均の推測という意味で違いがない。 では、分散はどうか。 定義通り標本の分散を(S^2 = frac{1}{n}{ (X_1-bar{X})^2 + (X_2-bar{X})^2 + cdots + (X_n-bar{X})^2 } )とすると、 (S^2)は母分散と等しくならない。不偏にならない。つまり、(E(S^2) ne sigma^2)。 その値が不偏であるか否かは、実際に期待値を式変形してみるとわかる。 結論を知っていないと出来ない変形ばかりだけども...。 begin{eqnarray} E(S^2) &=& Eleft[frac{1}{n} sum_{i=1}^n (x_i-bar{X})^2 right] \\ &=& frac{1}{n} Eleft[ sum_{i=1}^n (x_i-bar{X})^2 right] \\ &=& frac{1}{n} Eleft[ sum_{i=1}^n left( (x_i-mu)-(bar{X}-mu) right)^2 right] \\ &=& frac{1}{n} Eleft[ sum_{i=1}^n (x_i-bar{X})^2 -2sum_{i=1}^n(x_i-mu)(bar{X}-mu) + sum_{i=1}^n (bar{X}-mu)^2 right] \\ &=& frac{1}{n} Eleft[ sum_{i=1}^n (x_i-bar{X})^2 -2n (bar{X}-mu) +n (bar{X}-mu)^2 right] \\ &=& frac{1}{n} Eleft[ sum_{i=1}^n (x_i-bar{X})^2 - n(bar{X}-mu)^2 right] \\ &=& frac{1}{n} sum_{i=1}^n Eleft[ (x_i-mu)^2 right] - Eleft[ (bar{X}-mu)^2 right] \\ &=& frac{1}{n} sum_{i=1}^n V(x_i) - V(bar{X}) \\ &=& sigma^2 - frac{1}{n} sigma^2 \\ &=& frac{n-1}{n} sigma^2 end{eqnarray} ということで、(E(S^2)ne sigma^2)。不偏でない。 では、どうすれば不偏な標本分散を得られるのか。 (E(S^2)=frac{n-1}{n} sigma^2)から、(frac{n}{n-1}E(S^2)=sigma^2)なので、(s^2=frac{n}{n-1}E(S^2))とすれば、 (s^2=sigma^2)ということになり、(s^2)は不偏となる。(s^2)を不偏分散という。 begin{eqnarray} s^2 = frac{n}{n-1} { (X_1-bar{X})^2 + (X_2-bar{X})^2 + cdots + (X_n-bar{X})^2 } end{eqnarray} 100個データがあって、10個データをとったときと、100個データをとったときの (E(S^2))の母分散とのズレは以下の通り。10個のとき(E(S^2))をそのまま計算してしまうと、 その値は母分散から10%もズレてしまう。100個にしても1%ずれる。 begin{eqnarray} E(S_{10}^2) &=& frac{9}{10}sigma^2 \\ E(S_{100}^2) &=& frac{99}{100}sigma^2 \\ end{eqnarray}

[mathjax] 統計のど真ん中。大数の法則と中心極限定理。 確かに奇跡的な美しさを感じる...。 同じ確率分布に従う独立な確率変数(X_1,X_2,cdots,X_n)について、(n)が大きければ (bar{X}=frac{X_1+X_2+cdots+X_n}{n})と置いたときに、(E(bar{X})=E(X_n))として扱ってよかった。 （(bar{X})は(mu)に確率収束した。(n)が次第に大きくなるにつれて(V(bar{X)})が0に収束した。) begin{eqnarray} lim_{n rightarrow infty}P(|X_n-mu|ge epsilon) rightarrow 0 end{eqnarray} なんとなく、最頻値を峰として(n)の増加に伴って峰が険しくなっていきそうだけど、 実際、(X_i)がどんな確率分布に従っていたとしても、 (S_n=X_1+X_2+cdots+X_n)は平均(nmu)、標準偏差(sqrt{n}sigma)の正規分布に従う（と考えて良い）らしい。 (bar{X})は平均(mu)、標準偏差(frac{sigma}{sqrt{n}})の正規分布に従う。 これが中心極限定理(Central limit theorem)。 (n)が大きくなったとき、ランダムウォーク(S_n)は平均(nmu)、標準偏差(sqrt{n}sigma)の分布が正規分布になる、 というのは以下のようにかけるらしい。 右辺は標準正規分布の確率密度関数の定積分。 左辺は既に正規分布であることを前提にしているような... begin{eqnarray} Pleft(a le frac{S_n -nmu}{sqrt{n}sigma} le bright) rightarrow int_a^b frac{1}{sqrt{2pi}}e^{-frac{x^2}{2}}dx end{eqnarray} (frac{S_n-nmu}{sqrt{n}sigma}=frac{n(frac{1}{n}S_n-mu)}{sqrt{n}sigma}=frac{bar{X}-mu}{frac{sigma}{sqrt{n}}}) みたいに変形すると、 以下のようにすることもできる。 正規分布を標準正規分布に変換するために（つまり平均=0、標準偏差=1にするために)、 標準偏差で割る（つまり標準偏差何個分か？）変換。 begin{eqnarray} Pleft(a le frac{bar{X}-mu}{frac{sigma}{sqrt{n}}} le b right) rightarrow int_a^b frac{1}{sqrt{2pi}}e^{-frac{x^2}{2}}dx end{eqnarray} 統計学入門に奇跡的な証明が書いてあるけど完全に写経になるので終わり。

[mathjax] 大数の法則は、標本数が大きくなったときに標本平均が母平均に確率収束することを数学的に表す。 begin{eqnarray} lim_{n rightarrow infty} P(|bar{X}-mu|ge epsilon) rightarrow 0 end{eqnarray} 母集団の確率分布がわかっていれば、母平均、母分散、標準偏差などの統計量は定義済みだから、 標本数が多いという仮定があれば、標本から母平均を見立てて計算し、 見立てた母平均が\"どの程度\"真の母平均に収束しているかを言うことで、 標本がどの程度母集団を言い当てているかを言うことができる。 1か0が起こるベルヌーイ試行を100回やって、1が80回、0が20回起こったとする。 この100回の標本は\"たまたま\"取得した結果だけれども、 これが無限集団の一部分だったとして、今回得られた標本から無限集団の平均を得られないか。 まず、確率変数(X_1,X_2,cdots,X_{100})は以下の確率分布に従う。 begin{eqnarray} P(X_i=1) &=& frac{80}{100} = 0.8 \\ P(X_i=0) &=& frac{20}{100} = 0.2 end{eqnarray} さらに確率変数(X_1,X_2,cdots,X_{100})は(p=0.8)のベルヌーイ分布に従う。 その平均は(p=0.8)、分散は(p(1-p)=0.16)。 これは\"たまたま\"取得した100個の標本を使って作った母平均と母分散。 では、この\"たまたま\"の母平均と母分散を使って、 (bar{X}=frac{X_1+X_2+cdots+X_n}{n})に関する(E(bar{X}))、(V(bar{X}))を計算する。 (E(bar{X})=mu)、(V(bar{X})=frac{sigma^2}{n})であるから、 (E(bar{X})=0.8)、(V(bar{X})=0.0016)、（(sqrt{V(bar{X})}=0.04 )) 大数の法則（弱法則）により、\"たまたま\"の100回という試行において、 (bar{X})は0.8に確率収束し、その標準偏差は(0.04)であることがわかる。 試行回数を増やしていくと、標準偏差が0に収束する。 100回という試行回数において、1が起こる回数は平均80回で、上下に(pm 4% )ばらける。 無限大まで増やしていく途中、100回という区切りではこうだった、ということで、 さらに増やしていくと、上下の幅が(pm 4% )からさらに狭くなる。 全ての標本を得ないでも、手元の標本だけから、 上下の幅という曖昧さを残すことで推測できる。 大数の法則 (n)個の確率変数(X_1,X_2,cdots,X_n)が平均(mu)、分散(sigma^2)の分布に独立に従うとき、 和を(n)で割った平均の分布(bar{X}=frac{X_1+X_2+cdots+X_n}{n})について、平均、分散はそれぞれ、 begin{eqnarray} E(bar{X}) &=& frac{nmu}{n} = mu \\ V(bar{X}) &=& frac{nsigma^2}{n^2} = frac{sigma^2}{n} end{eqnarray} (n)の極限をとったとき、標本平均は母平均(mu)に収束し、分散はゼロに収束する。 begin{eqnarray} lim_{n rightarrow infty} E(bar{X}) &=& mu \\ lim_{n rightarrow infty} V(bar{X}) &=& 0 \\ end{eqnarray} (n)が少ない間は、(mu)の近傍に散らばって分布するが、 (n)が増えるに従って、(mu)の近傍の散らばりの範囲が狭くなっていく。 昔、写経した証明は以下。 だいぶ普通に読めるようになった。 [clink url=\"https://ikuty.com/2018/07/06/law_of_large_numbers/\"] 上で書いたように(V(bar{X})=frac{sigma^2}{n})だから、(sigma_n^2)は、(sigma_n^2=V(bar{X}_n)=frac{sigma^2}{n}) 以下のチェビシェフの不等式において、 begin{eqnarray} P(|X_n−mu|ge ksigma_n) &le& frac{1}{k^2} end{eqnarray} (ksigma_n=epsilon)とおく。(k=frac{epsilon}{sigma_n}=sqrt{n} frac{epsilon}{sigma}) だから、 begin{eqnarray} P(|X_n−mu|ge epsilon) &le& frac{sigma^2}{nepsilon^2} end{eqnarray} 極限を取ると、 begin{eqnarray} lim_{n rightarrow infty} P(|X_n−mu|≥epsilon) &=& lim_{n rightarrow infty }frac{1}{epsilon^2} cdot frac{sigma^2}{n} \\ &rightarrow& 0 end{eqnarray}

[mathjax] 620kgまで積載できるエレベータ。定員は何人でしょうか..? ネガティブ思考の大家としては、 太っている人が乗ることを考えて少なめな方が...と常日頃思っていたのだけど、 案外、悪いことばかり考えるのも客観的でないことがわかる。 数を打つと、大きいもとと小さいものが打ち消しあう。 良いことと悪いことが打ち消しあって平均に収まる。 確率変数が大量((n)個)にあったとき、それらが全て同じ分布に独立に従うとして、 その分布の平均が(mu)、分散が(sigma^2)だとする。それらの平均、分散は、 begin{eqnarray} E(X_1) = E(X_2) = cdots = E(X_n) = mu \\ V(X_1) = V(X_2) = cdots = V(X_n) = sigma^2 \\ end{eqnarray} 和の平均、分散は、標準偏差は、 begin{eqnarray} E(X_1+X_2+cdots + X_n) = nmu \\ V(X_1+V_2+cdots + X_n) = nsigma^2 \\ s = sqrt{n} sigma end{eqnarray} 人間1人の体重が平均65.0kg、標準偏差が6.0kgだとすると、 9人の平均は人数に比例して580kg。人数の増加にしたがって順当に見積もりが増える。 どれぐらい上にマージンが必要かというと..標準偏差は(sqrt{9}*6=18.0)kg。 最大積載重量が600kgであったとしても(+1sigma)の範囲に収まる！ 620kgあれば、(+2sigma)も収まる。9人はいける。 振れ幅は(sqrt{n})に比例する。 (n)が大きくなっても大して振れ幅は大きくならない。 ちなみに、法令上（建築基準法施行令第129条）は、 1人あたり65kgとして頭数で割るだけが決まっているそうで、 体重の振れ幅については結局謎のまま。 床面積あたりにするとキツキツになる模様。 http://www.mitsubishielectric.co.jp/elevator/inquiry/faq/category01/qa3.html 最大積載重量が実際のハードリミットでそこでブザーがなるんだから、 定員、というソフトリミットは目安でしかないんだな。