一般化パレート分布の平均と分散をレベルセットを使って計算する

数理統計

積分の計算をレベルセットの積分に変形する公式として, Layer Cake Representation が知られている.

Layer Cake Representation

定理 1. (Layer Cake Representation)

\nu を非負実数直線 [0,+\infty) 上のBorel 集合族の上の測度とし, 単調増大関数
\begin{equation}\label{eq_LCR_1}
\phi(\lambda) = \nu ([0,\lambda))
\end{equation}
が任意の \lambda>0 に対して有限であるとする. また, (\Omega,\Sigma, \mu) を \sigma -有限な測度空間とし, f を非負な可測関数とする.
このとき,
\begin{equation}
\int_\Omega \phi(f(x))\,d\mu(x)
= \int_0^\infty \mu\left(x\in \Omega; f(x)>\lambda\right)\,d\nu (\lambda)
\end{equation}
が成り立つ. 特に, a>0 に対して \phi(\lambda) = \lambda^a とすると,
\begin{equation}
\int_\Omega f(x)^a\,d\mu(x)
= a\int_0^\infty \lambda^{a-1}\mu\left(x\in \Omega; f(x)>\lambda\right)\,d\lambda
\end{equation}
が成り立つ.

注意 2.

  1. \mu は符号付き測度でもよい. このとき, \phi は有界変動な関数であり, 第2式は正負のいずれかが有限確定の場合に成り立つ.
  2. 第3式を確率論っぽく書くと, f の代わりに確率変数 X を用いて,

\begin{equation}
E[X^a]
= \int_\Omega X(\omega)^a\,dp(\omega)
= a\int_0^\infty \lambda^{a-1} P(X>\lambda)\,d\lambda
\end{equation}
と書ける. ただし P\left(\omega\in \Omega; X(\omega)>\lambda\right) を P(X>\lambda) と表している.

続きを読む

「割った余り」での割り算は、既約分数にしなくても答えは変わらない

アルゴリズム 整数

この記事では, P を素数とし, \textrm{mod } P の下での割り算について述べます.

整数 a,b に対する \textrm{mod } P の下での割り算 a\times b^{-1} では,
分数 a/b が既約分数になっていなくても計算結果は変わりません.

そのことを述べていきます.

復習

wakabame.hatenablog.com
でも触れましたが, まずは逆元とは何かを述べます.
P の倍数ではない a=0,1, \cdots に対して, a の法 P の下での逆元 a^{-1} を求めましょう.
まず,

x =0,1, \cdots , P-1 が a の法 P の下での逆元であるとは,
ある y_a \in \mathbb{N} が存在して,
\begin{equation}
ax+Py_a =1
\end{equation}
を満たすこと

と定式化できたのでした.
そのような x は a が P の倍数でない限り一意に存在し, 以後これを a^{-1} と表します.

このような a^{-1} がなぜとれるかについては,
拡張ユークリッドの互除法について述べられている @drken さんの記事
qiita.com
を参考にしてみてください.

示すべきこと

a\times b^{-1} \textrm{mod } P は既約でなくても計算結果が変わらないことを示すために,

P 未満の自然数 k に対して, ka\times (kb)^{-1} を計算して, その差を観察しましょう.

まずは, 拡張ユークリッドの互除法から b^{-1}, (kb)^{-1} はある整数 y_b, y_{kb} を用いて,
\begin{equation}
b^{-1} = \dfrac{1-Py_b}{b},\quad (kb)^{-1} = \dfrac{1-Py_{kb}}{kb}
\end{equation}
と表せたことに注意します.

これにより,
\begin{align*}
ab^{-1} - ka(kb)^{-1}
&= a\dfrac{1-Py_b}{b} -ka \dfrac{1-Py_{kb}}{kb}\\
&= \dfrac{a(y_{kb}-y_{b})}{b}P
\end{align*}

を得ます.

ここで, 左辺は整数だから右辺も整数であり, b は P の倍数ではないので
\dfrac{a(y_{kb}-y_{b})}{b}
は整数です.

したがって両辺は P の倍数であることから,
\begin{equation}
ab^{-1} \textrm{mod } P = ka(kb)^{-1} \textrm{mod } P
\end{equation}

が従います.

別の説明

