6. Những ước lượng tổng quát

Trở về Mục lục cuốn sách

Trong Mục 3, ta đã đề cập tới các phương pháp khác nhau để ước lượng tham số có thể nhận được ngay và đơn giản từ mẫu cho trước {x0, x1,  …, xn − 1}. Ở mục này, một phương pháp chung được xét đến; phương pháp này cho phép các ước lượng nhận được với các tham số tuỳ ý của phân bố xác suất. Phương pháp này được dựa theo nguyên tắc hợp lý cực đại, thể hiện trong Mục 6.1. Nguyên tắc này có thể được mở rộng cho việc mô hình hoá dữ liệu, trong đó thường một mẫu ba số {(x0, y0, σ0),  (x1, y1, σ1),   …, (xn − 1, yn − 1, σn − 1)} được cho trước. Nói chung xi là các điểm số liệu biểu diễn cho một tham số điều khiển nào đó mà ta có thể chọn trong mô phỏng, như là nhiệt độ của một chất khí. Ta coi như tất cả các giá trị xi đều khác nhau. Hệ quả là, mô phỏng được thực hiện ở n giá trị khác nhau của tham số điều khiển. Các điểm số liệu yi là trị trung bình của những lần đo (chẳng hạn là mật độ của khí) thu được từ những mô phỏng ứng với giá trị cố định xi của tham số điều khiển. Các giá trị σi là những thang sai số tương ứng. 1 Việc mô hình hoá dữ liệu nghĩa là ta muốn xác định được quan hệ y = y(x). Thông thường ta có sẵn kiến thức hoặc một giả sử nào đó về mối quan hệ này, tức là ta có một hàm kiểm định được gắn tham số, yθ⃗(x). Do vậy, tập hợp các tham số θ⃗ phải được điều chỉnh sao cho hàm yθ⃗(x) khớp “nhất” với mẫu. Việc làm này được gọi là khớp số liệu và sẽ được giải thích ở Mục 6.2. Phương pháp này cũng dùng được để so sánh nhiều hàm khớp khác nhau nhằm quyết định xem hàm nào đặc trưng cho mô hình thích hợp nhất.

6.1 Hợp lý cực đại

Ở đây, ta xét bài toán sau: Với một mẫu cho trước {x0, x1,  …, xn − 1} và một phân bố xác suất biểu diễn bởi một hàm khối xác suất pθ⃗(x) hoặc hàm mật độ xác suất fθ⃗(x), ta muốn xác định các tham số \vec{\theta}=(\theta_1,\ldots,\theta_{n_{\rm p}}) sao cho hàm khối hoặc hàm mật độ xác suất biểu diễn được số liệu “một cách tốt nhất”. Chữ này được đặt vào dấu nháy, bởi vì không có định nghĩa duy nhất nào khẳng định “tốt nhất” là gì, hoặc thậm chí một cách toán học để rút ra một chỉ tiêu phù hợp. Nếu coi như không hề có kiến thức gì trước về tham số, thì ta có thể dùng nguyên tắc sau:

Định nghĩa 26. Nguyên tắc hợp lý cực đại phát biểu rằng những tham số θ⃗ cần được chọn sao cho sự hợp lý của bộ số liệu, cho trước các tham số này, là lớn nhất.

Trong trường hợp một biến ngẫu nhiên rời rạc, nếu có thể coi rằng các điểm số liệu riêng là độc lập nhau, thì sự hợp lý của số liệu chỉ đơn giản là cho bởi tích những xác suất của từng điểm số liệu riêng lẻ. Từ đó có định nghĩa sau về hàm hợp lý

L(θ⃗) ≡ pθ⃗(x1)pθ⃗(x2)…pθ⃗(xn − 1) = ∏ i = 0n − 1pθ⃗(xi)

Với trường hợp liên tục, xác suất để trong một phép thử ngẫu nhiên nhận được đúng một mẫu nhất định sẽ bằng 0. Tuy vậy, với một tham số bất định ε nhỏ, thì xác suất để nhận được một giá trị trong khoảng [\tilde x-\epsilon,\tilde x+\epsilon]\text{Pr}(\tilde x-\epsilon \le X < \tilde x+\epsilon)= \int_{\tilde x-\epsilon}^{\tilde x+\epsilon} f_{\vec{\theta}}(x)\,dx \approx f_{\vec{\theta}}(\tilde x)2\epsilon. Vì 2ε tham gia chỉ như một hệ số, nên nó không liên quan đến việc xác định giá trị cực đại. Vì vậy, với trường hợp liên tục, ta xét hàm hợp lý sau đây

