Osaka の住宅価格の推定

Copyright (c) 2016, The MathWorks, Inc.
このデモでは、大阪の住宅価格を推定する回帰モデルを作成します。

データの読み込み

大阪市オープンデータポータルサイトから、大阪市の地価データをダウンロード
位置情報、部屋のサイズ、最寄り駅からの距離などの情報が住宅価格と一緒に入った Excel ファイルを Table 型として読み込みます。なお、住宅価格は一番最後の列の MedianValue という値です。
Housing = importfile_Osaka('data\Housing_Osaka_preprocessed.csv');
文字列の値が入っている ward (区)、 neighborhood (地域)、building (建物の建築材料) をカテゴリカル型に変更します。
カテゴリカル型にすることで、メモリの節約をすることができ、より大規模データを扱うことができます。
Housing.ward = categorical(Housing.ward);
Housing.neighborhood = categorical(Housing.neighborhood);
Housing.building = categorical(Housing.building);
Housing.shape_flag = categorical(Housing.shape_flag);

1. 線形回帰

fitlm 関数には一つの引数しか入れていません。デフォルトでは、テーブル型の変数 Housing の最後の変数が自動的に目的変数と認識されます。
カテゴリカル型の変数 (今回の場合は ward, neighborhood など) が含まれる際にはダミー変数を作る必要がありますが、fitlm などの回帰の関数は、自動的にカテゴリカル型の変数を認識してダミー変数を作成します。ダミー変数は数値を量的変数として扱うのではなく、Yes/No 2つの値のうちのどちらかしか取らないような質的変数です。今回の場合は、例えば ward_* というダミー変数が作成されています。
myFit = fitlm(Housing)
myFit =
線形回帰モデル: value ~ 1 + longitude + latitude + ward + ChangeRate + acreage + shape + road + transportation + neighborhood + shape_flag + usage_floor + building + build_house + build_store + build_office + build_warehouse + build_factory 推定された係数: Estimate SE tStat pValue ___________ __________ __________ __________ (Intercept) 1.382e+09 6.329e+08 2.1836 0.029728 longitude -8.7035e+06 4.7412e+06 -1.8357 0.067343 latitude -5.8181e+06 5.2051e+06 -1.1178 0.26451 ward_住之江 -8.4864e+05 4.5374e+05 -1.8703 0.062363 ward_住吉 -6.1866e+05 4.377e+05 -1.4134 0.15851 ward_北 -2151.9 2.1677e+05 -0.0099268 0.99209 ward_城東 1.6448e+05 3.5088e+05 0.46876 0.63957 ward_大正 -2.805e+05 4.2171e+05 -0.66515 0.50644 ward_天王寺 -4.6703e+05 2.8662e+05 -1.6295 0.10421 ward_平野 1.0591e+05 5.0899e+05 0.20807 0.83531 ward_旭 4.0191e+05 4.112e+05 0.97743 0.32911 ward_東住吉 -1.5645e+05 4.3445e+05 -0.36012 0.719 ward_東成 1.9176e+05 3.6259e+05 0.52887 0.59727 ward_東淀川 3.8433e+05 4.3228e+05 0.88908 0.37464 ward_此花 -4.0666e+05 4.2166e+05 -0.96442 0.33557 ward_浪速 -6.0681e+05 2.7675e+05 -2.1926 0.029065 ward_淀川 -2.357e+05 3.8871e+05 -0.60637 0.5447 ward_港 -6.9573e+05 3.6666e+05 -1.8975 0.058675 ward_生野 95599 3.5737e+05 0.26751 0.78925 ward_福島 -6.2707e+05 3.0434e+05 -2.0604 0.040177 ward_西 -6.3419e+05 2.3002e+05 -2.7571 0.0061709 ward_西成 -5.0041e+05 3.7467e+05 -1.3356 0.18264 ward_西淀川 -3.5762e+05 4.4067e+05 -0.81152 0.41768 ward_都島 -55551 3.3176e+05 -0.16744 0.86713 ward_阿倍野 -5.9797e+05 3.5119e+05 -1.7027 0.089609 ward_鶴見 3.6607e+05 4.2145e+05 0.86858 0.38573 ChangeRate 1.4472e+05 33205 4.3585 1.7735e-05 acreage -6.26 14.017 -0.4466 0.65547 shape -60186 51865 -1.1604 0.24675 road 8420.8 4216.7 1.997 0.046679 transportation -203.1 122.03 -1.6643 0.09704 neighborhood_住居地域 -88088 1.4904e+05 -0.59104 0.55492 neighborhood_商業地域 -6.8185e+05 2.2904e+05 -2.9769 0.0031365 neighborhood_工業地域 -4.2808e+05 2.3722e+05 -1.8046 0.072092 neighborhood_近隣商業地域 -3.9766e+05 2.9626e+05 -1.3423 0.18047 shape_flag_1 -2.0585e+05 2.1258e+05 -0.96834 0.33362 shape_flag_2 3.8982e+05 1.7214e+05 2.2645 0.024221 usage_floor 1.7729e+05 18338 9.6674 1.5467e-19 building_鉄筋コンクリート造 -6.2122e+05 1.5246e+05 -4.0746 5.8347e-05 building_鉄骨造 -1.7505e+05 1.1419e+05 -1.533 0.12627 build_house -5.8775e+05 1.4651e+05 -4.0117 7.5309e-05 build_store 3.2299e+05 1.2937e+05 2.4967 0.013046 build_office -6.9204e+05 1.5104e+05 -4.5819 6.6421e-06 build_warehouse -1.5451e+05 2.3704e+05 -0.65185 0.51497 build_factory 5.3827e+05 3.4877e+05 1.5433 0.12375 観測数: 361、誤差の自由度: 316 二乗平均平方根誤差: 7.37e+05 決定係数: 0.576、自由度調整済み決定係数 0.517 F 統計量と一定のモデルの比較: 9.75、p 値は 1.69e-37 です
myFit.NumPredictors
ans = 17
予測に使用した変数の数は17で、自由度調整済み決定係数は 0.517と、あまり良くない結果になっています。