積 kb についての逆元は,
\begin{equation}
(kb)^{-1} \equiv b^{-1}k^{-1} \pmod P
\end{equation}
を満たします.

これにより,
\begin{equation}
ka(kb)^{-1} \equiv kab^{-1}k^{-1} \equiv ab kk^{-1} \equiv ab\pmod P
\end{equation}
を示すこともできます（指摘されるまで気が付きませんでいた。。）.

プログラミングコンテストにおける C++ での標準入出力のまとめ

C++ 競技プログラミング

入力

1行に1要素の場合

#include <bits/stdc++.h>
using namespace std;

int main() {
  int a;
  
  cin >> a;
}

入力が小数、文字列の場合は, int の部分をそれぞれ double, string にすればよい

1行にスペース区切りで複数要素の場合

  • 2つの数字を別の変数 N, M に
#include <bits/stdc++.h>
using namespace std;

int main(){
  int N, M;
  
  cin >> N >> M ;
}

入力が小数、文字列の場合は, int の部分をそれぞれ double, string にすればよい

  • 複数の数字をベクトル vec に
#include <bits/stdc++.h>
using namespace std;
// int N が既に与えられているとする. 

int main(){
  vector<int> vec(N);
  for (int i = 0; i < N; i++) {
    cin >> vec.at(i);
  }
}

C++ の標準入力においては半角空白も改行も同列に区切り文字と見なされるため, 入力が

0 1 2

であっても, 

0
1
2

であっても動作に影響はない.

array を使う方法もあるが、AtCoderは初心者はvector を使うほうがいいと述べている.
https://atcoder.jp/contests/apg4b/tasks/APG4b_n

出力

1要素の出力

  • 文字列または数字 A を改行つきで表示
#include <bits/stdc++.h>
using namespace std;

int main(){
  
  cout << A << endl;
}

複数要素の出力

  • ベクトル vec を1要素ごとに改行して出力
#include <bits/stdc++.h>
using namespace std;

int main(){
  
  for(auto v: vec){
    cout << v << endl;
  }
}
  • リスト vec を空白区切りで出力
#include <bits/stdc++.h>
using namespace std;

int main(){
  
  for (int i = 0; i < vec.size()-1; i ++){
    cout << vec.at(i) << " ";
  }
  cout << vec.back() << endl;
}

みんなのプロコン2019 D-Ears をPython で解いてみる

AtCoder python アルゴリズム 動的計画法DP

問題文が迷わせに来ていますね.

全ての要素は非負ですので, りんごさんの要望を \left\{A_i\right\}_{i=1}^L, すぬけさんの石の置き方を \left\{B_i\right\}_{i=1}^L としたときのコスト
\begin{equation}
\sum_{i=1}^L |A_i-B_i|
\end{equation}
を最小化すればいいことがわかります.

すぬけさんの動き方はりんごさんが好きに決められるので,

\left\{A_i\right\}_{i=1}^L が与えられたとき, \left\{B_i\right\}_{i=1}^L を最適に選んだときの最小コスト

を求めればよいことが分かります.

つまり,

りんごさんの要望 \left\{A_i\right\}_{i=1}^L が先に与えられていて, すぬけさんの石の置き方を \left\{B_i\right\}_{i=1}^L はそれに応じてよいものを選べる.

という問題設定です.

ただし, \left\{B_i\right\}_{i=1}^L に対する制約は結構複雑です.
すぬけさんの移動に関しては, 順番を巻き戻しても石の配置は変わらないので出発点の座標 S, 終了点の座標 Tの間には
\begin{equation}
S\le T
\end{equation}
なる関係があるとして一般性を失いません.

そして, 出発点と終了点の間の座標に置かれている石はかならず奇数個になっています.
すぬけさんは同じ座標を行ったり来たりすることで, 任意の奇数を実現できることに注意しましょう.
これにより, A_i が奇数ならば 0 までコストを下げられ, 偶数ならば 1 までコストを下げられます.

また, すぬけさんは出発点よりも左に進むことができます.
すぬけさんが訪れる左端の座標を D とすると, D\le S であり, その区間に置かれている石はかならず 0 ではない偶数個になっています.
これにより, A_i が奇数ならば 0 までコストを下げられ, 0 でない偶数ならば 1 までコストを下げられます.
0 であれば, この区間の最小コストは 2 となります.