L(θ⃗) ≡ fθ⃗(x1)fθ⃗(x2)…fθ⃗(xn − 1) = ∏ i = 0n − 1fθ⃗(xi)

Để tìm ra cực đại của một hàm hợp lý L(θ⃗) theo cách giải tích, ta phải tính các đạo hàm bậc nhất lần lượt theo tất cả những tham số, rồi bắt chúng bằng 0. Vì việc tính đạo hàm của một tích cần phải áp dụng quy tắc tương ứng nên sẽ tiên hơn nếu ta xét hàm hợp lý loga

l(θ⃗) ≡ logL(θ⃗) .  (82)

Cách này đã chuyển tích của các hàm khối hoặc hàm mật độ xác suất của những điểm số liệu riêng lẻ thành một tổng, từ đó dễ tính được đạo hàm hơn. Ngoài ra, do loga là hàm đơn điệu, nên cực đại của hàm hợp lý cũng chính là cực đại của hàm hợp lý loga. Vì vậy, bộ tham số nào phù hợp “nhất” sẽ được xác định theo cách hợp lý cực đại từ hệ phương trình sau đây

\frac{\partial l(\vec{\theta})}{\partial \theta_k} \stackrel{!}{=} 0  (k=1,\ldots,n_{\rm p})  \label{eq:pd} (83)

Lưu ý rằng việc các đạo hàm bậc nhất bằng 0 chỉ đảm bảo rằng điểm cực trị đã đạt được. Hơn nữa, những phương trình này thường có nhiều nghiệm. Vì vậy, bạn phải kiểm tra hẳn ra xem những nghiệm nào là cực đại đại phương, và nghiệm nào là giá trị lớn nhất thực sự. Lưu ý rằng các ước lượng hợp lý cực đại, vì bản thân là hàm số của các mẫu, nên cũng là những biến ngẫu nhiên {\rm{}ML}_{\theta_k,n}(X_0,\ldots,X_{n-1}).

Một ví dụ để thử chơi là hãy xét phân bố luỹ thừa với hàm mật độ xác suất cho bởi PT (39). Nó có một tham số μ. Hàm hợp lý loga cho một mẫu {x0, x1,  …, xn − 1} trong trường hợp này là

\begin{aligned}  l(\mu) & = & \log \prod_{i=0}^{n-1}f_{\mu}(x_i) \\  & = & \sum_{i=0}^{n-1} \log  \left\{ \frac{1}{\mu} \exp\left( -\frac{x_i}{\mu} \right) \right\} \\  & = & \sum_{i=0}^{n-1}\left( \log\left\{\frac{1}{\mu}\right\}  -\frac{x_i}{\mu} \right)\\  & = & n \log\left\{\frac{1}{\mu}\right\} -\frac{n}{\mu} \overline{x}\end{aligned}

Lấy đạo hàm theo μ ta được:

\begin{aligned}  0 \stackrel{!}{=}\frac{\partial L(\vec{\theta})}{\partial \mu}  = n \frac{-1}{\mu^2} \mu - \frac{-n}{\mu^2} \overline{x}  = \frac{-n}{\mu^2} (\mu-\overline{x})\end{aligned}

Điều này ngụ ý \mu=\overline{x}. Rất dễ kiểm chứng rằng giá trị này tương ứng với một điểm cực đại. Vì giá trị kì vọng của phân bố lũy thừa đơn giản là \text{E}[X]=\mu, điều này phù hợp với kết quả ở Mục 3, ở đó đã cho thấy rằng trị trung bình mẫu là một ước lượng không chệch cho giá trị kì vọng.

