首頁 / 技術文章 / 使用 MATLAB® 工具箱

NAG MATLAB®工具箱應用範例

使用 NAG MATLAB® 工具箱 - Part 3

本系列的文章可參考 http://tw.nag-gc.com/IndustryArticles/usingtoolboxmatlab.asp 以及 http://tw.nag-gc.com/IndustryArticles/usingtoolboxmatlabpart2.asp；本文中，我們將持續探索我們的工具箱， 並展示如何在 MATLAB 中使用我們的 NAG 程式庫，並利用 MATLAB 的圖形功能顯示結果。

注意： 本文中的範例皆是由我們的展示程式中截取的，若您直接以文中的 MATLAB 程式執行，也許無法完整的獲得如文中的結果。文中範例的完整的程式，請您自 此處下載。

G13 函數 - 時間序列分析

所謂時間序列指的是在某段時間中，收集潛在時間相依流程中的觀測點集合。NAG 的 G13 函數系列包含了許多調查與建構時間序列的統計模型； 由這些函數建構出的模型可以讓我們更了解資料的特性，或者是從時間序列中進行預測 (也就是對未來的行為進行預測)。 例如：一般說的 自我迴歸整合移動平均 (autoregressive integrated moving average) (ARIMA) 模型能夠預測時間序列，請看以下例子。

了解時間序列特性的第一步，是去計算它的 自相關函數，它描述了在不同時間點所存在的相關性 (依賴程度)。 將兩個不同時間段進行分割稱之為 lag，而自相關函數通常是依不同的 lag 值，以 自相關係數 來表示結果。 g13ab 函數能夠依許多基本的統計量，如平均值與變異數來計算這個係數。以下是程式片段：

    
% 時間序列資料 (來自太陽黑子觀測值)
x = [5;  11;  16;  23;  36; 
    58;  29;  20;  10;   8;
     3;   0;   0;   2;  11;
    27;  47;  63;  60;  39;
    28;  26;  22;  11;  21;
    40;  78; 122; 103;  73;
    47;  35;  11;   5;  16;
    34;  70;  81; 111; 101;
    73;  40;  20;  16;   5;
    11;  22;  40;  60;   80.9];

% nk 表示 lags 值，自相關函數所需
nk = int32(40);

% 呼叫 NAG 函數
[xm, xv, coeff, stat, ifail] = g13ab(x, nk);

因為 lag 是一個不連續的變數，自相關比較適合用直方圖表示 (或者稱 autocorrelogram)，如下圖：


圖 1：計算時間序列的自相關函數

自相關函數包含了潛在時間相依流程的定量與定性訊息；本例中，當以 ARIMA 模型去擬合時間序列，可看出大約具有 11 周期的擺動。 autocorrelation 圖形的形狀可以顯現出模型參數。圖 1 中，可以了解到這個時間序列為 非週期性的，因此還需進一步探索。 如果相關係數在前面幾個 lags 中的值很高，然後快速下降，便是我們一般稱的 移動平均 (MA) 序列， 若結果以波形方式呈現，便為 自迴歸 (AR) 序列。大部分的例子中，通常需要使用 ARIMA 模型 (同時擁有 AR 與 MA) 去擬合序列的資料。

除了自相關函數，我們也可以自 偏自相關 (partial autocorrelation) 函數 的圖形中，獲得其他更多的訊息。 我們可以在以上的例子中以 g13ac 函數取代 g13ab 函數的呼叫。

本範例的 MATLAB 腳本位於 NAGToolboxDemos/Time_series_analysis/g13ab_demo.m，也可下載 壓縮檔。

D01 函數 - 積分

在數值分析中，定積分中的一維或多維度的求值是非常普遍的。D01 函數提供了非常多種的演算法求解； 個別函數的呼叫全然依積分的形式而定。如果函數的形式可以分析出來，如被積函數包含了奇異值 (singularities)，尤其是代數或對數型的樣式，那麼 d01aj 函數最為合適。 其他較為特殊的函數，例如：如果被積函數為震盪函數，可以採用 d01ak；已知點為不連續的被積函數，則使用 d01al 函數。

d01aj 函數採用的是 自適應 演算法 - 意思是它會將要積分的函數，切分成許多的子區間， 並反覆切割，直到滿足所設定的精確度條件為止 (以下程式中以 epsabs 與 epsrel 變數來指定)。

以下程式呼叫 NAG 函數求解積分

    
epsabs = 0.0;
epsrel = 0.000001;

% lw 與 liw 是輸出矩陣 w 與 iw 的矩陣大小
% lw 必須介於 800 與 2000 之間；liw 為 lw 的四分之一
lw = int32(1000);
liw = lw/4;
iw = zeros(liw, 'int32');
w = zeros(lw, 'double');

% 呼叫 NAG 函數並驗證誤差
[result, abserr, w, iw, ifail] = d01aj(func, a, b, epsabs, epsrel);
if ifail ~= 0
    disp([' Non-zero ifail after d01aj call ', num2str(ifail)]);
    return;
end;