D >0 のときは, D よりも左にはすぬけさんは移動していないことになるため, その区間に置かれている石は 0 個です.
これにより, 区間 i におけるコストは A_i そのものになります.

今までの議論は終了点の右に対しても行えます.
すぬけさんが到達できる右端の座標を U とすると, これまでの議論は,

  • D \le S \le T \le U は自由に取ることができる.
  • i = 0,\cdots, D ならばその区間のコストは A_i.
  • i = D+1,\cdots , S ならばその区間のコストは A_i \mod 2.
  • i = S+1,\cdots , T ならばその区間のコストは (A_i+1) \mod 2.
  • i = T+1,\cdots , U ならばその区間のコストは A_i \mod 2.
  • i = U+1,\cdots , N ならばその区間のコストは A_i.

という風にまとめられます.
それぞれの区間の長さは 0 であってもよいものとします.

貪欲法が上手く働かない理由

例えば L = 5 として,
\begin{equation}
A = [1, 0,0,0,1]
\end{equation}
であれば右端は見捨てた B =[1,0,0,0,0] が最適（の１つ）であるのに対して,
\begin{equation}
A = [100,0,0,0,100]
\end{equation}
のような状況であれば B = [100,1,1,1,100 ] とするのが最適になり,
そういったものをどう処理するかを指定するのが難しくなります.

さて, 計算量を減らそう.

D \le S \le T \le U を決定させればそのときの最小コストが求められることが分かりました.
しかしながら, それらをループで処理すると計算量は O(L^4) となり計算が間に合いません.

ですので, 入力の長さ L に対する動的計画法（ＤＰ）を考えます.
すなわち, \ell = 0,1,\cdots, L を任意に取るごとに, \left\{A_i\right\}_{i=1}^\ell が入力であったときの最小コストを求めましょう.

しかしながら, 安直なＤＰは働きません.

というのも, \left\{A_i\right\}_{i=1}^\ell に対する状態の分割が \ell+1 では通用しないためです.
これを確かめるために, 先ほど挙げた A = [100,0,0,0,100 ] の場合を考えましょう.
今考案したＤＰにおいて, 入力が \ell =4 のときまでを考えると, D=S=T=0, U=1 が最適な配置であるのに対して,
\ell = 5 のときは, D=0,S=1, T=4, U=5 が最良となってしまいます.

というわけで, 入力の長さ \ell だけでなく, 現時点での右端の"戦略"を 5 通り保持する必要があります.
ただし, \ell をインクリメントする際には, 右端の"戦略" よりも新しいものだけが選べることに注意します.
具体的には,

  • コストが A_i となる i = 0,\cdots, D なる区間を状態 0.
  • コストが A_i \mod 2.となる i = D+1,\cdots, S なる区間を状態 1.
  • コストが (A_i+1) \mod 2となる i = S+1,\cdots, T なる区間を状態 2.
  • コストが A_i \mod 2.となる i = T+1,\cdots, U なる区間を状態 3.
  • コストが A_i となる i = U+1,\cdots, N なる区間を状態 4.

と呼ぶことにすると, 例えば状態 0 で考えた場合には楽で,

先頭 \ell +1までの文字列のなかで, 右端の状態が 0 のときの最小コスト
＝ 先頭 \ell までの文字列のなかで, 右端の状態が 0 のときの最小コスト
　　+ \ell+1 文字目を状態 0 で実行したときのコスト

となります. さらに状態 1 を考える場合には,

先頭 \ell +1までの文字列のなかで, 右端の状態が 1 のときの最小コスト
＝ 先頭 \ell までの文字列のなかで, 右端の状態が 0,1 のときの最小コスト
　 + \ell+1 文字目を状態 1 で実行したときのコスト

となり, \ell での状態を複数個使う必要があります.
以上のことをまとめると, 一般に, 次のような漸化式が成り立ちます.
\begin{align}
\textrm{DP}[\ell + 1][j]
&= \min_{k = 0,\cdots, j } \textrm{DP}[\ell + 1][k] + \textrm{cost}(j , A_\ell) ,\\
\textrm{cost}(j , A_\ell)
&= \begin{cases}
A_\ell &\quad \textrm{if } j =0,4,\\
A_\ell \mod 2 &\quad \textrm{if } j =1,3,\\
(A_\ell +1) \mod 2 &\quad \textrm{if } j =2
\end{cases}
\end{align}
あとはこれを実装します. 