Nếu ta áp dụng nguyên tắc hợp lý cực đại đối với phân bố chuẩn có các tham số μσ2, ta thu được (ở đây sẽ không trình bày, bạn có thể xem, chẳng hạn [Dekking (2005)]) những ước lượng hợp lý cực đại, gồm trị trung bình mẫu \overline{x} (ứng với μ) và phương sai mẫu s2 (ứng với σ2). Điều này nghĩa là (xem PT (54)) ước lượng hợp lý cực đại cho σ2 bị chệch. May mắn là ta biết rằng độ chệch này dần triệt tiêu với n → ∞. Thật sự, có thể cho thấy rằng, dưới những điều kiện khá bình ổn đối với phân bố nội tại, tất cả ước lượng hợp lý cực đại {\rm{}ML}_{\theta_k,n}(X_0,\ldots,X_{n-1}) cho tham số θk đều tiệm cận không chệch, nghĩa là

\lim_{n\to \infty} \text{E}[{\rm{}ML}_{\theta_k,n}] = \theta_k

Trái ngược với các trường hợp phân bố lũy thừa và chuẩn, có nhiều trường hơp mà tham số hợp lý cực đại không trực tiếp gắp với một ước lượng mẫu tiêu chuẩn. Hơn nữa, thậm chí {\rm{}ML}_{\theta_k,n} thường không thể xác định được theo cách giải tích. Ở trường hợp này, ta phải tối ưu hóa hàm hợp lý loga bằng cách số trị, chẳng hạn, bằng cách dùng các phương pháp tương ứng trong thư viện khoa học GNU (GSL).

Ta hãy xét phân bố Fisher-Tippett làm ví dụ, xem PT (43), nhưng dịch chuyển ngang để có cực đại ở x0 thay vì ở 0. Như vậy, ta có hai tham số cần chỉnh là λx0. Hàm cần được tối ưu (hàm mục tiêu), nghĩa là hàm hợp lý loga trong trường hợp này, phải có một dạng đặc biệt khi dùng những hàm cực tiểu trong GSL. Đối số thứ nhất của hàm mục tiêu có chứa những tham số cần chỉnh của hàm mật độ xác suất, nghĩa là véc-tơ đối số chính của hàm mục tiêu. Đối số này phải có kiểu gsl_vector, một kiểu riêng của GSL dành cho véc-tơ. Ta cần phải kèm theo <gsl/gsl_vector.h> trong mã lệnh, để dùng được kiểu dữ liệu này. Những véc-tơ như vậy được tạo lập bằng gsl_vector_alloc(), đặt giá trị cho những phần tử bằng gsl_vector_set(), truy cập giá trị phần tử bằng gsl_vector_get() và xóa véc-tơ bằng gsl_vector_free(). Cách dùng những hàm này sẽ tự diễn giải được qua những ví dụ sau đây, song bạn cũng có thể tham khảo tài liệu GSL [Galassi và nnk. (2006)].

Đối số thứ hai của hàm mục tiêu có chứa một con trỏ đến dữ liệu phụ thêm phục vụ tính toán hàm mục tiêu, nghĩa là mẫu trong trường hợp này. Như vậy, mẫu phải được lưu trong một vùng bộ nhớ. Để làm điều này, ta dùng kiểu dữ liệu sau:

typedef struct
{
  int       n;       /* number of sample points; */
  double   *x;                         /* sample */
}
sample_t;

Vì thực ra gói GSL lại chứa các hàm cực tiểu, trong khi ta cần tính cực đại, nên hàm hợp lý loga ở đây trả lại số đối của giá trị hợp lý loga. Hàm hợp lý loga sẽ như sau:

Trước tiên, ta chuyển đổi những con trỏ được truyền đến làm đối số, thành dạng dữ liệu phù hợp (các dòng 8–10). Tiếp theo, giá trị hợp lý loga thực tế
l(λ, xo) = nlogλ − λ∑ i = 0n − 1(xi − x0) − ∑ i = 0n − 1exp( − λ(xi − x0))
được tính ở các dòng 12–15 và cuối cùng trả lại sau khi đảo dấu (dòng 17).