此程式中，d01aj 函數中所呼叫的 func 變數，指定的是 MATLAB m 檔案的函數字串名； 這個檔案必須是一個函數，將會傳回給定值的積分值。 以下是一個例子：

    
function [result] = myfunc(x)
    result = x^2*abs(sin(5*x));

(這是下圖中所計算的被積函數)。除了近似的積分估計值外，d01aj 函數還會輸出最後子區間的誤差值。 圖 2 顯示出每一個積分子區間的結果以及其相對應的誤差值。 從這個積分中也可以看出，在函數快速變化的區域，此演算法會窄化子區間的分割；且與較大誤差有關的子區間，其區間寬度也會較小 (函數接近於 0 的地方)。


圖 2：積分求值估計

本範例的 MATLAB 腳本位於 NAGToolboxDemos/Quadrature/d01aj_demo.m， 也可下載 壓縮檔。

G03 函數 - 多變量方法

多變量資料集乃是包含許多不同 目標 所觀察到的 變數。在科學分析中有許多這樣的範例資料，NAG 程式庫中的 G03 函數可以利用這些資料集來進行分析。 舉例說明，在英國的環境科學家，想要了解在他們所觀測的，分佈在英國與歐洲大陸間的 14 個區域中，進行觀測所取樣 300 個水鼠 (Arvicola 屬) 頭顱表現的 13 個特徵中，分屬於多少種的水鼠種類。 。自歐洲取得的資料已經歸類為二大類 (terrestis 與 sapidus)；科學家的任務便要來決定英國所取得的資料是屬於那一個種類的。

資料處理是以 14 個區域中所觀測的 13 個變數的平均值計算開始。這可以視為 13 維度空間中的 14 個點， 一般我們會用 主成分分析 (principal components analysis) - g03aa 函數 來分析這樣的資料。主成分分析是用來將原始資料維度簡化的一個方法；只要這些數值可以足夠表示原來的資料，那麼這些數值就可以取代原來的數值。 然而對於這個田鼠的資料而言，考量到前三個主成分只能解釋原始資料 65% 的變異，顯然是不適合的。

另一個可以應用在這個問題上的技術就是 metric scaling。第一步，便是建立一個 14 乘 14 的 相異矩陣 (dissimilarity matrix) ， 其中的每個元素就是其與原來 13 維度空間的距離。g03ea 函數可用來計算此相異矩陣。以下是程式片段：

    
% 每一個資料區域的名稱
label = char( 'Surrey', 'Shropshire', 'Yorkshire', ...
'Perthshire', 'Aberdeen', 'Eilean Gamhna', 'Alps', ...
'Yugoslavia', 'Germany', 'Norway', 'Pyrenees I', ...
'Pyrenees II', 'North Spain', 'South Spain' );

% 田鼠的資料。14 列 (每一個區域)，13 行 (每一個變數)
data0 = [
   [48.5  89.2   7.9  42.1  92.1 100.0 100.0 ...
    35.3  11.4 100.0  71.9  31.6   2.8];
   [67.7  67.0   2.0  23.0  93.0 100.0  86.0 ...
    44.0  14.0  99.0  97.0  31.0  17.0];
   [51.7  84.8   0.0  31.3  88.2 100.0  94.1 ...
    18.8  25.0 100.0  83.3  33.3   5.9];
   [42.9  50.0   0.0  50.0  77.3 100.0  90.9 ...
    36.4  59.1 100.0 100.0  38.9   0.0];
   [18.1  79.6   4.1  44.9  79.6 100.0  77.6 ...
    16.7  37.1 100.0  90.4   9.8   0.0];
   [65.0  81.8   9.1  31.8  81.8 100.0  59.1 ...
    20.0  30.0 100.0 100.0   5.0   9.1];
   [57.1  76.2  21.4  38.1  66.7  97.6  14.3 ...
    23.5   9.5 100.0  91.4  11.8  17.5];
   [26.7  53.1  23.5  38.2  44.1  94.1  11.8 ...
    11.8  18.2 100.0  94.9  12.5   5.9];
   [38.5  67.9  17.9  21.4  82.1 100.0  60.7 ...
    35.7  24.0 100.0  91.7  37.5   0.0];
   [33.3  83.3  27.8  29.4  86.1 100.0  63.9 ...
    53.8  18.8 100.0  83.3   8.3  34.3];
   [47.6  92.9  26.7  10.0  36.7 100.0  50.0 ...
    14.3   7.4 100.0  86.4  90.9   3.3];
   [60.0  90.9  13.6  68.2  40.9 100.0  18.2 ...
    100.0  5.0  80.0  90.0  50.0   0.0];
   [53.8  88.1   7.1  33.3  88.1 100.0  19.0 ...
    85.7   9.8  73.8  72.2  73.7   2.4];
   [29.2  74.0  16.0  46.0  86.0 100.0  18.0 ...
    88.0  16.3  72.0  80.4  69.6   4.0];
];

% 使用標準轉換，消除區間中不一致的變量差異的影響
data1 = asin(1.0 - 0.02*data0);