L = int(input())
A = [int(input()) for i in range(L)]

DP = [[0 for i in range(5)] for j in range(L+1)]

def cost(s, a):
    if s == 0 or s == 4:
        return a
    elif s == 1 or s == 3:
        if a > 0:
            return a%2
        else:
            return 2
    else:
        return (a+1)%2

for i in range(L):
    DP[i+1][0] = min(DP[i][:1]) + cost(0, A[i])
    DP[i+1][1] = min(DP[i][:2]) + cost(1, A[i])
    DP[i+1][2] = min(DP[i][:3]) + cost(2, A[i])
    DP[i+1][3] = min(DP[i][:4]) + cost(3, A[i])
    DP[i+1][4] = min(DP[i][:5]) + cost(4, A[i])

print(min(DP[L]))

これでめでたくＡＣしますね
https://atcoder.jp/contests/yahoo-procon2019-qual/submissions/4228991

AtCoder の DPまとめコンテストを python で解いてみる(F問題まで)

この記事は書きかけです.

動的計画法(DP)について, AtCoder でコンテストが開かれました.

atcoder.jp
Python で解けたものから, 自分なりに解説をつけていきたいと思います.

A - Frog 1

atcoder.jp

考えること

足場 n, n =1,\cdots, N にたどり着くまでの最小コストを考えます.
これを a_n と表すことにすると,
\begin{equation}
\begin{cases}
a_1 &= 0,\\
a_2 &= |h_1-h_2|,\\
a_n &= \min\{a_{n-1} + |h_n-h_{n-1}|, a_{n-2}+ |h_n-h_{n-2}| \},&\quad \textrm{if } n=3,4,\cdots, N
\end{cases}
\end{equation}
という風に一般的に表されます.
つまり, a_n を n=3, 4,\cdots, N のときもとめるには以前の計算結果を使いまわすことで結果を得られます.
a_N を求めたいのに a_n, n=1,2,\cdots, N-1 をわざわざ求めるのは直観的には遠回りに感じますが,
a_{n-2}, a_{n-1} のそれぞれの値を保存できていれば a_{n} を瞬時に求めることが出来るので, うまくいきます. 

N = int(input())
h = list(map(int, input().split()))
 
DP = [0, abs(h[0]-h[1])]
for i in range(N-2):
    DP += [min(DP[-2] + abs(h[i]-h[i+2]), DP[-1] + abs(h[i+1]-h[i+2]))]
 
print(DP[-1])
続きを読む

【ベイズ推論】ディリクレ分布の規格化定数の導出について

カテゴリ分布のベイズ推論をしたい！

K\in \mathbb{N} を与えられた自然数とし, 母数 \pi = (\pi_1,\pi_2,\cdots, \pi_K)\in \mathbb{R}^K,
\sum_{k=1}^K\pi_k =1 , \pi_k \ge 0 ( k=1,2,\cdots, K)となるカテゴリ分布を考えます.

すなわちどれか 1 成分ののみが 1, それ以外が 0 となるような K 次元値の確率変数 S に対するカテゴリ分布
\begin{equation}
\mathrm{Cat}(S|\pi) = \prod_{k=1}^K \pi_k^{S_k}
\end{equation}
を考えましょう.

本記事では, これが指数型分布族に属することを手がかりに, 共役事前分布としてディリクレ分布が得られることを示し, ついでにディリクレ分布の規格化定数を求める計算を行います.

指数型分布族の定義と, 観測モデルから決定される共役事前分布については次の記事の記法を用います：
wakabame.hatenablog.com

ディリクレ分布については, 以下の記事群が参考になります：
ディリクレ分布 - Wikipedia
ディリクレ分布の意味と正規化，平均などの計算 | 高校数学の美しい物語

規格化定数の導出については正田 備也氏の
http://www.cis.nagasaki-u.ac.jp/~masada/DirDistNorm.pdf
を参考に行いました.

続きを読む