GSL đã lập sẵn một số thuật toán cực tiểu. Chúng được xếp vào một trong hai nhóm. Nhóm thứ nhát là những thuật toán yêu cầu hàm mục tiêu và các đạo hàm thứ nhất của nó. Nhóm thứ hai gồm các thuật toán chỉ cần hàm mục tiêu. Ở đây ta dùng thuật toán simplex, thuộc nhóm thứ hai. Nó hoạt động bằng cách khai triển một simplex,2 lượng giá hàm mục tiêu tại các đỉnh của simplex này, rồi lặp lại việc điều chỉnh simplex đến khi nó rất nhỏ và chứa nghiệm số. Lưu ý rằng thuật toán chỉ có thể tìm được cực tiểu địa phương, và chỉ một cực tiểu như vậy. Nếu tồn tại nhiều cực tiểu thì việc chọn các tham số ban đầu sẽ ảnh hưởng mạnh đến kết quả; trường hợp đó thì có thể ta phải thử nhiều tham số. Để biết thêm chi tiết, hãy xem . Ở đây ta chỉ trình bày cách dùng công cụ tìm cực tiểu. Bản thân nó được lưu trong một cấu trúc dữ liệu đặc biệt có kiểu gsl_multimin_fminimizer. Hàm mục tiêu phải được đặt vào trong một biến “vây quanh” thuộc kiểu gsl_multimin_function. Hơn nữa, ta cần hai biến gsl_vector để lưu giá trị ước tính hiện thời của cực trị (chỉ vị trí của simplex) và để lưu kích thước của simplex. Đồng thời, par được dùng ở đây để chỉ số chiều của đối số hàm mục tiêu (2) và sample để lưu mẫu.

Các biến này được khai báo như sau:

  int num_par;                        /* number of parameters */
  sample_t sample;                                  /* sample */

  gsl_multimin_fminimizer *s;           /* the full mimimizer */
  gsl_vector *simplex_size;        /* (relative) simplex size */ 
  gsl_vector *par; /* params to be optimized = args of target */
  gsl_multimin_function f;  /* holds function to be optimized */

Việc cấp phát bộ nhớ và khởi tạo những biến này có thể như sau:

  sample.n = 10000;                          /* initilization */
  sample.x = (double *) malloc(sample.n*sizeof(double));
  num_par = 2;

  f.f = &ll_ft;                    /* initialize minimization */
  f.n = num_par;
  f.params = &sample;
  simplex_size = gsl_vector_alloc(num_par);  /* alloc simplex */
  gsl_vector_set_all(simplex_size, 1.0);      /* init simplex */
  par = gsl_vector_alloc(num_par);  /* alloc + init arguments */
  gsl_vector_set(par, 0, 1.0);
  gsl_vector_set(par, 1, 1.0);
  s =
    gsl_multimin_fminimizer_alloc(gsl_multimin_fminimizer_nmsimplex, 
                                  num_par);
  gsl_multimin_fminimizer_set(s, &f, par, simplex_size);

Cách thiết lập đối tượng cực tiểu hóa gồm hai bước, trước hết là cấp phát bộ nhớ bằng gsl_multimin_fminimizer_alloc(), sau đó là khởi tạo bằng gsl_multimin_fminimizer_set() cùng việc truyền vào hàm mục tiêu, điểm khởi đầu par và kích thước simplex (ban đầu).3 Mảng sample.x[] phải được điền kín bằng mẫu cụ thể (ở đây sẽ không chỉ ra).

Vòng lặp cực tiểu hóa có dạng như sau:

  do                                  /* perform minimization */
  {
    iter++;
    status = gsl_multimin_fminimizer_iterate(s);  /* one step */
    if(status)    /* error ? */ 
      break;
    size = gsl_multimin_fminimizer_size(s);    /* converged ? */
    status = gsl_multimin_test_size(size, 1e-4);
  }
  while( (status == GSL_CONTINUE) && (iter<100) );
        

Công việc chính được thực hiện trong gsl_multimin_fminimizer_iterate(). Sau đó, nó được kiểm tra xem liệu có lỗi xảy ra không. Tiếp theo, kích thước của simplex được tính ra và cuối cùng là kiểm tra xem kích thước này có giảm xuống nhỏ hơn một giới hạn nào đó, ở đây là 10−4.

Những giá trị ước tính của tham số có thể nhận được bằng gsl_vector_get(s->x, 0)gsl_vector_get(s->x, 1). Lưu ý rằng sau cùng tất cả bộ nhớ đã huy động phải được giải phóng:

