∑ ∑ ∑ ∑

●【課題 4】の解説
１．理
論（地域システム論 A 講義資料 No.3 より）
ハフモデルの式に多重積形の誤差項をくわえ、その対数をとって、線形回帰モデル
に適用する。
OBS
ハフの原モデルを取り上げよう。行先選択確率の観測値を Pij
PRD
る予測値を Pij
、ハフモデルによ
ε
、誤差項を e ij とおく。m 居住地、n 商業地の場合を取り上げよう。
Sj
PijOBS = PijPRD e
ε ij
=
Tijλ εij
e
n S
j
(1.1)
∑Tλ
j =1 ij
式(1.1)の両辺の対数をとると、
Sj
log PijOBS
Tijλ εij
e
= log n
Sj
∑Tλ
j =1 ij
log PijOBS
log
n
= − log(∑
Sj
λ
j =1 Tij
PijOBS
) + log S j − λ log Tij + ε ij
= −α i − λ log Tij + ε ij , i = 1,..., m, j = 1,..., n
Sj
n
ここで、 i のみに依存する定数 log(
(1.2)
(*)
S
∑ T λj ) を、 α i とおいている。
j =1 ij
式(1.2)の最後の式(*)は、以下のようにして(単)線形回帰モデルに帰着できる。
(a) 自由度の考慮
OBS
確率 Pij
は、 j に関する総和が 1 であるので、自由に動ける数(自由度)は n − 1
である。そこで、最後のデータを一つ取り除く。
(b) 居住地 i を固定し、居住地ごとの回帰を実行する方法
居住地を i = 1 に固定する。変数を次のように対応させる。
Y j = log
P1OBS
j
Sj
, X j = log T1 j
パラメータ a1 , b1 を次のように対応させる。
a1 = −α1 , b1 = −λ
以上のもとで、居住地 1 に対する式(1.2)は次の単回帰モデルに帰着する。
Y j = a1 + b1 X j + ε1 j , j = 1,..., n − 1
同様に、他の居住地 i に対しても、個別の単回帰モデルを以下のように適用で
きる。
Y j = ai + bi X j + ε ij , j = 1,..., n − 1, for i
(1.3)
居住地別の単回帰モデル適用の問題点
•
距離抵抗の係数 λ は、理論上はすべての居住地で同一であるのに、居住
地ごとに異なってしまう。
次のモデルはこの問題点を解決したものである。
(c) ハフモデルの推定 (居住地を一括した単回帰モデルの適用)
居住地別の単回帰モデルの式(1.3)において、m コの居住地のサンプルをプール
することを考える。これは次のように、サンプルの番号付けを変えればよい。
(i, j ) → k = (n − 1)(i − 1) + j , i = 1,..., m, j = 1,..., n − 1
(1.4)
すなわち、変数の対応を次のようにつける。
Y( n −1)(i −1)+ j = log
PijOBS
Sj
, X ( n −1)(i −1)+ j = log Tij , i = 1,..., m, j = 1,..., n − 1 (1.5)
また、パラメータ bi は、居住地によらず一定で、 bi = b = −λ とおく。パラメー
タ ai は、居住地ごとの定数項に対する係数であったから、1,0 の値をとるダミ
i
ー変数と呼ばれる、次の変数 E , i = 1,..., m を導入する。
1 if [(k − 1) (n − 1)] + 1 = i
Eki = 
(1.6)
otherwise
0
以上のもとで、式(1.2)の最後の式である次式(1.7)は
log
PijOBS
Sj
= −α i − λ log Tij + ε ij , i = 1,..., m, j = 1,..., n − 1
(1.7)
式(1.4)、(1.5)、(1.6)を用いて、次の線形重回帰モデル(1.8)に帰着する。
Yk = a1Ek1 + " + am Ekm + bX k + ε k , k = 1,..., m(n − 1) (1.8)
２．プログラム
課題 4 では、(1.8)式から最小 2 乗法でパラメータを推定するプログラムを示している。
PROC IML;
RESET NOPRINT;
RESET LOG;
/** INPUT DATA **/
S={5000 15000 20000};
T={10 20 15,5 5 5,30 10 15};
POBS={0.3 0.22 0.48,0.13 0.36 0.51,0.03 0.61 0.36};
*仮設例の数値を入力している。①
*地区 1 ダミー;
E1={1 1 1, 0 0 0, 0 0 0};
*地区 2 ダミー;
E2={0 0 0, 1 1 1, 0 0 0};
*地区 3 ダミー;
E3={0 0 0, 0 0 0, 1 1 1};
*E1～E3 までは地区ダミー変数を定義している。②
YY=POBS[,1:2];
SS=REPEAT(S,3,1);
SS=SS[,1:2];
TT=T;
TT=TT[,1:2];
E1=E1[,1:2];
E2=E2[,1:2];
E3=E3[,1:2];
*ここまでは、選択確率の和が 1 であることを考慮して、最後のデータを 1 つ取り除いてい
る。③
YY=SHAPE(YY,1);
YY=YY`;
SS=SHAPE(SS,1);
SS=SS`;
TT=SHAPE(TT,1);
TT=TT`;
*ここでは、データを縦に並べている。④
Y1=LOG(YY/SS);
*上では、
PijOBS
Sj
の log をとっている。⑤
X1=LOG(TT);
*上では、 T の log をとっている。⑥
E1=SHAPE(E1,1);
E1=E1`;
E2=SHAPE(E2,1);
E2=E2`;
E3=SHAPE(E3,1);
E3=E3`;
*ダミー変数を縦に並べている。⑦
XX=E1||E2||E3||X1;
 PijOBS
 Sj