2. ステップワイズ回帰

ステップワイズ回帰はどの変数がモデルを改善するのに重要であるかを 見つける特徴選択の回帰方法です。
(独立)変数を追加したときに 精度が良くなれば、その変数を追加します。
また、変数を抜いたときに、 抜かない場合と結果がほぼ変わらないのであれば その変数を削除します。
変数が重要だとわかった場合、重要な変数同士の交差項も追加されます。例えば2項目の ChangeRate:road は2つの変数の交差項も見ています。
myFit2 = stepwiselm(Housing)
1。usage_floor, FStat = 193.1921, pValue = 1.240299e-35 を追加中 2。building, FStat = 4.7571, pValue = 3.0985e-07 を追加中 3。ChangeRate, FStat = 13.1902, pValue = 2.97773e-06 を追加中 4。ChangeRate:usage_floor, FStat = 45.2047, pValue = 7.10298e-11 を追加中 5。build_store, FStat = 9.6853, pValue = 0.0020085 を追加中 6。ChangeRate:build_store, FStat = 27.5056, pValue = 2.70544e-07 を追加中 7。road, FStat = 6.0615, pValue = 0.014295 を追加中 8。road:usage_floor, FStat = 41.2926, pValue = 4.27345e-10 を追加中 9。build_house, FStat = 7.747, pValue = 0.0056721 を追加中 10。build_office, FStat = 16.6113, pValue = 5.68442e-05 を追加中 11。ChangeRate:build_office, FStat = 19.0337, pValue = 1.69462e-05 を追加中 12。usage_floor:build_house, FStat = 31.9461, pValue = 3.31069e-08 を追加中 13。build_house:build_office, FStat = 29.1907, pValue = 1.22303e-07 を追加中 14。usage_floor:build_office, FStat = 28.1639, pValue = 1.99871e-07 を追加中 15。ChangeRate:build_house, FStat = 8.9836, pValue = 0.0029222 を追加中 16。ChangeRate:road, FStat = 5.4449, pValue = 0.020203 を追加中 17。road:building, FStat = 5.8067, pValue = 0.0033128 を追加中 18。usage_floor:building, FStat = 4.4693, pValue = 0.012139 を追加中 19。shape_flag, FStat = 3.2109, pValue = 0.041558 を追加中 20。shape_flag:build_store, FStat = 3.7356, pValue = 0.024858 を追加中 21. ChangeRate:usage_floor, FStat = 0.022844, pValue = 0.87995 を削除しています
myFit2 =
線形回帰モデル: value ~ 1 + ChangeRate*road + ChangeRate*build_house + ChangeRate*build_store + ChangeRate*build_office + road*usage_floor + road*building + shape_flag*build_store + usage_floor*building + usage_floor*build_house + usage_floor*build_office + build_house*build_office 推定された係数: Estimate SE tStat pValue ___________ __________ ________ __________ (Intercept) -2.0991e+05 3.7317e+05 -0.56251 0.57414 ChangeRate 1.4068e+05 36878 3.8147 0.00016223 road -27804 18791 -1.4796 0.13991 shape_flag_1 -2.3565e+05 2.0102e+05 -1.1722 0.24193 shape_flag_2 79657 1.6245e+05 0.49034 0.62422 usage_floor 2.7304e+05 1.3897e+05 1.9647 0.050269 building_鉄筋コンクリート造 3.8532e+05 3.8424e+05 1.0028 0.31667 building_鉄骨造 2.8598e+05 3.409e+05 0.83889 0.40213 build_house 5.4124e+05 1.8135e+05 2.9846 0.0030479 build_store 41465 1.0135e+05 0.40912 0.68271 build_office 1.2121e+05 2.1453e+05 0.56501 0.57244 ChangeRate:road 4872.7 1417.8 3.4367 0.00066273 ChangeRate:build_house -1.0576e+05 36348 -2.9098 0.0038579 ChangeRate:build_store 46533 27638 1.6837 0.093173 ChangeRate:build_office -1.7766e+05 38109 -4.662 4.5246e-06 road:usage_floor 6490.2 1021.9 6.351 6.9452e-10 road:building_鉄筋コンクリート造 -38831 19837 -1.9574 0.051123 road:building_鉄骨造 -8557.4 18907 -0.4526 0.65113 shape_flag_1:build_store 4.8171e+05 3.0802e+05 1.5639 0.11879 shape_flag_2:build_store 5.6452e+05 2.4531e+05 2.3012 0.021991 usage_floor:building_鉄筋コンクリート造 -17034 1.3703e+05 -0.12431 0.90115 usage_floor:building_鉄骨造 -88154 1.3783e+05 -0.63958 0.52288 usage_floor:build_house -2.937e+05 36731 -7.9959 2.0932e-14 usage_floor:build_office -2.0699e+05 39089 -5.2955 2.1494e-07 build_house:build_office 1.508e+06 2.1091e+05 7.1499 5.4547e-12 観測数: 361、誤差の自由度: 336 二乗平均平方根誤差: 5.3e+05 決定係数: 0.767、自由度調整済み決定係数 0.75 F 統計量と一定のモデルの比較: 46、p 値は 3.62e-91 です
myFit2.NumPredictors
ans = 8
関数 stepwiselm を実行した結果は、予測に使用した変数の数は8で、自由度調整済み決定係数は 0.75と先ほどより少ない変数の数で良い結果を得られています。

