論理の流刑地

地獄の底を、爆笑しながら闊歩する

繰り返し比例調整法とその周辺について(その1)

Introduction

繰り返し比例調整法(iterative proprtional fitting, IPF)という方法がある。
古くは1940年の論文に遡る由緒正しき手法であるが、その汎用さゆえか割と近年の最新手法にも部分的に使われていたりしている。

IPFは二変数(カテゴリ変数)A, Bがありそのクロス表が調査などから実測値として得られているが、
A,Bの結合分布のありかたはそのままに、A,Bの各カテゴリの周辺度数自体は別の分布に合わせたい、というときに使う手法である。

仮に元のクロス表をS, Sの分布を維持して周辺度数を別の分布に合わせたテーブルをTとすると、任意のセルの選び方に対してオッズ比がSとTで保持される

自分も昔からこの方法は知っていたけど、
そもそものところとして「何で比例調整の反復が上のような目的を達成するのか」という仕組みがわからなかったので、少し調べたことをまとめとく*1

参考URL・書籍

  1. Deming , Stephan : On a Least Squares Adjustment of a Sampled Frequency Table When the Expected Marginal Totals are Known ※PDF注意, 1940年の論文
  2. Stephan : An Iterative Method of Adjusting Sample Frequency Tables When Expected Marginal Totals are Known※PDF注意, Stephan本人による1942年のレビュー&改良論文
  3. 読書日記: 読了:Little & Wu (1991) 標本から得たAxBクロス表を既知の周辺分布に合わせたい、標本にはバイアスがあることがわかっている、さあどうするか※"GOD"小野滋氏のブログより類似手法の比較手法のレビュー論文のメモ。
  4. 小野氏がレビューしてる論文(Little and Wu 1991) ※フリーではないので注意
  5. 山下(2015)『非線形計画法』朝倉書店

非線形計画法 (応用最適化シリーズ)

非線形計画法 (応用最適化シリーズ)

Deming & Stephan(1940)のアイディア

わかりやすいように表記は若干かえつつも、私の関心にとって重要なところだけメモする。

問題の定式化

二つのR×Sのクロス表P , \Piがある。
オッズ比はPの値を保持したままで、周辺度数を\Piに合わせるのが、行いたいことである。

各クロス表の各セルの(相対)度数をp_{ij} , \pi_{ij}とあらわすとすると、この問題は以下の重み付き二乗和を最小化するような\pi_{ij}を見つける問題となる。
 WSQ =\sum_{i=1}^R \sum_{j=1}^S \frac{(\pi_{ij} - p_{ij})^2 }{p_{ij} }
ただし、この最小化は以下の周辺度数に関する制約を有する。この制約は自由度だけ、すなわちR+S-1個存在する。
 \sum_{j=1}^S \pi_{ij} = \pi_{i*} \ , \sum_{i=1}^R \pi_{ij} = \pi_{*j}

これは線形制約をもつ非線形最適化問題であり、いわゆる制約付き凸二次計画問題(山本 2015: 91-102)としてとらえられる。

[凸二次計画問題の一般形]

目的:\frac{1}{2} \boldsymbol{x}^T Q\boldsymbol{x}+\boldsymbol{q} \boldsymbol{x}  の最小化
条件: (\boldsymbol{a}^i)^T \boldsymbol{x} - b_i = 0 , i \in E
    (\boldsymbol{a}^j)^T \boldsymbol{x} - b_j \leqq 0 , j \in G
※行列Qは正定値対照行列⇒fは凸関数

ラグランジュ乗数の導入とそこからの式変形

上の問題は、制約付き凸二次計画問題のなかでも不等式制約を持たないシンプルなものに属する。
等式制約を充足しつつ最小化の対象となる二乗和の勾配がゼロとなるようなx(ここで言えば\pi_{ij})を求める問題となるので、
これは(あの有名な)ラグランジュの未定乗数法で解けるはずである。

制約式に乗数\lambdaをかけてたした(ここではStephanとDemingの表記に合わせて符号は負にしてあるが)
f(\pi_{ij}, \lambda_{i}, \lambda_{j}) = WSQ - \lambda_{i} (\sum_{j=1}^S \pi_{ij} - \pi_{i*})- \lambda_{j} (\sum_{i=1}^R \pi_{ij} - \pi_{*j})
を、\pi_{ij}, \lambda_{i},  \lambda_{j}偏微分してゼロになるようなパラメータ推定を行えばよいのである。