Xét một ví dụ, trong đó n = 10000 điểm số liệu được phát sinh theo một phân bố Fisher-Tippett có các tham số λ = 3.0, x0 = 2.0. Với những tham số trên, phép cực tiểu hóa sẽ hội tụ về các giá trị \hat\lambda=2.995\hat x_0=2.003 sau 39 lần lặp.

6.2 Khớp số liệu

Ở mục trước, các tham số của một phân bố xác suất được chọn sao cho phân bố mô tả được đúng đắn nhất số liệu hiện có. Ở đây, ta xét một trường hợp tổng quát hơn, gọi là mô hình hóa số liệu. Như đã giải thích ở trên, ở đây ta có mẫu gồm ba số {(x0, y0, σ0),  (x1, y1, σ1),   …, (xn − 1, yn − 1, σn − 1)}. Thông thường yi là các giá trị đo được từ một mô phỏng với tham số kiểm soát nào đó (chẳng hạn nhiệt độ) cố định tại những giá trị xi riêng biệt; σi là thang sai số tương ứng với yi. Ở đây, ta muốn xác định các tham số \vec{\theta}=(\theta_1,\ldots,\theta_{n_{\rm p}}) sao cho hàm cho trước gắn tham số yθ⃗(x) phải khớp số liệu “nhất”, chẳng hạn ta muốn khớp hàm vào số liệu này. Tương tự như trường hợp khớp một hàm khối xác suất hoặc hàm mật độ xác suất, không có một quy tắc chung nào cho ý nghĩa của “tốt nhất”.

Ta hãy giả sử rằng yi là các biến ngẫu nhiên, nghĩa là so sánh giữa những lần mô phỏng khác nhau. Như vậy, các giá trị đo đạc sẽ phân tán xung quanh các giá trị “thực” yθ⃗(xi) tương ứng. Sự phân tán này có thể được mô tả gần đúng bằng một phân bố chuẩn với trị trung bình yθ⃗(xi) và phương sai σi2:

\begin{aligned}  \label{eq:distrMeasurements}  q_{\vec{\theta}}(y_i)\sim \exp\left(-\frac{  (y_i -y_{\vec{\theta}}(x_i))^2}{2\sigma_i^2} \right) \,.\end{aligned} (85)

Giả sử này thường đúng, chẳng hạn khi từng điểm mẫu yi bản thân nó là một trị trung bình mẫu thu được từ một mô phỏng thực hiện tại một giá trị tham số điều khiển xi, còn σi là thang sai số tương ứng. Hàm hợp lý loga của toàn bộ mẫu số liệu là

\begin{aligned}  l(\vec{\theta}) & = & \log \prod_{i=0}^{n-1} q_{\vec{\theta}}(y_i) \\  & \sim & - \sum_{i=0}^{n-1}\frac{1}{2}  \left( \frac{y_i -y_{\vec{\theta}}(x_i)}{\sigma_i}\right)^2\end{aligned}

Việc làm cực đại l(θ⃗) tương đương với làm cực tiểu  − 2l(θ⃗), do vậy ta làm cực tiểu hiệu số quân phương

\chi^2_{\vec{\theta}} = \sum_{i=0}^{n-1}  \left( \frac{y_i -y_{\vec{\theta}}(x_i)}{\sigma_i}\right)^2 (86)

Điều này nghĩa là các tham số θ⃗ được xác định sao cho hàm yθ⃗(x) đi theo các điểm số liệu {(x0, y0), …(xn − 1, yn − 1)} càng sát càng tốt, và độ lệch được đo bằng các thang sai số σi. Như vậy, những điểm số liệu có thang sai số nhỏ sẽ tham gia vào phương trình với trọng lượng lớn hơn. Toàn bộ quy trình này được gọi là khớp bình phương nhỏ nhất.

Hiệu số quân phương cần làm cực tiểu là biến ngẫu nhiên. Lưu ý rằng những số hạng không hề độc lập thống kê, vì chúng được liên hệ bằng n_{\rm p} tham số \vec{\hat \theta}, các tham số này được xác định bằng cách làm tối thiểu χθ⃗2. Do vậy, dạng phân bố của \chi^2_{\vec{\hat\theta}} được xấp xỉ bởi phân bố khi-bình phương (mà hàm mật độ xác suất cho bởi PT (45)) với nnp bậc tự do. Phân bố này có thể được dùng để đánh giá ý nghĩa thống kê của một phép khớp bình phương nhỏ nhất, như sau đây.

