波士頓房價數據分析 (1)

數據描述

作為一個剛接觸數據科學的化工學生
本篇算是一個學習過程的小練習
我會提到如何使用MATLAB作簡單的線性迴歸分析
以及如何將高維度的數據“降維”以直觀解釋數據的內在結構

我使用現在機器學習領域很熱門的經典波士頓房價數據(Boston housing dataset)集做分析
這個數據集的用意是統計各種環境以及人為因素對房價的影響
總共有例如房間數，房屋稅率等總共13項影響房價中位數(MEDV)的因素

  1. CRIM:每人平均犯罪率
  2. ZN:居住區中高於25000平方英尺的停車場所佔土地比例
  3. INDUS:非零售企業佔地英畝數
  4. CHAS:查爾斯河亞變數(1代表包圍河，0代表沒有）
  5. NOX:氧化氮濃度(單位為千萬分之一）
  6. RM:平均房間數
  7. AGE: 1940年前建造之房東自用比例(屋齡）
  8. DIS:到五個波士頓工作中心的加權平均距離(上班通勤難易度）
  9. RAD:到達放射狀高速公路的難易程度(城內移動的難易度）
  10. TAX:房屋稅
  11. PTRATIO:學生與教師比例(學區好壞）
  12. B:與黑人占比呈正相關的數字(種族因素）
  13. LSTAT:低端人口比例
  14. MEDV:房價中位數

線性迴歸分析

在數據分析之前往往都要考慮拿到的數據集有沒有缺項或混亂的狀況
但是這套數據很完善故可跳過整理的步驟
如果我們想建立一個可以預測房價的數學模型
最簡單的就是線性迴歸了
如果我們將所有因素影響力以矩陣A表示並將房價以向量b表示
我們只需要找到關係式Ax = b中的係數矩陣x
在MATLAB中x可以透過簡單的regress函數達成
然而為了避免關係式中出現“非零截距”
可以在A矩陣的最後一列(column)後面加上所有元素為1的列

%% Import the Boston housing dataset 
close all;clc;clear all;
Data = readtable('Boston.csv');
% The median value of owner-occupied homes
b = table2array(Data(:,14)); 
% Other variables
A = table2array(Data(:,1:13));
A = [A ones(size(A,1),1)];

一般來說機器學習(其實粗略來說就是某種迴歸模型）會需要將數據分成“訓練集”以及“測試集”
但是作為一個簡單的練習我就省略了這個步驟

去除高度相關的影響因子

一個好的線性迴歸模型有一個重要的步驟就是要確保所有影響因素之間是互相獨立的
我們可以計算數據的相關矩陣(correlation matrix)
當數字的絕對值越接近1就代表越相關
如果兩個影響因素相關性太高(一般可以取大於0.7)
我們就可以考慮把“多餘”的因素剔除

%% Correlation matrix
figure(1)
% Calculate the correlation coefficients
R = corrcoef(A(:,1:13));
k = Data.Properties.VariableNames(1,1:13);
x_label = string(k);
xvalues = x_label;
yvalues = x_label;
% Plot the correlation matrix using a heatmap
h = heatmap(xvalues,yvalues,R);
h.CellLabelFormat = '%.2f'
h.Title = "Correlation matrix";
Boston housing dataset各因素間的相關係數矩陣

藉由上圖我們可以發現tax和rad兩項因素的相關性高達0.91
同樣我們可以發現indus, nox, age, dis和前兩項因素相關性相當高
我最後在這六項高度相關的因素中只保留了tax

%% Remove redundant columns
A(:,[3 5 7 8 9]) = [];
% Perform a linear regression 
x = regress(b,A);
% x is calculated such that Ax = b
%% Plot the MHV against neighborhood

原始數據與模型預測比較

接著我們就可以利用建立好的模型比較模型預測與原數據測試準確度
我也利用房價做了簡單的排序
可以發現預測並非完美
但尚可接受

% Unsorted 
figure(2)

subplot(1,2,1)
plot(b,'k-o'); hold on
xlabel('Neighborhood');
ylabel('Median home value');
plot(A*x,'r-o');
legend('Raw','Regression')
% Sorted 
subplot(1,2,2)
[b,sortind] = sort(b);
plot(b,'k-o');hold on
plot(A(sortind,:)*x,'r-o')
xlabel('Neighborhood');
ylabel('Median home value');
legend('Raw','Regression')
hold off
Boston housing dataset總共包含506戶的資料，線性迴歸可以用來建立預測模型。黑色為原始資料分佈；紅色是線性迴歸模型擬合的結果。右圖僅是左圖依房價做排序

t-SNE:降維數據可視化

接著我想知道這組14維的數據呈現什麼樣的內在結構
利用降維 (Dimensionality reduction)的方法
我們可以把高維度的數據轉換為容易讓人類了解的二維或三維數據並畫成圖
一般基礎的統計學會利用奇異值分解(singular value decomposition)
擷取出數據間變異程度最高的組成作為新的維度，
即為主成份分析(principal component analysis)，或簡稱PCA
然而PCA常常會使數據在低維度的表示中“擠成一團”
十幾年前有一個稱為t-distributed stochastic neighbor embedding (t-SNE)的方法問世
簡單來說就是利用學生t分佈的尾端比一般的常態分佈還要平緩

$$q_{ij}=\frac{{(1+{\left \| \mathbf{y_{i}-y_{j}} \right \|}^{2})}^{-1}}{\sum _{k}\sum _{l\neq k}{(1+{\left \| \mathbf{y_{k}-y_{l}} \right \|}^{2})}^{-1}}$$

藉由保持高維度和低維度轉換前後數據間的“距離分佈”來得到不至於擠成一團的低維度表示
上式就是自由度為1的學生t分佈表達式
有趣的是我們可以控制數據點周遭有多少“鄰居”
藉由改變稱為"perplexity"的超參數
我們就可以在低維度選擇要呈現局部(local)點的關係或者廣域(global)的數據結構
太大或太小的perplexity都不好

由於各種因素單位不盡相同
為了避免降維過程演算法偏袒某個因素
我先將個別因素做標準化(就是平均值為零，標準差為1的變換）

%% t-SNE dimensionality reduction
figure(3)

% Rescale the attribute matrix A 
A_rescaled = normalize(A);
A_rescaled(:,end) = [];

% Specify the hyperparameter "perplexity" for t-SNE
perplexity = [10 50 200];
len_p = length(perplexity);

% Plot t-SNE projection in low dimensions
for i = 1:len_p
    
    subplot(1,len_p,i)
    [Y,loss] = tsne(A_rescaled,'Perplexity',perplexity(i),...
    'NumDimensions',2);
    scatter(Y(:,1),Y(:,2),[],b,'filled');hold on
    colorbar
    colormap default
    xlabel('tSNE1');
    ylabel('tSNE2');
    title("Perplexity = " + num2str(perplexity(i)))
    fprintf('KL divergence is: %.3f',loss)

end



Boston housing dataset的t-SNE低維投影(點擊圖片可放大）


把房價高低用顏色來標註區別
並把原本13維的數據在二維平面表示
可以發現在適中的perplexity(50)下
t-SNE圖中大略有兩個大集團
分別對應高房價與低房價
可見總共13項影響因素(去除多餘的剩8項)有區別高端低端的能力
注意在t-SNE的演算中我並沒有標記最後的房價
是之後為了視覺化才上色的
顯然t-SNE能夠在未監督(unsupervised)的情況下找尋數據內在結構
並有分類的潛力
雖然有一些高房價的數據在低端集團中分不開(我猜是搞不清楚行情的房東害的？）

大約在2018年有一個最新的降維方法稱作
Uniform manifold approximation and projection (UMAP)問世
主要是由於單細胞定序(single cell sequencing)技術所需要而發展起來的
單細胞定序可以讓我們知道不同基因的表現量
並用降維的方式判斷出是什麼樣的細胞
這在發展生物學的實驗中幾乎是黃金準則
然而單細胞定序原始數據是成千上萬維
一個合適的降維方法就非常重要
t-SNE有許多缺陷
像是運行速度非常緩慢，並且廣域結構難以在低維度被確保等等
本篇僅為練習性質故沒有採用效能更佳的UMAP
有興趣深入了解可以參照Google AI支援的解說


最近寫畢業論文發現data science這個近年來火紅的領域對一般科學領域研究生來說挺有用的
有些概念能夠幫助詮釋實驗數據
或者將冷冰的數據轉化成未來簡化實驗流程的武器
好處多多






Comments

Post a Comment

Popular posts from this blog

大聖塔芭芭拉(Santa Barbara)地區私房景點推薦

Image
在SB的三年多陸續造訪許多或許只有當地人會比較清楚的景點，不需要人擠人才是真享受！ 本文介紹個人足跡以及背後的一些故事。 More Mesa, Goleta 如何到達？ 自Hollister Avenue轉進Cottage Hospital旁叉路，續行穿越腳踏車道後左轉到底即可。 建議將車停在Shoreline Dr.和Orchid Dr.轉角處。 介紹 是一塊寬廣的台地，南方緊鄰太平洋，在峭壁下方有非常狹窄的沙灘/礫灘。是一個人煙稀少的地方，有時可以看到人們在這裡泡步運動，以及在峭壁旁邊玩飛行傘。可選擇台地西南側地勢下凹處沿著峭壁下至海灘（不易發現，需要沿著海岸小路仔細觀察），或使用東南側靠近Hope ranch（遠看是一座小丘）林間的簡易木製樓梯。下至海灘後可以發現高低落差非常明顯，由於一般人不太會下到沙灘，這裡也被歸類為天體海灘，也可見到其他附近海岸不太常見的螃蟹或水鳥。 通往沙灘的下凹處 Matilija Trail, Ojai 如何到達？ 距離SB車程約一小時的深山，從Hwy 101轉Hwy 33後過Ojaha會有指標沿著陡坡向上，續行到底即可。 介紹 算是一個人不會非常多的景點，步道本身貫穿私人果樹農場，入口處可以看到農家養羊咩咩。要注意水量多的時期部分路段會有小溪橫切，如果怕鞋子會被弄濕請做好準備。山谷非常美，可以看見層理地形，置身果樹當中有種世外桃源之感。筆者本身沒有走很遠，如果有時間可以考慮在後方紮營，攀登更深山裡的瀑布。 Bill Wallace trail, El Capitan 如何到達？ 距離SB車程約20分鐘。沿著Hwy 101至117出口，延Calle Real走一段後到El Capitan Terrace Ln.入口處停車場停車 介紹 是個人最喜歡的輕鬆健行路線，雖然全線幾乎無遮蔽（請做好防曬），但是景緻極好可遙望101公路海岸線，春天三至四月會有小黃花可看。環狀路線下山會經過一處露營區。 North campus natural reserve, Goleta 如何到達？ El Colegio Rd.西側West Campus停車，再步行抵達 介紹 我之前住家附近一處許多人運動的寬廣區域，冬季多雨期潟湖會注滿水相當美麗，可賞鴨，清晨及傍晚可以在小路上看見出來覓食的小兔子。注意此處是保育類水鳥Snowy plover的保護區，請勿跨
Read more

[Upstream bioprocess development] Some notes about recent upstream process characterization approach

Image
This article provides a very brief introduction of some new aspects of biopharmaceutical process development. Implementation of quality by design (QbD) strategy has long been regarded as the paradigm in the development of biopharmaceutical production processes. An essential part of that is to identify the critical quality attributes (CQAs) and their relations between the critical process parameters (CPPs). Intensified design of experiments (iDoE) is a fairly new strategy (proposed during 2017) to efficiently reduce the number experiments required to characterize a bioprocess by bridging different experimental conditions in a more homogeneous way, although cautions might be needed to consider potential cell memory effects and metabolic shifts in mammalian cells. Once the process dynamics is captured, advanced process control techniques such as model predictive control (MPC) can then be adopted for optimal trajectory planning and closed-loop noise rejection during the production phase.
Read more

程序設計被忽略的一角：最大選擇率

Image
  (作者就讀於加州大學聖塔芭芭拉分校化工系博士班) 系上的Mike Doherty教授透過系辦寄了一封信給全體研究生要大家去參加他給大四學生上的一堂課。對他有點研究的人就會知道他是台大化工系美籍教授Jeff博士班的指導教授（也就是我的師公了)。Mike在化學製程設計的領域有不少貢獻，連著名的“Perry's Chemical Engineering Handbook”中關於蒸餾(Distillation)的章節也是他主導撰寫的，由他開設的大四程序設計課自然勾起我的好奇心。 令我訝異的是教室裡面也坐了許多其他化工系的教授，顯然Doherty要講的內容是許多教師們也不清楚的東西了。由於是程序設計課，他開頭就引用了20世紀中被拋出的一個問題: “有沒有辦法預測要如何設計一個有化學反應器和分離裝置的程序，使得想要的產物的選擇率(Selectivity)達到最大？”   對反應工程有點基礎概念的人就會知道根據有反應動力學的數據可以知道要如何選擇反應器以增加選擇率（例如要選用柱塞流PFR或連續攪拌CSTR或兩者的串聯並聯)，但這仍然只是一個粗略的估計，人們要的是找到一個“最佳”的配置。許多人包括Mike在內在為了解決第一個問題運用了一些幾何方面的論述想要解決這個問題，然而到了80年代中期學界對這個問題逐漸沒那麼在乎了，因為真正重要的問題這時浮現了： 有沒有辦法好，還要更好？ 於是上面的問題有了一個變體： “可不可以在只知道會有哪些反應發生的情況下，預測想要產物的“最大”選擇率，而不需要知道反應器和分離裝置如何安排的細節？” 精明一點的人看到這可能就會聯想到這個概念跟熱力學的 卡諾循環效率 有異曲同工之妙， 你可以理論上知道所有熱機(Heat Engine)能達到的最高能量利用效率(卡諾效率)，但你不會知道要怎麼達到那樣的效率。 以選擇率的問題來說，如果能計算出它理論上的最大值，或許就可以提供一個程序整合的準則，讓設計程序的人可以知道自己的設計離“最好”還差多遠。俄亥俄州立大學大學教授 Martin Feinbrug 在2000年左右在Industrial&Engineering Chemistry Research上發表了一篇名為"General Kinetic Bounds on Productivity and Selectivity in Re
Read more