日々の学び

dbt Chartsで自作の財務分析基盤を可視化してみた話

dbt chartsが利用可能になった、ということなので、実際に使ってみた。 データソースは自作の財務分析基盤。EDINET報告書/JPX適時開示・日足終値を収集・加工している。 ClaudeProの$20プランでSonnet5.5を使った。自作の基盤はDuckDB/Prefect/dbtで作っていて、 同じdbtファミリーとなるが、Claudeはdbtの上に来るため正直dbtファミリーとしてdbt-chartsが あるメリットをあまり感じなかった。Claudeがmanifest相当の内容を理解してしまうし、 DuckDBにアクセスしてデータを読むことで、ymlのメタデータを超える理解をしてしまう。 むしろ、メタデータに嘘が書かれていてClaudeを惑わしてしまい、途中基盤側を修正したりした。 財務や適時開示といったメジャーなオープンデータなので指標の集計はLLMの知識で解決できたが、 これが企業の中のクローズドなデータだったら、説明する層が必要になると思う。 dbt-charts は Vega-Lite という描画ライブラリをベースとしている。 一般的なチャートは用意されていたが、カスタマイズ性は並かな...と思った。 あまり思った通りの表示にならない。 Vegaという、描画要素を全てjsonで書く低レイヤーの描画ライブラリの上に、 より高水準の Vega-Lite があり、その過程である程度機能が落ちているように見える。 「少ない記述で、よく使う種類のグラフを宣言的に書く」という趣旨らしい。 dbt chartsの構成ファイルとしてymlを書いて、それが宣言的な記述になって、 Vega-Liteにマップされる、という感じ。 [clink implicit=\"false\" url=\"https://github.com/ikuty/finance-dwh\" imgurl=\"https://ikuty.com/wp-content/uploads/2026/10/github.jpg\" title=\"finance-dwh\" excerpt=\"個人運用の財務データ分析基盤におけるウェアハウス層。金融庁のEDINET API v2 (https://api.edinet-fsa.go.jp/api/v2/)から入手した有価証券報告書・四半期報告書 等の開示書類データをデータソースとするレイク層(finance-lake、生データのまま保存) を入力とし、raw / landing / cleansed / intermediate / mart の層からなる DWH を構築する。 将来的には、AI(エージェント・LLM経由の問い合わせ)・BI(ダッシュボード等)双方から指標定義の 一貫性を保てるよう、セマンティックレイヤーの構築を目指している(未着手、mart層と 同様「利用目的が固まってから」着手する方針)。\"] [arst_toc tag=\"h4\"] type=lineで終値比較 まずは終値比較。 NTT, KDDI, Softbank, 楽天 について、1年前の終値を100とし、終値を重ねてみてください。 type=barで通期のバリュエーション推移 dbt charts で、通信4社(NTT・KDDI・ソフトバンク・楽天)の通期(2021〜2025年)の 主要指標を、実データで横棒グラフにしてください(PNG出力)。 - 構成: 縦軸=決算期、指標ごとのパネル(1段目: 時価総額・PER・PBR、2段目: ROE・ROA・EPS)、 4社のグループ棒、棒に数値ラベル、順序は NTT → KDDI → ソフトバンク → 楽天 - 時価総額・PER・PBR: 開示日(submit_date_time の日付)の翌日以降で最初の取引日の値 - EPS: 開示値 × period_end_cum_adj(分割調整後) - ROE・ROA: % 表示(ROA は net_income / total_assets) type=areaで時価総額積み上げ比較 dbt charts で、通信4社(NTT・KDDI・ソフトバンク・楽天)の時価総額の日次推移を 積み上げ面グラフにしてください(2021年〜、実データ、PNG出力)。 - データ: mart__jpx_edinet__daily_valuation_indicators の market_cap(兆円) - 銘柄: jpx_code 9432 / 9433 / 9434 / 4755 - 積み順: 下から NTT → KDDI → ソフトバンク → 楽天 - 横軸: 年ラベル(2022〜2026)を横書きで毎年表示 type=scatterで全上場企業のPBRxROE分布 株価純資産倍率(PBR)、つまり、時価総額 ÷ 純資産で、市場が会社の純資産の何倍で評価しているか。 自己資本利益率(PER)、つまり、純利益 ÷ 純資産で、株主のお金でどれだけ稼いだか。 両者は PBR = PER x ROE で繋がっている。PERが同程度ならROEが高いほどPBRも高くなるはず。 会社の資産価値(解散価値)に対して株価が割安になっている、 つまり株価が1株当たり純資産を下回っている状態の会社が多数あることで有名。 PBRを1倍より大きくすることが要請されている。 PBRを対数軸とすると、PBR1倍割れの銘柄のPERの分布がわかる。 基本的にはPBR1倍割れの銘柄はROEも低いが、左上にROEが高いものもある。 雪印メグミルク(2270)、大和ハウス工業(1925)、石原産業(4028)、サンフロンティア不動産(8934)などが 低PBRだが比較的ROEが高い。北陸電力(9505)、西日本鉄道(9031)も同じグルーピングだが、 自己資本比率が低くROEが高めに出やすい。 dbt charts で、全上場会社の PBR × ROE の散布図を、実データで1枚の画像(PNG)に 上下2パネルで出してください。上: PBR 線形軸、下: PBR 対数軸。 - 1社1点。PBR は daily_valuation_indicators の最新日、ROE は各社の最新の通期の値 - 点の色: 会計基準(日本基準=淡い灰青、IFRS=オレンジ、US GAAP=濃い青)。多数派を先に描く - 点の大きさ: 時価総額(億円の常用対数 + 1) - 外れ値除外: PBR 0以下・10超、ROE -50〜60% の外 - ROE は % 表示、点の透明度 0.6、タイトルは大文字を変換しない type=pie,donutでPBRレンジ構成比,会計基準ごとの時価総額 円グラフとドーナツ。ドーナツは円グラフが流用されている。 dbt charts で、全上場会社の構成比を、ドーナツと円グラフの2枚を左右に並べた1枚の画像(PNG)にしてください(実データ)。 - 左(ドーナツ): 最新日の PBR レンジ別の会社数(PBR 0 超)。区分は 0.5未満 / 0.5〜1倍 / 1〜2倍 / 2〜3倍 / 3倍以上、 各区分に「社数」ラベル、中央に合計会社数。1倍割れは橙系、1倍以上は青系の色。 - 右(円グラフ): 会計基準別の時価総額(兆円)。会計基準は各社の最新の通期の値。 色は 日本基準=淡い灰青、IFRS=オレンジ、US GAAP=濃い青。各スライスに「基準名 + 兆円」ラベル。 type=kpiで市場サマリー 自力でラベルを並べなくても、1つのチャートで値を表示できる。 dbt charts で、全上場会社の市場サマリーを KPI 3枚を1行に並べたボードにして、 実データで PNG 出力してください。 - 対象日: daily_valuation_indicators の最新日。補足に「前年同日比」(1年前以前で直近の取引日との差) - KPI 1: PBR 1倍割れの会社の割合(PBR 0 超の会社のうち PBR 1 未満、% 表示、差はポイント) - KPI 2: PBR 中央値(PBR 0 超、小数2桁) - KPI 3: PER 中央値(PER 0 超〜100 以下、小数1桁) type=tableで時価総額1000億円以上のPBR1倍割れROE10%超え企業一覧 これは説明不要。何も考えずに多用しそう。 - 並び: ROE の降順、全行を表示(ページ分割しない) - 列: コード / 会社名 / 会計基準 / PBR(倍)/ ROE(%)/ PER(倍)/ 時価総額(億円)/ 自己資本比率(%) type=heatmapでPBRレンジxROEレンジ会社数分布 dbt charts で、全上場会社の PBR × ROE の会社数ヒートマップを、実データで PNG 出力してください。 - 縦軸: 最新日の PBR(上が高い)。12 階級: 0.25未満 / 0.25-0.5 / 0.5-0.75 / 0.75-1 / 1-1.25 / 1.25-1.5 / 1.5-2 / 2-2.5 / 2.5-3 / 3-4 / 4-5 / 5以上(PBR 0 超の会社) - 横軸: 各社の最新の通期 ROE(左が低い)。12 階級: -10%未満 / -10〜0% / 0〜2.5% / 2.5〜5% / 5〜7.5% / 7.5〜10% / 10〜12.5% / 12.5〜15% / 15〜20% / 20〜25% / 25〜30% / 30%以上 - マスの値: 会社数(色が濃いほど多い)。配色は橙〜赤(dbt-seq-rust) type=pointmapで本店所在地都道府県を日本地図にマップする むちゃくちゃ苦戦した。日本地図は外部から持ってこないと表示されない。 しかも、webサーバでホストしないといけない。 もっと工夫すると出来るのかもしれないが、少なくとも公式に説明はない。 dbt charts で、本店所在地の都道府県別の時価総額合計を、日本地図の上に点の大きさで表した PNG を、実データで出してください(type: point_map)。 - 集計: 最新日の market_cap を、各社の最新の通期の prefecture(本店所在地)ごとに合計(兆円) - 点の位置: 47 都道府県庁所在地の緯度・経度(SQL に埋め込み)。点の大きさ = 時価総額合計 - 背景: 日本の都道府県境界(amay077/JapanPrefGeoJson 由来を簡略化したもの)。灰色の塗り、白い枠線 - 描画は scripts/dct_render.py 経由で行い、境界データは手元の HTTP サーバ(ポート 8765)で配信する まとめ 自作の財務分析基盤をdbt chartsで可視化してみた。 JavaScriptを書けない環境でもdbtファミリーで可視化まで完結させられることを確認した。 ymlで宣言的に図を書けるのでシンプルな構成に良い。

Snowflake SAML2.0 Security Integrationを使用したSP/IdP Initiated SAML Federationと構成の詳細

AWS IAM Identity Center を IdP として使用し、Snowflake SAML2.0 Security Integration を構成する。 構成の際に多くのパラメタの設定が必要だが、設定可能なパラメタの意味について深掘りしてみる。 この記事は自分の学び用なので、事実の確認、説明用画像作成のために生成AIを使用するが、 記事の作成、校正には使用しない。 [arst_toc tag=\"h4\"] SP起点(SP initiated) flow SP側にSSOボタンなどを配置して、SSOボタン押下でSSO認証とSPログインを開始するフロー。 SP,User,IdP間の呼び出しシーケンスは下図の通り。 IdP起点(IdP initiated) flow IdP側にログインボタンを配置して、ログインボタン押下でIdP認証とSPログインを開始するフロー。 SP,User,IdP間の呼び出しシーケンスは下図の通り。 ログアウト SP側のセッションと、IdP側のセッションは独立している。SP起点、IdP起点のいずれにおいても、 基本的には、片方をログアウトしたからといってもう片方が勝手にログアウトしたりしない。 ChromeでSP起点でフェデレーションログインした後、ChromeでSPのセッションをログアウトした場合、 IdP側のセッションはまだ生きているので、Chromeで再度フェデレーションを開始したとき、 IdP側の認証は走らず、SPにログインできる。 ChromeでSP起点でフェデレーションした後、SafariでSP起点でフェデレーションしたとき、 Chromeでログアウトしたとしても、Safariのセッションはログアウトしないため操作を続行できる。 ただし、EntraIDのみ、IdP起点で開始したとき、「グローバルログアウト,シングルログアウト,SLO」が サポートされており、IdPからログアウトすると、全てのセッションからログアウトする。 たしかに、EntraIDが気持ち悪い動作するな、と言う時、これが動いている時がありそう。 セッションタイムアウト SP(Snowflake)のセッションがタイムアウトした場合、ユーザはIdPを介して再度認証が必要。 IdPでCancel操作をすると、そこでセッションを終了できる。 IdPのセッションがタイムアウトした場合、Snowflakeセッションに影響しない。 その時点でアクティブなSP(Snowflake)セッションは生きたままとなる。 識別子優先ログイン 組織ごとにそれぞれのIdPとSAML連携したいといった場合がある。複数のintegrationを持てる。 また、ユーザによってはフェデレーション連携させずに、Snowflake認証だけにしたいケースもある。 ユーザによって認証に必要な入力が異なるため、全てのユーザに対して最小公倍数的に入力項目を 出してしまうと、ユーザによって不要な項目が並んでいるように見える。 識別子優先ログインをONにすると、SP(Snowflake)側の認証入力が多段階となる。 つまり、1段階目で識別子(ユーザ名、またはメールアドレス) を入力させ、 入力された識別子がどの認証方式に紐づいているかを判定したのち、 ユーザに適した2段階目の認証入力画面(PW入力、SSOボタン)を表示する。 画面遷移が増えるが、不要な入力項目が現れなくなる。 SAML2.0 Security Integration これを作るだけでSnowflakeにSAML2.0 Federationを追加できる。 CREATE [ OR REPLACE ] SECURITY INTEGRATION [ IF NOT EXISTS ] <name> TYPE = SAML2 ENABLED = { TRUE | FALSE } { METADATA_URL = \'<string_literal>\' | <idp_parameters> } [ ALLOWED_USER_DOMAINS = ( \'<string_literal>\' [ , \'<string_literal>\' , ... ] ) ] [ ALLOWED_EMAIL_PATTERNS = ( \'<string_literal>\' [ , \'<string_literal>\' , ... ] ) ] [ SAML2_SP_INITIATED_LOGIN_PAGE_LABEL = \'<string_literal>\' ] [ SAML2_ENABLE_SP_INITIATED = TRUE | FALSE ] [ SAML2_SNOWFLAKE_X509_CERT = \'<string_literal>\' ] [ SAML2_SIGN_REQUEST = TRUE | FALSE ] [ SAML2_REQUESTED_NAMEID_FORMAT = \'<string_literal>\' ] [ SAML2_POST_LOGOUT_REDIRECT_URL = \'<string_literal>\' ] [ SAML2_FORCE_AUTHN = TRUE | FALSE ] [ SAML2_SNOWFLAKE_ISSUER_URL = \'<string_literal>\' ] [ SAML2_SNOWFLAKE_ACS_URL = \'<string_literal>\' ] [ COMMENT = \'<string_literal>\' ] パラメタの意味を1つずつ解釈していく。 METADATA_URLとidp_parameters まず、METADATA_URL について、IdPが証明書の更新などのIdP構成設定を動的に 取得・同期できるように、URLを公開している場合がある。 その場合、METADATA_URL に指定することができる。ない場合は idp_parameters を手動で設定する。 idp_parameters として設定できる項目は以下の通り。 SAML2_ISSUER IdPメタデータの entityID(IdP Entity ID)。IdP(AWS IAM Identity Center)を一意に識別する文字列。Snowflakeはこの値をIdPからの応答(SAML Response)の Issuer と照合し、正規のIdPからの応答かを検証する。 SAML2_SSO_URL IdPメタデータの SingleSignOnService Location。SP-initiatedログイン時、SnowflakeがブラウザをリダイレクトさせるIdPのログインエンドポイント。 SAML2_PROVIDER IdPの種別。OKTA、ADFS、Customのいずれか。専用値がないIdP(AWS IAM Identity Centerを含む)はすべて Custom を指定する。 SAML2_X509_CERT IdPメタデータの X509Certificate。IdPがSAML Responseに付与する署名を検証するための公開鍵証明書。ヘッダー(-----BEGIN CERTIFICATE-----)・フッター・改行を除いたBase64本文のみを指定する。 多層防御 ALLOWED_USER_DOMAINS SAML2セキュリティ統合で認証できるメールドメインの許可リスト。IdPから返ってきた識別子(login_name/emailとして扱われる値)の @ 以降がこのリストに含まれていないと、SAMLレスポンス自体は正当でも認証を拒否する。IdP設定ミスや、共有IdPで意図しないドメインのユーザーがログインできてしまう事故を防ぐための追加防御層。ただし、NameIDにメールアドレスではなくusernameを使う場合、@domain部分が存在せず、このパラメータで意味のある制限をかけることができない。メールアドレス形式のNameIDに切り替える場合は設定を検討する。 ALLOWED_EMAIL_PATTERNS SAML2セキュリティ統合で認証するメールアドレスが一致すべき正規表現のリスト。ALLOWED_USER_DOMAINS のドメイン単位の許可リストより柔軟な条件(例: .*@example.com$ や特定のローカルパートのみ許可など)を指定したい場合に使う。こちらも同様に、NameID(メールアドレス形式ではない)には適用できない。 SP_INITIATED flowの設定 SAML2_SP_INITIATED_LOGIN_PAGE_LABEL Snowflakeのログイン画面の「Log In With」ボタンの後に表示するラベル文字列。任意の名称でよい。 SAML2_ENABLE_SP_INITIATED ログインページに「Log In With」ボタンを表示するかどうか。TRUEでSP-initiatedログインが有効になる。FALSEの場合はIdP-initiatedのみ許可。 SAMLメッセージの署名 SAML2_SIGN_REQUEST SnowflakeがIdPに送るSAMLリクエスト(AuthnRequest)に署名するかどうか。署名させたい場合は TRUE にし、AWS側アプリ設定でリクエスト署名検証を有効化する必要がある。TRUEにするメリット: IdPからSPへの応答(SAMLレスポンス)はSAML2_X509_CERTによる署名検証で既に真正性が保証されているが、逆方向(SPからIdPへのAuthnRequest)は署名なしだと改ざん・なりすましを検知できない。TRUEにすることでIdP側がリクエストの送信元(Snowflake)と内容の完全性を検証できるようになり、偽装されたAuthnRequestでユーザーを不正なフローに誘導される攻撃への耐性が上がる。特にセキュリティ要件が厳しい環境や監査対応が必要な場合はTRUEを推奨。 SAML2_SNOWFLAKE_X509_CERT 指定しない場合Snowflakeが自己署名した証明書が使われるが、SAML2_SNOWFLAKE_X509_CERT により、第三者(CA)による証明書を設定できる。これにより、SPからIdPへのAuthnRequestを第三者(CA)による証明書(の秘密鍵)で署名するようになる。 SPからIdPへのAuthnRequestを署名する限界と価値 そもそも「証明書」とは「公開鍵」+「(誰かの)秘密鍵による発行者の署名」という構造となっている。第三者が署名するといっても、その第三者が何者か、についてはより上位の第三者が署名することで証明される。最終的には、ルートCAと呼ばれるトラストストアに登録されている事前登録済みの発行者に突き当たり、ルートCAは自己署名を行っている。ルートCAは特権的に自己署名が信頼される様子。WebサイトのSSL/TLS証明書と検証は、まさにこの証明書チェーンの検証を行なっている。 SAML連携において、SPとIdPはこうしたチェーン検証を行わず「指定された証明書を信頼する」という形で個別に相手の証明書を事前登録するだけなので、CA発行か自己署名かどうかは実質的にセキュリティ上の差を生まない。例えば次のようなシナリオでCA発行証明書が意味を持つ。 既に証明書の発行・更新・監査の仕組みが既にある場合、Snowflake固有の自己署名証明書を個別管理する対象として放置せず、他のシステムと同じ棚卸し、アラート、更新プロセスに組み込める。「Snowflakeだけの特殊な例外」として管理外になるリスクを排除できる。 ISO27001/SOC2などの監査で「証明書は組織の承認されたCAから発行されたものであること」という統制が求められている場合、自己署名証明書だと指摘対象になり得る。CA発行にしておけば、発行記録・Subject検証済みという裏付けを監査人に示せる。 NameID規格化 SAML2_REQUESTED_NAMEID_FORMAT SPからIdP方向の AuthnRequest 送信時、「私(SP)はこの形式のNameIDが欲しい」というIdPへのリクエストを送る。このフォーマットを SAML2_REQUESTED_NAMEID_FORMAT で指定する。 一方、IdPからSP方向の SAML Response受信時、Snowflakeが受け取ったNameIDの値をLOGIN_NAMEと文字列比較する処理において、単に文字列として比較されるだけなので、SAML2_REQUESTED_NAMEID_FORMAT に一致していなくて良い。 ログイン・ログアウト時の挙動 SAML2_POST_LOGOUT_REDIRECT_URL Snowflake Webインターフェースの「Log Out」ボタンをクリックした後にSnowflakeがユーザーをリダイレクトするエンドポイント。未設定の場合はSnowflakeの標準ログイン画面に戻る。IdPのログアウトページや社内ポータルなど任意のURLへ誘導したい場合に使う。ただしこれはSAML Single Logout(SLO、IdP側セッションも含めた一括ログアウト)ではなく、あくまでSnowflake側のログアウト完了後にブラウザを1回リダイレクトさせるだけの機能の様子。 SAML2_FORCE_AUTHN ユーザーが最初の認証フロー中に、Snowflakeへアクセスするための再認証を強制されるかどうか。TRUEにすると、ユーザーがAWSアクセスポータルに既にログイン済み(SSOセッションが生きている)でも、Snowflakeへのログインのたびに毎回IdPでの認証をやり直させる。NameIDの形式に関係なく設定可能。特に理由がなければデフォルトのままでよい。 SAML2_SNOWFLAKE_ISSUER_URL SnowflakeサービスプロバイダーのEntityID/Issuerを明示的に指定するパラメータ。省略時はSnowflakeが https://.snowflakecomputing.com をデフォルトとして自動生成する。DESC SECURITY INTEGRATION では、明示指定の有無にかかわらず現在の値が確認できる。 SAML2_SNOWFLAKE_ACS_URL IdPがSAML認証応答を送り返すSnowflakeのAssertion Consumer Service URLを明示的に指定するパラメータ。省略時はSnowflakeが https://.snowflakecomputing.com/fed/login をデフォルトとして自動生成する。DESC SECURITY INTEGRATION では、明示指定の有無にかかわらず現在の値が確認できる。 おわりに Snowflake SAML2.0 Security Integrationを使用したSP/IdP Initiated SAML Federationと構成について ドキュメントをまとめてみた。諸々煩雑な構成が必要だが、Snowflakeが綺麗にラップしていて、パラメタを理解して 設定するだけで、うまくSAML2.0フェデレーションを構成できることがわかった。

Streamlit appをrestricted caller’s rightsで動作させる場合にコンテナインスタンスが必須となる背景を考えた話

行レベルセキュリティに頼ったマルチテナントは、実装ミスによる認可破綻が起こりやすいので怖い。 RBAC・インフラレベルで認可する仕組みができれば、appは認可コードを一切かかずに、 Snowflakeに認可の安全性を移譲できる。 しかし、これまでStreamlit in SnowflakeではOwner\'s rightsでしか動作せず、 実現ができなかった。6月1日に「Restricted caller\'s rights」が一般提供(GA)され、 caller\'s rightsで Streamlit を動作させられるようになった。 ただし、コンテナインスタンスが必須となる。 どういう仕組みで機能するのか気になったので調べてみた。 [arst_toc tag=\"h4\"] Restricted caller\'s rightsが一般提供(GA)された これまで、ストアドプロシージャ、SPCSサービス、Streamlit in SnowflakeアプリはOwner\'s role、 すなわち、リソースの所有者の権限でしか動作させることができなかった。 2026年6月1日に「Restricted caller’s rights」がGAされたことで、 これらのリソースをCaller\'s role、すなわち呼び出し元権限で動作させることが可能となった。 呼び出し元の権限次第で、Snowflake側のガバナンスが全く効かない、という世界線は存在せず、 「呼び出し元の権限」に対して「別のロールによる許可」で制限した権限でのみ、これらを実行できる。 混乱する代理人問題とCaller\'s right Caller\'s rightで動作する、ということは、所有者が書いたコードを、閲覧者の権限で実行するということ。 閲覧者の強い権限により「閲覧するだけのつもりだったがDROPできてしまった」みたいなことになる。 これを混乱する代理人問題と言い、権限を持つ閲覧者が所有者のコードに意図せず権限を貸している。 アプリが悪意を持っていたりアホだったりした場合に被害が拡大する要因となる。 GRANT CallerとCaller\'s rightの権限波及の仕組み 実行主体が所有者(Owner)から呼び出し元(Caller)に移るためSnowflake側の権限波及が大きく変わる。 管理者が別の管理者に MANAGE CALLER GRANTS することで、別の管理者は CALLER GRANTS できる。 別の管理者は CALLER USAGE, CALLER SELECT 等で「このアプリ(Owner\'s role)が閲覧者の代理として 使って良い権限」をホワイトリスト形式で指定する。 結果、閲覧者がたまたま ACCOUNTADMIN でも 管理者はアプリが許可した対象へのSELECTしか 行使できない。つまり、閲覧者がDROPやALTERを持っていたとしてもアプリ経由では実行できない。 閲覧者の身元・ポリシー(行/列制御)は活かしつつアプリが行使できる権限の上限は管理者が固定する、 が実現できるようになる。 なぜウェアハウスランタイムではCaller\'s right動作ができないのか Restricted caller’s rights and Streamlit in Snowflake 公式によると、ウェアハウスランタイムではCaller\'s right動作できない。 By default, all Streamlit in Snowflake apps run with the privileges of the owner, not the privileges of the caller. The Streamlit app developer can define whether a container-runtime app runs with owner’s rights or restricted caller’s rights. Restricted caller’s rights aren’t supported in warehouse runtimes. Restricted caller’s rights requires Streamlit version 1.53.1 or later. 以下が核心。Streamlit in SnowflakeアプリはSQLをOwner権限のストアドプロシージャ相当のサンドボックスで動く。だから、DESCRIBE/SHOW/LIST 制限などの「ストアドプロシージャ相当の制約」が付く。 Runtime environments for Streamlit apps Warehouse runtimes provide an on-demand, personal instance of the Streamlit app for each viewer. When a viewer opens the app, a new instance of the app is created for that viewer. Each viewer has their own isolated environment, which increases user load times. While both runtimes execute SQL queries using the owner’s privileges, apps using warehouse runtimes are subject to similar restrictions as owner’s rights stored procedures. For more information, see Owner’s rights stored procedures. ストアドプロシージャ準拠の権限モデルにはCaller\'s rightが想定されていなかったため、 同じ理由でStreamlit in SnowflakeでもCaller\'s rightができなかった。 Caller\'s right の本質は、所有者とは別の閲覧者としてのセッションを張る、ということ。 これは、アプリプロセスが複数の認証コンテキスト(OwnerとCaller)を同時に保持・維持する必要がある。 ステートレスなアプリのホスティングでセッション変数が揮発する問題の解決に、 ウェアハウスにステートを持たせるという解決策をとらないのは、それは当たり前だし、 そういう用途なら、常時稼働でステートを持てるコンテナウェアハウスが適当だよ、ということかな。 ストアドプロシージャにおけるOwner\'s right、Caller\'s rightについては以下があたる。 [clink implicit=\"false\" url=\"https://docs.snowflake.com/en/developer-guide/stored-procedure/stored-procedures-rights\" imgurl=\"https://www.snowflake.com/content/snowflake-site/global/ja/lp/snowflake-ai-data-predictions/_jcr_content/root/responsivegrid/container_882109576/container/flexible_column_cont/flexible_column_content_container_1/container/image.coreimg.svg/1768274028304/snowflake-logo-blue.svg\" title=\"Understanding caller’s rights and owner’s rights stored procedures\" excerpt=\"A stored procedure runs with either the caller’s rights or the owner’s rights. It cannot run with both at the same time. This topic describes the differences between a caller’s rights stored procedure and an owner’s rights stored procedure.\"] まとめ Streamlit in Snowflakeのクエリ実行は背後でストアドプロシージャが動作していることを理解した。 ストアドプロシージャはセキュリティのため、制限されたサンドボックスで動作している。 混乱する代理人問題の対応として、SnowflakeではCaller\'s rightsによる動作を許可していない。 このモデルでは、Caller\'s rights の入り込む余地がないか、既存の破壊的変更が必要となる。 そこで、常時稼働、ステートフルを前提としたコンテナインスタンスに Caller\'s rights を ホワイトリスト式で制限する「Restricted caller\'s rights」が導入された。

やりなおし統計

Fellegi-Sunterモデルに基づく確率的名寄せパッケージ Splinkを試してみる

[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}

標本調査に必要なサンプル数の下限を与える2次関数

[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] 各地点において関数の値を最大にするベクトル((frac{partial f}{partial x_0},frac{partial f}{partial x_1}))を全地点に対して計算したものを勾配とかいう。 ある地点において、このベクトルの方向に向かうことにより最も関数の値を大きくする。 で、今後のために正負を反転して関数の値を最小にするベクトルを考えることにした。 関数の値を小さくする操作を繰り返していけば、いずれ\"最小値\"が見つかるはず。 というモチベを続けるのが勾配降下法。学習率(eta)を使って以下みたいに書ける。。 begin{eqnarray} x_0 = x_0 - eta frac{partial f}{partial x_0} \\ x_1 = x_1 - eta frac{partial f}{partial x_1} end{eqnarray} ということで(f(x_0,x_1)=x_0^2+x_1^2)の最小値を初期値((3.0,4.0))、 学習率(eta=0.1)に設定して計算してみる。 import numpy as np def numerical_gradient(f, x): h = 1e-4 grad = np.zeros_like(x) for idx in range(x.size): tmp_val = x[idx] x[idx] = tmp_val + h fxh1 = f(x) x[idx] = tmp_val - h fxh2 = f(x) grad[idx] = (fxh1 - fxh2) / (2*h) x[idx] = tmp_val return grad def gradient_descent(f, init_x, lr=0.01, step_num=100): x = init_x for i in range(step_num): grad = numerical_gradient(f,x) x -= lr * grad return x def function2(x): return x[0]**2 + x[1]**2 init_x = np.array([-3.0, 4.0]) v = gradient_descent(function2, init_x=init_x, lr=0.1, step_num=100) v # array([-6.11110793e-10, 8.14814391e-10]) ((0,0))に収束した。 ニューラルネットワークの勾配 損失関数を重みパラメータで微分する。以下みたいな感じ。 損失関数の大小を見るとして、例えば(w_{11})以外の重みを固定したとして(w_{11})をわずかに 増やしたときに損失関数の値がどれだけ大きくなるか。 損失関数の値はパラメータ(W)と入力(x)から決まるベクトルだけれども、それぞれ乱数と入力値が設定されている。 begin{eqnarray} W= begin{pmatrix} w_{11} & w_{12} & w_{13} \\ w_{21} & w_{22} & w_{23} end{pmatrix}, frac{partial L}{partial W}= begin{pmatrix} frac{partial L}{partial w_{11}} & frac{partial L}{partial w_{12}} & frac{partial L}{partial w_{13}} \\ frac{partial L}{partial w_{21}} & frac{partial L}{partial w_{22}} & frac{partial L}{partial w_{23}} end{pmatrix} end{eqnarray} 重み(W)が乱数で決まるネットワークがあるとする。このネットワークは入力と重みの積を出力 として返す。出力はSoftmaxを経由するとする。 ネットワークの出力と教師データのクロスエントロピー誤差を誤差として使う。 その前に、数値微分関数を多次元対応する。 普通、配列の次元が(n)個になると(n)重ループが必要になるけれども、 Numpy.nditer()を使うと(n)乗ループを1回のループにまとめることができる。 下のmulti_indexが((0,0),(0,1),(0,2),(1,0),(1,1),(1,2))みたいに イテレータが(n)次のタプルを返す。反復回数はタプルの要素数の直積。 Numpy配列にそのタプルでアクセスすることで晴れて全ての要素にアクセスできる。 def numerical_gradient_md(f, x): h = 1e-4 grad = np.zeros_like(x) it = np.nditer(x, flags=[\'multi_index\'], op_flags=[\'readwrite\']) while not it.finished: idx = it.multi_index tmp_val = x[idx] x[idx] = tmp_val + h fxh1 = f(x) # f(x+h) x[idx] = tmp_val - h fxh2 = f(x) # f(x-h) grad[idx] = (fxh1 - fxh2) / (2*h) x[idx] = tmp_val # 値を元に戻す it.iternext() return grad 初期値(x=(0.6,0.9))、教師データ(t=(0,0,1))をネットワークに入力する。 predict()は(1 times 3)を返す。 それをSoftmax()を通して、(t)とのクロスエントロピー誤差を求めたものが以下。 import numpy as np def cross_entropy_error(y, t): if y.ndim == 1: t = t.reshape(1, t.size) y = y.reshape(1,y.size) batch_size = y.shape[0] delta = 1e-7 return -np.sum( t * np.log( y + delta)) / batch_size def softmax(x): c = np.max(x) return np.exp(x-c) / np.sum(np.exp(x-c)) import sys, os sys.path.append(os.pardir) import numpy as np class simpleNet: def __init__(self): self.W = np.random.randn(2,3) def predict(self, x): return np.dot(x, self.W) def loss(self, x, t): z = self.predict(x) y = softmax(z) loss = cross_entropy_error(y, t) return loss net = simpleNet() x = np.array([0.6, 0.9]) p = net.predict(x) t = np.array([0, 0, 1]) net.loss(x, t) # 0.9463818740797788 このlossを(W)で微分したのが以下。 あえてパラメータ(W)を引数にとり損失関数の値を計算する(f(W))を定義することで、 数値微分が何と何の演算なのかをわかりやすくしている。 実際は(f(W))は(W)とは関係なく(x)と(t)だけから結果を返すけれども、 損失関数(f(W))を(W)で微分するという操作が自明になるようにコードを合わせている。 def f(W): return net.loss(x, t) dW = numerical_gradient_md(f, net.W) dW # array([[ 0.07627371, 0.49923236, -0.57550607], # [ 0.11441057, 0.74884853, -0.8632591 ]]) 結果の解釈 上記の(w),(W),(t)から(frac{partial L}{partial W})が求まった。 損失関数が何か複雑な形をしているという状況で、 (frac{partial L}{partial w_{11}})は(w_{11})がわずかに動いたときに損失関数の値が変化する量を表している。 それが(w_{11})から(w_{32})まで6個分存在する。 begin{eqnarray} frac{partial L}{partial W} = begin{pmatrix} frac{partial L}{partial w_{11}} & frac{partial L}{partial w_{21}} & frac{partial L}{partial w_{31}} \\ frac{partial L}{partial w_{12}} & frac{partial L}{partial w_{22}} & frac{partial L}{partial w_{32}} end{pmatrix} = begin{pmatrix} 0.07627371 & 0.49923236 & -0.57550607 \\ 0.11441057 & 0.74884853 & -0.8632591 end{pmatrix} end{eqnarray}

勾配の可視化

[mathjax] 2変数関数(f(x_0,x_1))を各変数で偏微分する。 地点((i,j))におけるベクトル((frac{partial f(x_0,j)}{partial x_0},frac{partial f(i,x_1)}{partial x_1}))を全地点で記録していき、ベクトル場を得る。 このベクトル場が勾配(gradient)。 (f(x_0,x_1)=x_0^2+x_1^2)について、(-4.0 le x_0 le 4.0)、(-4.0 le x_1 le 4.0)の範囲で、 勾配を求めてみる。また、勾配を可視化してみる。 まず、2変数関数(f(x_0,x_1))の偏微分係数を求める関数の定義。 ((3.0,3.0))の偏微分係数は((6.00..,6.00..))。 def numerical_gradient(f, x): h = 10e-4 grad = np.zeros_like(x) for idx in range(x.size): tmp_val = x[idx] x[idx] = tmp_val + h fxh1 = f(x) x[idx] = tmp_val - h fxh2 = f(x) grad[idx] = (fxh1 - fxh2) / 2*h x[idx] = tmp_val return grad def function2(x): return x[0]**2 + x[1]**2 p = np.array([3.0,3.0]) v = numerical_gradient(function2, p) v # array([6.e-06, 6.e-06]) (-4.0 le x_0 le 4.0)、(-4.0 le x_1 le 4.0)の範囲((0.5)刻み)で偏微分係数を求めて、 ベクトル場っぽく表示してみる。matplotlibのquiver()は便利。 各地点において関数の値を最も増やす方向が表示されている。 w_range = 4 dw = 0.5 w0 = np.arange(-w_range, w_range, dw) w1 = np.arange(-w_range, w_range, dw) wn = w0.shape[0] diff_w0 = np.zeros((len(w0), len(w1))) diff_w1 = np.zeros((len(w0), len(w1))) for i0 in range(wn): for i1 in range(wn): d = numerical_gradient(function2, np.array([ w0[i0], w1[i1] ])) diff_w0[i1, i0], diff_w1[i1, i0] = (d[0], d[1]) plt.xlabel(\'$x_0$\',fontsize=14) #x軸のラベル plt.ylabel(\'$x_1$\',fontsize=14) #y軸のラベル plt.xticks(range(-w_range,w_range+1,1)) #x軸に表示する値 plt.yticks(range(-w_range,w_range+1,1)) #y軸に表示する値 plt.quiver(w0, w1, diff_w0, diff_w1) plt.show() 値が大きい方向に矢印が向いている。例えば((-3.0,3.0))における偏微分係数は((-6.0,6.0))。 左上方向へのベクトル。 参考にしている本にはことわりが書いてあり、勾配にマイナスをつけたものを図にしている。 その場合、関数の値を最も減らす方向が表示されることになる。 各地点において、この勾配を参照することで、どちらに移動すれば関数の値を最も小さくできるかがわかる。

おっさんが数値微分を復習する

引き続き、ゼロDの写経を続ける。今回は、学生の頃が懐かしい懐ワード、数値微分。 全然Deepに入れないけれどおっさんの復習。解析的な微分と数値微分が一致するところを確認してみる。 昔と違うのは、PythonとJupyterNotebookで超絶簡単に実験できるし、 こうやってWordPressでLaTeXで記事を書いたりできる点。 [mathjax] まず、微分の基本的な考え方は以下の通り。高校数学の数3の範囲。 begin{eqnarray} frac{df(x)}{fx} = lim_{hrightarrow infty} frac{f(x+h)-f(x)}{h} end{eqnarray} 情報系学科に入って最初の方でEuler法とRunge-Kutta法を教わってコードを 書いたりレポート書いたりする。懐すぎる..。 または、基本情報の試験かなんかで、小さい値と小さい値どうしの計算で発生する問題が現れる。 ゼロDにはこの定義を少し改良した方法が載っている。へぇ。 begin{eqnarray} frac{df(x)}{fx} = lim_{hrightarrow infty} frac{f(x+h)-f(x-h)}{2h} end{eqnarray} 写経なので、がんばって数値微分を書いて動かしてみる。 簡単な2次関数(f(x))。 begin{eqnarray} f(x) &=& x^2 - 5x +3 \\ f\'(x) &=& 2x - 5 end{eqnarray} def numerical_diff(f, x): h = 10e-4 return (f(x+h) - f(x-h)) / (2*h) (f(x))と、(x=2.5)のところの(f\'(x))をmatplotlibで書いてみる。懐い... import matplotlib.pyplot as plt import numpy as np def f(x): return x**2 - 5*x + 3 x = np.arange(-10, 10, 0.1) y = f(x) dy = numerical_diff(f,x) plt.plot(x, y) x1 = -2.5 dy1 = numerical_diff(f, x1) y1 = f(x1) # y-y1 = dy1(x-x1) -> y = dy1(x-x1) + y1 j = lambda x: dy1 * (x-x1) + y1 plt.plot(x,j(x)) plt.xlabel(\'x\') plt.ylabel(\'y\') plt.grid() plt.show() 偏微分 2変数以上の関数の数値微分は以下の通り。片方を止める。 数値微分の方法は上記と同じものを使った。 begin{eqnarray} frac{partial f(x_0,x_1)}{partial x_0} &=& lim_{hrightarrow infty} frac{f(x_0 +h,x_1)-f(x_0-h,x_1)}{2h} \\ frac{partial f(x_0,x_1)}{partial x_1} &=& lim_{hrightarrow infty} frac{f(x_0,x_1+h)-f(x_0,x_1-h)}{2h} end{eqnarray} ((x_0,x_1)=(1,1))における(x_0)に対する偏微分(frac{partial f(x_0,x_1)}{x_0})、(x_1)に対する偏微分(frac{partial f(x_0,x_1)}{x_1})を求めてみる。 ちゃんと(frac{partial f(x_0,1.0)}{x_0}=2.00..)、(frac{partial f(1.0,x_1)}{x_1}=2.00..)になった。 import matplotlib.pyplot as plt import numpy as np from mpl_toolkits.mplot3d import Axes3D def f(x): return x[0]**2 + x[1]**2 X = np.meshgrid(np.arange(-5., 5., 0.2),np.arange(-5., 5., 0.2)) Z = f(X) fig = plt.figure(figsize=(6, 6)) axes = fig.add_subplot(111, projection=\'3d\') axes.plot_surface(X[0],X[1], Z) f0 = lambda x: x**2 + 1.0**2 f1 = lambda x: 1.0**2 + x**2 df0 = numerical_diff(f0, 1.0) df1 = numerical_diff(f1, 1.0) print(df0) # 2.0000000000000018 print(df1) # 2.0000000000000018 plt.show()