Trong tường hợp ta muốn mô phỏng hàm phân bố chưa biết của một mẫu, như ở Mục 6.1, chẳng hạn, cho một phân bố liên tục, thì về nguyên tắc cũng dùng được phương pháp bình phương nhỏ nhất. Ở đây ta sẽ khớp hàm mật độ xác suất được gắn tham số vào một histogram hàm mật độ xác suất, mà cũng theo dạng mẫu như trên {(xi, yi, σi)}. Mặc dù nguyên tắc bình phương nhỏ nhất được phát triển dựa trên nguyên tắc hợp lý cục đại, nhưng thường thì những tham số khác nhau sẽ nhận được nếu ta khớp một hàm mật độ xác suất vào với một histogram hàm mật độ xác suất so với việc nhận những tham số này trực tiếp từ phương pháp hợp lý cực đại. Thường thì phương pháp hợp lý cực đại sẽ cho kết quả chính xác hơn [Bauke (2007)]. Do đó ta nên dùng phép khớp bình phương nhỏ nhất cho việc khớp một hàm không phải là hàm khối xác suất hoặc mật độ xác suất lên một bộ số liệu.

Thật may là, để thực sự khớp bình phương nhỏ nhất, bạn không cần phải tự viết những hàm khớp, vì đã có những mã lệnh sẵn có. Cả hai chương trình ở Mục 4 là gnuplotxmgrace, đều cho phép khớp hàm bất kì. Bạn nên dùng gnuplot, vì nó linh hoạt hơn trong việc khớp hàm và cho bạn biết nhiều thông tin có ích hơn để ước tính chất lượng của phép khớp.

Lấy ví dụ, giả sử ta muốn khớp một hàm đại số theo dạng f(L) = e + aLb với bộ số liệu trong file sg_e0_L.dat. Trước hết, bạn phải định nghĩa hàm rồi cung cấp một số ước tính sơ bộ (khác 0) về các tham số chưa biết. Lưu ý rằng toán tử lũy thừa được kí hiệu bởi * và đối số tiêu chuẩn cho một lời định nghĩa hàm là x, nhưng cách đặt tên đối số này là tùy bạn chọn:

gnuplot> f(x)=e+a*x**b
gnuplot> e=-1.8
gnuplot> a=1
gnuplot> b=-1

Việc khớp được thực hiện nhờ lệnh fit. Chương trình dùng thuật toán bình phương nhỏ nhất phi tuyến của Levenberg-Marquardt [Press và nnk. (2007)], cách này cho phép ta khớp số liệu với gần như bất kì hàm nào. Để thực thi lệnh, bạn phải chỉ rõ hàm khớp, bộ số liệu và những tham số cần chỉnh. Trước hết, ta xét trường hợp chỉ có hai cột số liệu được dùng hoặc sẵn có (trong trường hợp này, gnuplot cho rằng σi = 1). Với ví dụ hiện tại, bạn có thể nhập vào:

gnuplot> fit f(x) "sg_e0_L.dat" via e,a,b

Tiếp theo gnuplot viết những thông tin ghi chép (log) vào bộ phận đầu ra, để mô tả quá trình khớp. Sau khi việc tính toán khớp đã hội tụ, chương trình sẽ in ra thông tin sau, đối với ví dụ đang xét:

After 17 iterations the fit converged.
final sum of squares of residuals : 7.55104e-06
rel. change during last iteration : -2.42172e-09

degrees of freedom (ndf) : 5
rms of residuals      (stdfit) = sqrt(WSSR/ndf)      : 0.00122891
phương sai of residuals (reduced chisquare) = WSSR/ndf : 1.51021e-06

Final set of parameters            Asymptotic Standard Error
=======================            ==========================

e               = -1.78786         +/- 0.0008548    (0.04781%)
a               = 2.5425           +/- 0.2282       (8.976%)
b               = -2.80103         +/- 0.08265      (2.951%)

correlation matrix of the fit parameters:

               e      a      b      
e               1.000 
a               0.708  1.000 
b              -0.766 -0.991  1.000 
 