すなわち、
\frac{\partial }{\partial \pi_{ij}} f(\pi_{ij}, \lambda_{i}, \lambda_{j}) = \frac{1}{p_{ij}}(\pi_{ij}- p_{ij}) -  \lambda_{i} -  \lambda_{j} = 0
\frac{\partial }{\partial \lambda_{i}} f(\pi_{ij}, \lambda_{i}, \lambda_{j})  = \sum_{j=1}^S \pi_{ij} - \pi_{i*} = 0
\frac{\partial }{\partial \lambda_{j}} f(\pi_{ij}, \lambda_{i}, \lambda_{j})  = \sum_{i=1}^R \pi_{ij} - \pi_{*j} = 0
三条件をみたすような推定値を求めることになる。

一つ目の条件を変形することで、\pi_{ij} = p_{ij}( 1 + \lambda_{i} + \lambda_{j} )...(A)が得られる。
またここで、\sum_{i=1}^R p_{ij} = p_{*j}, \ \sum_{j=1}^S p_{ij} =p_{i*} と表記することにすると、

 \pi_{i*} = \sum_{j=1}^S \pi_{ij} = \sum_{j=1}^S p_{ij}( 1 + \lambda_{i} + \lambda_{j} )
 \pi_{i*}= p_{i*} + \lambda_{i} p_{i*} + \sum_{j=1}^S p_{ij} \lambda_{j}
\therefore \ \pi_{i*} - p_{i*} = \lambda_{i} p_{i*} + \sum_{j=1}^S p_{ij} \lambda_{j}
ここから\lambda_{i}を求めると
\lambda_{i} = \frac{1}{p_{i*}}(\pi_{i*}  - \sum_{j=1}^S p_{ij} \lambda_{j}) -1

これを(A)に代入すると以下の式が得られる
\pi_{ij} = p_{ij}(1 + \lambda_{i} + \frac{1}{p_{i*}}(\pi_{j*} -\sum_{j=1}^S p_{ij} \lambda_{j}) -1)
\pi_{ij} = p_{ij}( \frac{\pi_{i*}}{p_{i*}} + \lambda_{i} - \frac{1}{p_{i*}} (\sum_{j=1}^S p_{ij} \lambda_{j}) )

ここで、\frac{1}{p_{i*}} (\sum_{j=1}^S p_{ij} \lambda_{j}) )というのはi行における\lambda_{j}の重み付き平均であるから、
それを\lambda_{j}^{i-AVG}と記すことにすると、
 \pi_{ij} =p_{ij}(\frac{\pi_{i*}}{p_{i*}} + \lambda_{i}+\lambda_{j}^{i-AVG})

上のようなi行に関する操作は、j列に関しても同様な操作が可能であるから、
 \pi_{ij} =p_{ij}(\frac{\pi_{*j}}{p_{*j}} + \lambda_{j}+\lambda_{i}^{j-AVG})
も成り立つ。

繰り返し比例調整法のlogic

 \pi_{ij} =p_{ij}(\frac{\pi_{i*}}{p_{i*}} + \lambda_{i}+\lambda_{j}^{i-AVG})
 \pi_{ij} =p_{ij}(\frac{\pi_{*j}}{p_{*j}} + \lambda_{j}+\lambda_{i}^{j-AVG})

が導出できたところで、Deming and Stephan(1940:439-441)によれば、
上2式における第2・3項を無視した
 \pi_{ij} =p_{ij}(\frac{\pi_{i*}}{p_{i*}})
 \pi_{ij} =p_{ij}(\frac{\pi_{*j}}{p_{*j}})
の手続きを繰り返すことで、誤差(第2・3項の占める割合)が徐々に減少していく、という論理である。
※元論文ではこれがなぜ収束するかまでは書いていない。。

とりあえず、反復比例調整法で実測値と期待値の周辺度数の比率をかけていくことが、
非線形計画法(というかラグランジュの未定乗数法)のロジックに基づくことがわかった。

別アプローチ:連立一次方程式として解く