学習データに対する予測の描画

figure;
plot(Housing.value, 'r.');
hold on
[ypred_sw,yci_sw] = predict(myFit2, Housing);
plot(ypred_sw, 'LineWidth', 2);
plot(yci_sw(:,1), '-', 'color', [0.7 0.7 0.7])
plot(yci_sw(:,2), '-', 'color', [0.7 0.7 0.7])
xlabel('x');
ylabel('y');
legend({'data','predictions', 'confidence interval (95%)'},'Location','Best');
hold off;
title('学習データに対する予測 (ステップワイズ回帰)')
価格がとても高い場所については学習サンプル数も少ないので推定があまりうまくいっていないことが分かります。この推定を改善する方法はないでしょうか。

予測結果を地図上に描画

地図描画に使用するファイル READHGT のあるフォルダをパスに追加します。
addpath(fullfile(pwd, 'readhgt'))
このファイルは MATLAB Central File Exchange で提供されています。
取得したい画像の緯度経度の範囲を指定します。
latlim = [34.55 34.8];
lonlim = [135.4 135.6];
画像を取得して描画します。
readhgt(latlim(1):latlim(2), lonlim(1):lonlim(2))
xlim(lonlim)
ylim(latlim)
h = gca;
delete( h.Children(1:end-1) )
h.Visible = 'on';
hold on
idx = ~isnan(myFit2.predict);
predict_res = myFit2.predict;
s = scatter(Housing(idx,:).longitude, Housing(idx,:).latitude, 70, predict_res(idx), 'filled');
図の設定の変更
hold off
alpha(s, 0.7)
colorbar('eastoutside')
title('Osaka の住宅価格 (円/㎡)')
h.Title.FontName = 'Meiryo UI'; % フォントの設定
colormap('jet'); % カラーマップの設定