Các dòng cần quan tâm nhất là những dòng tại đó \vec{\hat\theta} cho các tham số cùng với thang sai số chuẩn.4 Ngoài ra, chất lượng của phép khớp có thể được ước tính bằng thông tin chứa trong ba dòng bắt đầu bằng “degree of freedom”. Dòng thứ nhất trong số chúng cho biết số bậc tự do, vốn chỉ bằng nnp. Hiệu số quân phương \chi^2_{\vec{\hat\theta}} được kí hiệu là WSSR trong kết quả của gnuplot. Một độ đo chất lượng phép khớp là xác suất Q mà giá trị của hiệu số quân phương bằng hoặc lớn hơn so với giá trị của phép khớp hiện thời, với giả sử rằng các điểm số liệu được phân bố như trong PT (85) [Press và nnk. (1995)]. Giá trị của Q càng lớn, chất lượng của phép khớp càng cao. Như đã đề cập từ trước, Q cso thể được tính từ một phân bố khi-bình phương với nnp bậc tự do. Để tính Q bằng gnuplot output bạn có thể dùng chương trình nhỏ Q.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_cdf.h>

int main(int argc, char *argv[])
{
  double WSSRndf;
  int ndf;

  if(argc != 3)
  {
    printf("USAGE %s <ndf> <WSSR/ndf>\n", argv[0]);
    exit(1);
  }
  ndf = atoi(argv[1]);
  sscanf(argv[2], "%lf", &WSSRndf);
  printf("# Q=%e\n", gsl_cdf_chisq_Q(ndf*WSSRndf, ndf));

  return(0);
}

Chương trình này dùng hàm gsl_cdf_chisq_Q() từ GSL. Chương trình được gọi bằng lệnh theo mẫu Q <ndf> <WSSR/ndf>, vốn có thể được lấy được từ đầu ra của gnuplot. Lưu ý rằng trong trường hợp này ta thu được Q = 1, đây là giá trị quá lớn, vì σi = 1 đã được dùng, như dưới đây.

Để xem kết quả của phép khớp cùng với số liệu gốc, bạn hãy gõ vào

gnuplot> plot "sg_e0_L.dat" w e, f(x)

Hình 29. Cửa sổ Gnuplot cho thấy kết quả của một lệnh khớp cùng với số liệu đầu vào.

Kết quả được xuất hiện trên Hình 29. Lưu ý rằng sự hội tụ thì phụ thuộc vào lựa chọn giá trị ban đầu các tham số. Thuật toán có thể dẫn tới việc quá trình dò tìm bị giam vào một cực tiểu địa phương trong trường hợp những tham số cách quá xa nhữn giá trị tốt nhất. Hãy thử các giá trị ban đầu e=1, a=-3b=1 xem sao! Hơn nữa, không phải mọi tham số cho hàm đều bắt buộc phải tham gia quá trình khớp. Bạn cũng có thể ấn định kết quả cho một vài tham số nhất định đồng thời bỏ qua chúng trong danh sách via ở cuối lệnh fit. Hãy nhớ rằng trong ví dụ trên tất cả các điểm số liệu tham gia vào kết quả đều có đóng góp như nhau, nghĩa là coi như σi = 1 ∀ i. Bạn có thể ra lệnh cho thuật toán phải xét đến các thang sai số, chẳng hạn được xếp ở cột thứ ba, bằng cách gõ vào

gnuplot> fit f(x) "sg_e0_L.dat" using 1:2:3 via e,a,b

Như vậy, các điểm số liệu với các thang sai số lớn hơn thì ít ảnh hưởng đến số liệu. Trong trường hợp này, một kết quả khác với giá trị Q sẽ xuất hiện (Bạn hãy thử làm xem!).

Sau cùng, bạn có thể hạn chế các điểm số liệu được xét đến trong việc khớp, vốn chỉ áp dụng được nếu một tập con của mẫu tuân theo quy luật hàm số mà bạn xét. Điều này có thể được làm theo cách giống như việc hạn chế phạm vi các giá trị được vẽ biểu đồ, chẳng hạn bằng cách gõ vào

gnuplot> fit [5:12] f(x) "sg_e0_L.dat" using 1:2:3 via e,a,b

Các thông tin khác về cách dùng lệnh fit, chẳng hạn như khớp số liệu nhiều chiều, có thể tìm được bằng cách dùng trợ giúp trực tuyến của gnuplot qua việc gõ lệnh help fit.