update = 'I';

% 距離為歐幾里得平方
dist = 'S';

% 對變數不做尺度轉換
scal = 'U';

% 觀察值數量 (收集資料的區域)
nobs = int32(14);

% 變數的數量 (每一區域所量測的數值)
nvar = int32(13);

% 哪一個變數並不會包含在距離計算中
isx = ones(1, nvar, 'int32');

% 每一個變數的縮放調整參數 (未使用，但是必須提供)
s = ones(1, nvar, 'double');

% 輸入距離矩陣 (未使用，必須提供)
d = zeros(nobs*(nobs-1)/2, 1, 'double');

% 呼叫 NAG 函數計算相異矩陣
[sOut, dOut, ifail] = g03ea(update, dist, scal, nobs, data1, ...
    isx, s, d, 'm', nvar);

第二個步驟便是找出一個新的矩陣能夠表示較少維度 (本例中是 3) 的資料距離；我們使用 metric scaling 將舊的矩陣投射到新的矩陣中。 接下來的 g03fa 函數，便是用來處理 metric scaling：

    

% 依變數數量將矩陣做標準化處理
dOut = dOut ./ double(nvar);

% 每個區域中有不同觀察數量，修正偏差
ns = [19 50 17 11 49 11 21 17 14 18 16 11 21 25];

jend = nobs-1;
for j = 1:jend
    istart = j+1;
    for i = istart:nobs
        index = (i-1)*(i-2)/2 + j;
        dOut(index) = abs(dOut(index) - (1.0/ns(i) + 1.0/ns(j)));
    end
end

% 進行 nvar 到 ndim 的 metric scaling 計算
roots = 'L';
ndim = int32(3);
[x, eval, ifail] = g03fa(roots, nobs, dOut, ndim);

產生的三維資料可用散佈圖顯示：

    
% 畫出 3 維資料的散佈圖
scatter3(x(1:nobs), x(nobs+1:2*nobs), x(2*nobs+1:3*nobs), ...
    40.0, colours, 'filled');

% 加上文字、顏色等屬性
dims = size(label);
for i = 1 : nobs
    text(x(i), x(nobs+i), x(2*nobs+i), label(i,1:dims(2)), ...
        'Color', [colours(i), colours(nobs+i), colours(2*nobs+i)], ...
        'FontSize', 13);
end

結果顯示在圖 3，圖中可以觀察出英國地區的資料 (黑點部分) 很接近藍點 (terrestis 種)，再來才是紅點 (sapidus 種)， 此意味著他們是屬於同一類的。


圖 3：田鼠資料度量尺度的散佈圖

本範例的 MATLAB 腳本位於 NAGToolboxDemos/Multivariate_methods/g03fa_demo.m， 也可下載 壓縮檔。

G05 函數 - 亂數產生器

在許多的計算領域中，亂數序列的可靠性是一件很重要的任務，例如蒙地卡羅模擬。G05 函數包含了許多不同的演算法，我們這裡說明其中的兩種演算法。

在亂數產生的領域中有兩個很明顯的區別，一個為 偽隨機 (pseudo)，一個為 准隨機 (quasi)。 所謂的偽隨機為統計屬性上較趨近於 "真實" 的隨機數 - 也就是說，這些亂數是由原始的亂數處理所取得 (計數器中所經過的時間頻率等等)。 例如：在偽隨機序列中每一個連續數字間的相關性是微乎其微的。相反的，准隨機則沒有這樣的特性 - 他們會在空間中提供一個較為平均的分佈值，這樣的特性很適合應用於蒙地卡羅模擬中，我們給定序列的長度，會獲得比偽隨機更正確的估計值。

我們的例子可以更清楚的說明。以下的程式會使用 g05fa 函數產生兩個偽隨機序列：

    
% (x,y) 產生的數量
n = int32(10000);

% 指定產生偽隨機亂數
g05za('O');

% 在 [a, b] 區間中，產生長度為 n 的二個向量 
a = 0;
b = 1;
[x] = g05fa(a, b, n);
[y] = g05fa(a, b, n);

以下程式中採用 g05yd 函數的 Faure 演算法，產生准隨機的二維亂數序列：

    
% 指定亂數向量的維度
idim = int32(2);

% 初始 Faure 演算法
[iref, ifail] = g05yc(idim);

% 產生二維向量
[quasi, irefOut, ifail] = g05yd(n, iref);

如果每個序列可以被解釋為二維空間中的個別點，那麼我們可用散佈圖表示 (如圖 4)，其中可以很明顯的顯示出兩種不同類型的亂數產生器的差異處。 雖然，他們似乎都可以隨機的分佈在二維的空間中，但是准隨機序列以更一致的方式呈現 (技術上來說，這是因為產生序列中的每一點時，演算法是以 最大化避免 的方式避免於其它數值接近)。


圖 4：兩種不同的亂數類別

本範例的 MATLAB 腳本位於 NAGToolboxDemos/Random_number_generation/g05fa_demo.m， 也可下載 壓縮檔。