3. ガウス過程回帰による住宅価格予測

: それぞれの点に対するガウスカーネルの尖り具合。大きいと隣接する訓練点が互いに強く影響を及ぼし合うため、精度は下がるが汎化性能は上がる。小さい、特に訓練点の間隔より小さいと、学習データにのみフィットする不自然な結果になる。
: 誤差の許容度。大きいと学習データに対して従おうとするので過学習しやすくなる。小さいと、学習データに対する精度は低くなるが汎化性能は上がる。
gprMdl = fitrgp(Housing, 'value');
[ypred_gp,~,yci_gp] = predict(gprMdl, Housing);

学習データに対する予測の描画

figure;
plot(Housing.value, 'r.');
hold on
plot(ypred_gp, 'LineWidth', 2);
plot(yci_gp(:,1), '-', 'color', [0.7 0.7 0.7])
plot(yci_gp(:,2), '-', 'color', [0.7 0.7 0.7])
xlabel('x');
ylabel('y');
legend({'data','predictions', 'confidence interval (95%)'},'Location','Best');
hold off;
title('学習データに対する予測 (ガウス過程回帰)')

新規データに対する住宅価格予測

住宅情報の入ったテーブル型変数 Housing1を読み込みます。
load Housing1
1行目は大阪駅周辺、2行目は放出駅の緯度経度情報です。その他の項目は ward (区) 情報以外同じです。
市の中心部の駅(大阪駅)と、中心街から少し離れた駅(放出駅)で値段の違いがあるか見てみましょう。
OsakaStation = predict(gprMdl, Housing1(1,:));
round(OsakaStation)
ans = 738778
HanatenStation_close = predict(gprMdl, Housing1(2,:));
round(HanatenStation_close)
ans = 710173
1平方メータあたりの価格は、大阪駅周辺でおよそ73.9万円、放出駅周辺でおよそ71.0万円で、大阪駅周辺の方が高いことがわかります。
2行目と3行目は同じ放出駅ですが、駅からの距離が違います。駅からの距離が、2行目は50m、3行目は1000mになります。
HanatenStation_far = predict(gprMdl, Housing1(3,:));
round(HanatenStation_far)
ans = 664242
放出駅から50mの場合は71.0万円、1km の場合は 66.4万円となっており、駅からの距離が近いほうが価格が高くなっていることが分かります。