Bài tập

(8) Khớp bình phương nhỏ nhất

Hãy sao chép chương trình từ ví dụ (2) thành một chương trình mới rồi sửa nó đi sao cho các giá trị của phân bố Fisher-Tippett dịch chuyển ngang, với các tham số λ cùng vị trí đỉnh x0, được phát sinh ra. Những số này phải được lưu vào trong một histogram và một hàm mật độ xác suất được ghi ra bộ phận đầu ra chuẩn.

  • Hãy chọn các tham số cho histogram (phạm vị, phạm vi ngăn) sao cho các histogram phù hợp với số liệu được phát sinh.
  • Chạy chương trình để phát sinh ra n = 105 giá trị ứng với các tham số x0 = 2.0λ = 3.0. Hãy dẫn kết quả histogram hàm mật độ xác suất tới một file (chẳng hạn bằng cách gõ thêm > ft.dat ở cuối câu lệnh).
  • Vẽ đồ thị kết quả bằng gnuplot.
  • Định nghĩa hàm mật độ xác suất cho phân bố Fisher-Tippett trong gnuplot rồi khớp hàm số này với số liệu, dùng x0λ làm các tham số chỉnh được. Hãy chọn một khoảng thích hợp để thực hiện việc khớp.
  • Vẽ đồ thị số liệu cùng với hàm khớp.
  • Kết quả so với phép khớp hợp lý cực đại được đề cập ở Mục 6.1 như thế nào?
  • Liệu phép khớp (đặc biệt là với λ) có được cải thiện không nếu bạn tăng số các điểm mẫu lên đến 106?

Gợi ý: Việc dịch chuyển hàm phân bố được thực hiện đơn giản bằng cách thêm x0 vào số ngẫu nhiên được phát sinh ra. Hãy dùng histogram từ Chương 3, hoặc lập nên một histogram tạm dùng bằng một mảng hist (đồng thời xem ở hàm main() của chương trình reject.c trình bày một phần ở Mục 2.4).

Lời giải: fitFT.gp fisher_tippett2.c

Tài liệu được trích dẫn

Bauke, H. (2007). Parameter estimation for power-law distributions by maximum likelihood methods, Eur. Phys. J. B 58, pp. 167-173.

Dekking, F. M., Kraaikamp, C., Lopuhaa, H. P., and Meester, L. E. (2005). A Modern Introduction to Probability and Statistics, (Springer, London).

Galassi M. et al (2006). GNU Scienti c Library Reference Manual, (Network Theory Ltd, Bristol), see also http://www.gnu.org/software/gsl/

Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P. (1995). Numerical Recipes in C, (Cambridge University Press, Cambridge).


  1. Đôi khi các điểm số liệu xi cũng là những đại lượng đo đạc đặc trưng bởi các thang sai số. Có thể dễ dàng khái quá hoá phương pháp đang xét cho trường hợp như vậy.
  2. Simplex là một tập lồi trong không gian n chiều phát sinh từ n + 1 đỉnh.
  3. Simplex này được giãn rộng bằng parn véc-tơ cho bởi par cộng với (0, …, 0,  simplex_size[i], 0, …, 0) với i = 1, …, n.
  4. Các “thang sai số” này được tính theo một cách mà chỉ đúng khi khớp các hàm tuyến tính, do vậy khi dùng chúng ta phải cẩn thận.
Advertisements

1 Phản hồi

Filed under Ngẫu nhiên và mô phỏng

One response to “6. Những ước lượng tổng quát

  1. Pingback: Ngẫu nhiên và xác suất trong mô phỏng máy tính | Blog của Chiến

Trả lời

Mời bạn điền thông tin vào ô dưới đây hoặc kích vào một biểu tượng để đăng nhập:

WordPress.com Logo

Bạn đang bình luận bằng tài khoản WordPress.com Log Out / Thay đổi )

Twitter picture

Bạn đang bình luận bằng tài khoản Twitter Log Out / Thay đổi )

Facebook photo

Bạn đang bình luận bằng tài khoản Facebook Log Out / Thay đổi )

Google+ photo

Bạn đang bình luận bằng tài khoản Google+ Log Out / Thay đổi )

Connecting to %s