*最小二乗法の公式に合うようにダミー変数と log 

 を結合して、X を作っている。⑧


Y=Y1;
X=XX;
BETA =INV(X`*X)*X`*Y;
−1
*最小 2 乗法の公式 βˆ = ( X ′X ) X ′Y より、βを推定する。⑨
PRINT BETA;
３．プログラムの流れ
プログラムの流れを行列を使ってみていこう。
下の番号は、上のプログラムの解説の部分にある番号と対応している。
① 仮設例の数値を入力している。
A
B
1  0.3 0.22
P OBS = 2 0.13 0.36
3 0.03 0.61
A
B
S = [5000 15000
A
B
C
0.48
0.51
0.36 
C
20000]
C
1 10 20 15
T = 2  5 5 5 
3 30 10 15
②
ダミー変数を定義している。
A B C
1 1 1 1 
E1 = 2 0 0 0 
3 0 0 0 
③
A B C
1 0 0 0 
E2 = 2 1 1 1 
3 0 0 0 
A B C
1 0 0 0
E3 = 2 0 0 0 
3 1 1 1 
選択確率の和が 1 であることを考慮して、最後のデータを 1 つ取り除く。
A B
A
B
A
B
P
OBS
1  0.3 0.22 
= 2 0.13 0.36 
3 0.03 0.61
5000 15000 
S = 5000 15000 
5000 15000 
1 10 20 
T = 2  5 5 
3 30 10 
A B
1 1 1 
E1 = 2 0 0 
3 0 0 
④
A B
1 0 0 
E2 = 2 1 1 
3 0 0 
データを縦に並べておく。
 0.3 
0.22 
A
B

1  0.3 0.22  


0.13
P OBS = 2 0.13 0.36  → 

0.36 

3 0.03 0.61
 0.03


 0.61
 5000 
15000 
A
B

5000 15000  


5000
S = 5000 15000  → 

15000 

5000 15000 
 5000 


15000 
10 
 20 
A B
1 10 20   
5
T = 2  5 5  →  
5
3 30 10   
30 
 
10 
⑤
PijOBS
Sj
の log をとって、Y と定義する。
Y = log
P OBS
Sj
A B
1 0 0
E3 = 2 0 0 
3 1 1 
 −9.7212 
 −11.1299 


 −10.5574 
=

 −10.6375 
 −12.0238 


 −10.1101 
⑥
T の log をとる
 2.3026 
 2.9957 


1.6094 
log T = 

1.6094 
 3.4012 


 2.3026 
⑦
ダミー変数を縦に並べている。
1 
1 
1 1 1   
0 
E1 = 2 0 0  →  
0
3 0 0   
0 
 
0 
A B
⑧
0 
0 
1 0 0  
1 
E2 = 2 1 1  →  
1
3 0 0   
0 
 
0 
A B
0 
0 
1 0 0  
0 
E3 = 2 0 0  →  
0
3 1 1   
1 
 
1 
A B
重回帰モデルの公式にあうように X と Y を作る。
1  0   0   2.3026 
 −9.7212 
1  0   0   2.9957 
 −11.1299 


   

OBS










10.5574
0
1
0
1.6094
−
P
Y = log
=
 X = E1 E2 E3 log T =       

Sj
 −10.6375 
0  1   0  1.6094 
 −12.0238 
0  0  1   3.4012 


   

 −10.1101 
0  0  1   2.3026 
−1
⑨ 最小 2 乗法の公式 βˆ = ( X ′X ) X ′Y より、βを推定する。出力される結果は、X の並べ
て順に出てくるので、βは最後の-1.82 がβの値である。
 −5.59 
 −7.66 

 = ( X `X ) −1 X `Y
 −5.86 


 −1.82 

Download Report