この論文は主に調整法の提案論文として知られているので、上記の調整法がメインの論文だと思われがちであるが、
実際にD&S論文を読んでみると、(コンピュータ普及前の1940年当時としては)ラグランジュ未定乗数法から得られる後述の正規方程式をとくのが大変だということで
比較的単純かつ簡略な"simplified procedure"として比例調整を提案している(p.439

というわけでsimplifiedでない方法ももちろん提案している。
上の式変形の途中で
\lambda_{i} = \frac{1}{p_{i*}}(\pi_{i*}  - \sum_{j=1}^S p_{ij} \lambda_{j}) -1 .... (C)
\pi_{*j} - p_{*j} = \lambda_{j} p_{*j} + \sum_{i=1}^R p_{ij} \lambda_{i} ....(D)
という式が得られていた。

(C)を(D)に代入すると、
\sum_i p_{ij} ( \frac{1}{p_{i*}} (\pi_{i*} - \sum_{j} p_{ij} \lambda_j  )-1 ) + p_{*j} \lambda_j = \pi_{*j} - p_{*j}
\Leftrightarrow -p_{*j} + \sum_i \frac{p_{ij}}{p_{i*}} \pi_{i*} - \sum_i \frac{p_{ij}}{p_{i*}} \sum_j p_{ij} \lambda_j + p_{*j} \lambda_j =  \pi_{*j} - p_{*j}
\Leftrightarrow - \sum_i \frac{p_{ij}}{p_{i*}} \sum_j p_{ij} \lambda_j + p_{*j} \lambda_j =  \pi_{*j} - \sum_i \frac{p_{ij}}{p_{i*}} \pi_{i*}

たとえば第s列の場合(j=s)について考えると、
- \sum_i \frac{p_{is} p_{i1}}{p_{i*}} \lambda_1 +\cdots+ (p_{*s} - \sum_i \frac{p_{is} p_{is}}{p_{i*}} ) \lambda_s+ \cdots + - \sum_i \frac{p_{is} p_{i(S-1)}}{p_{i*}} \lambda_{S-1} =  \pi_{*s} - \sum_i \frac{p_{is}}{p_{i*}} \pi_{i*}
という式が得られるが、この式がS-1個だけ得られ、また変数\lambda_jもS-1個ある。

これは連立一次方程式として解けるので、\lambda_jの推定値が得られ、さらにそれを(C)に代入することで\lambda_iも推定できる。
それを(A)の\pi_{ij} = p_{ij}( 1 + \lambda_{i} + \lambda_{j} )に代入すると、\pi_{ij}も求まる。

後世に伝わったのは、上記の比例調整法だが、一応このような正規方程式を解くことでも推定できるとしている論文であったのである。

あくまでもマシンパワーが限られた時代に苦肉の策として生み出されたのがIPFであるので、現在は素直にこっちの方程式を解いた方がいいんじゃね?と思ったりするのだが、あまりそういうのは見たことがないのが不可思議である。

検算

一応ここまで読んでしまったので、ほんとに方程式を解いても求まるかRで書いて検算してみる。

以下が関数

solve_IPF <- function( tbl1 , tbl2){
  
  ##---型キャスト---##
  if("data.frame" %in% class( tbl1 ) ){
    tbl1 <- as.matrix(tbl1) #tbl1
  }
  if( "data.frame" %in% class( tbl2 )  ){
    tbl2 <- as.matrix(tbl2) #tbl2
  }
  
### ==== 前準備 ==== ####
#-- 相対頻度への変換  
from_n <- sum( tbl1)
to_n <- sum( tbl2)
p_tbl1 <- tbl1 / from_n
p_tbl2 <- tbl2 / to_n

#-- 期待度数など  
to_Mrgn_R <- rowSums(p_tbl2)
to_Mrgn_C <- colSums(p_tbl2)
from_Mrgn_R <- rowSums(p_tbl1)
from_Mrgn_C <- colSums(p_tbl1)

from_RateInRow <- diag(1/from_Mrgn_R) %*%  p_tbl1
  
### ==== 対角成分の取得 ==== ####
diag_vec <- calc_diagIPF( from_Mrgn_C,from_Mrgn_R, p_tbl1)
### ==== 非対角成分の取得 === ####
nondiag_mat <- calc_nondiagIPF( from_Mrgn_C,from_Mrgn_R, p_tbl1 )  

### ==== 従属ベクトルの取得 === ####
b_vec <- calc_vecb( to_Mrgn_C, to_Mrgn_R, from_RateInRow)

#連立方程式の形にして解く
mat_A <- nondiag_mat
diag(mat_A) <- diag_vec

r <- nrow(p_tbl1)
s <- ncol(p_tbl1)
matA_v2 <- mat_A[1:(s-1),1:(s-1)] 
b_vec_v2 <- b_vec[1:(s-1)]

#列に対応するラグランジュ乗数を求める
lambda_j <-  c(solve(matA_v2) %*% matrix(b_vec_v2, ncol=1),0)
#行に対応するラグランジュ乗数を求める
lambda_i <- calc_lambda_i( lambda_j , p_tbl1 , from_Mrgn_R, to_Mrgn_R)

#ラグランジュ乗数をつかって、各セルを推定する
multi_mat <- matrix(rep( 1, s*r ) , nrow=r ) + matrix(rep( lambda_i, s), nrow=r)+ matrix(rep(lambda_j, r) , byrow=T, nrow=r)

est_p_mat <- p_tbl1 *multi_mat
ret_tbl <- sum(tbl2 )* est_p_mat

return(list( p1 = p_tbl1 , p2 = p_tbl2 , p1_Mrgn_R = from_Mrgn_R , from_RateInRow = from_RateInRow, b = b_vec, diag =diag_vec ,nondiag = nondiag_mat,lambda_j = lambda_j,lambda_i = lambda_i, multi=multi_mat, IPF = ret_tbl))  

} #function


#--補助関数③:右辺のベクトルの計算--#
calc_vecb <- function(mgn_C, mgn_R ,rateRow){
ret_vec1 <- mgn_C #要素S個のベクトル
ret_vec2 <- as.vector(t(rateRow) %*% matrix(mgn_R, ncol=1)) #要素 S×R に R×1で S×1にしたのをベクトル化
return( ret_vec1 - ret_vec2)
} #func


#--補助関数④:対角成分の計算--#
calc_diagIPF <- function(from_mgn_C ,from_mgn_R, from_tbl){
  ret_vec1 <- from_mgn_C #第一項
  from_tbl_sq <- from_tbl **2 #R ×S
  denom <-   1/ from_mgn_R # R個
  ret_vec2 <- as.vector(t(from_tbl_sq) %*% matrix( denom , ncol=1))
  return(ret_vec1 - ret_vec2)
} #func

#--補助関数⑤:非対角成分の計算--#
calc_nondiagIPF <-  function(from_mgn_C ,from_mgn_R, from_tbl){
 denom <- 1/ from_mgn_R
 from_tbl_ij <- t(diag(denom) %*%  from_tbl) %*% from_tbl #要素(i,j)に実測値テーブルのi列とj列の掛け合わせが入っているもの(S×S)
 return(-1 * from_tbl_ij)
} #func

#--補助関数⑥:行に対応するラグランジュ乗数の推定--#
calc_lambda_i <- function(lambda_j, from_tbl ,from_mgn_R ,to_mgn_R){
  denom <- 1 /from_mgn_R #分母にあたるもの
  lambda_i <- denom *(to_mgn_R- as.vector(from_tbl %*% matrix(lambda_j, ncol=1))) - rep(1, length(to_mgn_R)) # R×S とS×1
  return(lambda_i)
} #func


当然このように計算したものは、比例調整法による推定値と一致.....しなかった
あれ?って思って小野氏のブログを読んで気づいたのだが、実はD&S(1940)は間違っているとStephanは後に訂正している。

it was stated that "the final results coincide with the least squares solution."
This statement was incorrect , for although the adjusted values satisfy the condition equations,
they do not satisfy the normal eqquations and hence provide only an approximation to the solution.
( Stephan 1942 :166)

簡略に訳すと、「D&S論文では間違ったことを言いました。比例調整法によって得られる推定値は、正規方程式の解とは一致しませんでした。近似でしかないです」とのこと。
いやいや.....*2

なぜ一致しないか?

Littele and Wu(1991)が後に指摘している通り、反復比例調整法で得られる推定値は、以下のモデルを想定している。
ln(\frac{\hat{\pi}_{ij}}{p_{ij}} )= \hat{\mu} + \hat{\alpha_i}+  \hat{\beta_j}.....(M1)
(ここで\hat{\alpha_i},\hat{\beta_j}は先ほどの説明におけるラグランジュ乗数の推定値\lambda_i, \lambda_jにあたる。\hat{\mu}は全体効果)

いっぽうDeming and Stephen(1940)の正規方程式の解は、先にみたように以下のモデルを想定している
\frac{\hat{\pi}_{ij}}{p_{ij}}= \hat{\mu} + \hat{\alpha_i}+  \hat{\beta_j}......(M2)

簡単な検算で確認できるが、対数をとったモデルM1は元の実測テーブルと推定値の間でオッズ比は保存されるが、
方程式から得られるM2では、元の実測テーブルと推定値の間でオッズ比が変動してしまう。

実際にM2に基づいてオッズ比の変動倍率を計算して元のオッズ比にかけると、正規方程式の解によって推定されたクロス表におけるオッズ比と一致した。

Little and Wu(1991:87)のレビューによれば、実は反復比例調整法(rakingともいうらしい)は、
上の重み付き二乗和ではなく、以下の値を最小化しているのだとIreland and Kullback(1968)で指摘されているらしい。

I(\boldsymbol{\pi} ; p) = \sum_i \sum_j \hat{\pi_{ij} } ln(\frac{\hat{\pi_{ij}}}{p_{ij}} )

なんか見覚えがある形だ。
実は(論文の第二著者の名前からも推測できるように)これはいわゆるカルバック・ライブラー距離である。
反復比例調整法はKL距離の最小化を行っているアルゴリズムであったのだ!

まぁ、半世紀以上前の相当古い論文を読むときはその内容の妥当性にも気をつけようねという話でした....

というわけでEnjoy!!

*1:その2があるかはわからない

*2:でも当時は検算にも今より手間かかったからしゃーないのかな