2012-05-27 29 views
6

Tôi muốn tích hợp hàm mật độ xác suất từ ​​(-\infty, a] do không có sẵn tệp cdf ở dạng đóng. Nhưng tôi không chắc chắn làm thế nào để làm điều này trong C + +.Thường trình bậc hai cho mật độ xác suất

Nhiệm vụ này khá đơn giản trong Mathematica; Tất cả những gì tôi cần làm là xác định hàm,

f[x_, lambda_, alpha_, beta_, mu_] := 
    Module[{gamma}, 
    gamma = Sqrt[alpha^2 - beta^2]; 
    (gamma^(2*lambda)/((2*alpha)^(lambda - 1/2)*Sqrt[Pi]*Gamma[lambda]))* 
     Abs[x - mu]^(lambda - 1/2)* 
     BesselK[lambda - 1/2, alpha Abs[x - mu]] E^(beta (x - mu)) 
    ]; 

và sau đó gọi số NIntegrate Định kỳ để tích hợp số.

F[x_, lambda_, alpha_, beta_, mu_] := 
    NIntegrate[f[t, lambda, alpha, beta, mu], {t, -\[Infinity], x}] 

Bây giờ tôi muốn đạt được điều tương tự trong C++. Tôi sử dụng thường lệ gsl_integration_qagil từ thư viện số gsl. Nó được thiết kế để tích hợp các chức năng trong khoảng thời gian bán vô hạn (-\infty, a] mà chỉ là những gì tôi muốn. Nhưng tiếc là tôi không thể làm cho nó hoạt động được.

Đây là hàm mật độ trong C++,

density(double x) 
{ 
using namespace boost::math; 

if(x == _mu) 
    return std::numeric_limits<double>::infinity(); 

    return pow(_gamma, 2*_lambda)/(pow(2*_alpha, _lambda-0.5)*sqrt(_pi)*tgamma(_lambda))* pow(abs(x-_mu), _lambda - 0.5) * cyl_bessel_k(_lambda-0.5, _alpha*abs(x - _mu)) * exp(_beta*(x - _mu)); 

} 

Sau đó, tôi cố gắng và tích hợp để có được những lũy ​​bằng cách gọi thông thường GSL.

cdf(double x) 
{ 
gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000); 

    double result, error;  
    gsl_function F; 
    F.function = &density; 

    double epsabs = 0; 
    double epsrel = 1e-12; 

    gsl_integration_qagil (&F, x, epsabs, epsrel, 1000, w, &result, &error); 

    printf("result   = % .18f\n", result); 
    printf ("estimated error = % .18f\n", error); 
    printf ("intervals = %d\n", w->size); 

    gsl_integration_workspace_free (w); 

    return result; 

} 

Tuy nhiên gsl_integration_qagil trả về lỗi, number of iterations was insufficient.

double mu = 0.0f; 
double lambda = 3.0f; 
double alpha = 265.0f; 
double beta = -5.0f; 

cout << cdf(0.01) << endl; 

Nếu tôi tăng kích thước của không gian làm việc thì chức năng bessel sẽ không đánh giá.

Tôi đã tự hỏi liệu có ai có thể cho tôi bất kỳ thông tin chi tiết nào về vấn đề của tôi hay không. Một cuộc gọi đến hàm Mathematica tương ứng F ở trên với x = 0.01 trả về 0.904384.

Có thể mật độ được tập trung xung quanh một khoảng thời gian rất nhỏ (ví dụ: bên ngoài [-0.05, 0.05] mật độ gần như là 0, một ô được đưa ra dưới đây). Nếu có thể làm gì về điều này. Cảm ơn vì đã đọc.

density

+0

Chức năng này có vẻ đối xứng, nghĩa là 'cdf (0) = 1/2'.Hãy nhớ rằng cdf được đánh giá tại x là giống như tích phân từ 0 đến x cộng với cdf được đánh giá ở 0. Tất nhiên, tôi đang đi từ hình dạng của biểu đồ, nó có thể không thực sự chính xác đối xứng. –

+0

Xin chào @Ben, tôi ước rằng đó là trường hợp! – mark

Trả lời

1

Re: tích hợp lên +/- vô cực:

Tôi sẽ sử dụng Mathematica để tìm một ràng buộc thực nghiệm cho | x - μ | >> K, trong đó K đại diện cho "chiều rộng" xung quanh giá trị trung bình và K là hàm alpha, beta và lambda - ví dụ F nhỏ hơn và xấp xỉ bằng (x- μ) -2 hoặc ae -b (x- μ) hoặc bất kỳ thứ gì. Các chức năng này đã biết tích phân ra đến vô cùng, mà bạn có thể đánh giá theo kinh nghiệm. Sau đó, bạn có thể tích hợp số lượng ra K, và sử dụng xấp xỉ bounded để có được từ K đến vô cùng.

Tìm ra K có thể hơi phức tạp; Tôi không quen thuộc với các hàm Bessel nên tôi không thể giúp bạn nhiều ở đó.

Nói chung, tôi nhận thấy rằng để tính toán số không rõ ràng, cách tốt nhất là thực hiện nhiều phép toán phân tích như bạn có thể trước bạn thực hiện đánh giá bằng số. (Loại giống như một camera lấy nét tự động - lấy nó gần nơi bạn muốn, sau đó để máy ảnh làm phần còn lại.)

+0

Xin chào @ Jason, vì tôi biết hành vi tiệm cận của quyền sở hữu, tôi có thể phải áp dụng cách tiếp cận này. Nhưng tôi chỉ không hiểu làm thế nào toán học đã làm việc ra tích phân. Có lẽ nó đang làm một cái gì đó tương tự đằng sau hậu trường. – mark

+0

Bạn có bản sao Bí quyết số không? Tôi nghĩ có lẽ họ thảo luận về các tích phân mở rộng đến vô cùng, không chắc chắn cách họ xử lý nó. –

+0

Vuốt qua nó ngay bây giờ! – mark

1

Tôi đã không cố gắng ++ mã C, nhưng bằng cách xem các chức năng trong Mathematica, nó dường như vô cùng lên đến đỉnh điểm xung quanh mu, với sự lây lan của các đỉnh xác định bởi các thông số lambda, alpha, beta .

Điều tôi sẽ làm là tìm kiếm sơ bộ pdf: xem bên phải và bên trái của x = mu cho đến khi bạn tìm thấy giá trị đầu tiên bên dưới dung sai đã cho. Sử dụng chúng như là giới hạn cho cdf của bạn, thay vì vô cực âm.

Pseudo code sau:

x_mu 
step = 0.000001 
adaptive_step(y_value) -> returns a small step size if close to 0, and larger if far. 

while (pdf_current > tolerance): 
    step = adaptive_step(pdf_current) 
    xtest = xtest - step 
    pdf_current = pdf(xtest) 

left_bound = xtest 

//repeat for left bound 

Với lên đến đỉnh điểm cách chặt chẽ chức năng này có vẻ là, thắt chặt các giới hạn có lẽ sẽ giúp bạn tiết kiệm rất nhiều thời gian máy tính mà hiện đang bị lãng phí tính toán số không. Ngoài ra, bạn có thể sử dụng thói quen tích hợp bị chặn, thay vì - \ infty, b.

Chỉ cần một ý nghĩ ...

PS: Mathematica mang lại cho tôi F [0.01, 3, 265, -5, 0] = 0.884505

+0

Xin chào @Adriano, Cảm ơn bạn đã nhập. Tôi có thể phải tiếp cận vấn đề của mình theo kiểu này. – mark

1

Tôi tìm thấy mô tả đầy đủ về glsl này có http://linux.math.tifr.res.in/manuals/html/gsl-ref-html/gsl-ref_16.html, bạn có thể tìm thấy thông tin hữu ích.

Vì tôi không phải là chuyên gia GSL, tôi không tập trung vào vấn đề của bạn từ quan điểm toán học, mà đúng hơn là tôi đã nhắc bạn một số khía cạnh quan trọng về lập trình điểm nổi.

Bạn không thể đại diện chính xác các số sử dụng chuẩn IEEE 754. MathLab làm ẩn thực tế bằng cách sử dụng một logic biểu diễn số vô hạn, để cung cấp cho bạn các kết quả không có lỗi, đây là lý do tại sao nó chậm so với mã gốc.

Tôi mạnh mẽ recommand liên kết này cho bất cứ ai tham gia vào tính toán khoa học sử dụng một FPU: http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html

Giả sử bạn đã rất thích bài viết đó, tôi đã nhận thấy điều này vào liên kết GSL trên: "Các thói quen sẽ không hội tụ nếu giới hạn lỗi quá nghiêm ngặt ".

Giới hạn của bạn có thể quá nghiêm ngặt nếu chênh lệch giữa phần trên và dưới nhỏ hơn giá trị tối thiểu có thể biểu diễn của số double, đó là std :: numeric_limits :: epsilon() ;. Ngoài ra, hãy nhớ rằng, từ liên kết thứ 2, đối với bất kỳ trình biên dịch C/C++ nào, chế độ làm tròn mặc định là "cắt ngắn", điều này giới thiệu các lỗi tính toán tinh tế leeding đến các kết quả sai. Tôi đã có vấn đề với một clipper dòng Liang Barsky đơn giản, thứ tự đầu tiên! Hãy tưởng tượng sự lộn xộn trong dòng này:

return pow(_gamma, 2*_lambda)/(pow(2*_alpha, _lambda-0.5)*sqrt(_pi)*tgamma(_lambda))* pow(abs(x-_mu), _lambda - 0.5) * cyl_bessel_k(_lambda-0.5, _alpha*abs(x - _mu)) * exp(_beta*(x - _mu)); 

Theo nguyên tắc chung, nó là khôn ngoan trong C/C++, để thêm biến thêm giữ kết quả trung gian, vì vậy bạn có thể gỡ lỗi từng bước, sau đó nhìn thấy bất kỳ lỗi làm tròn, bạn không nên cố gắng nhập biểu thức như thế này trong bất kỳ ngôn ngữ lập trình gốc nào. Người ta không thể tối ưu hóa các biến tốt hơn so với một trình biên dịch.

Cuối cùng, như một quy tắc chung, bạn nên nhân mọi thứ rồi chia, trừ khi bạn tự tin về hành vi động của phép tính của mình.

Chúc may mắn.

+0

Xin chào @ Xin chào, tôi ghét lỗi tắt! Để diễn giải Moore, "sự khắt khe bị mất khi chúng ta cần nó nhất". Ông đề cập đến các nhà toán học thực hiện chúng là các thuật toán sử dụng số học dấu chấm động. Do đó sự ra đời của số học khoảng. Tôi cho một cái nhìn về phía trước để ngày này trở thành một tính năng phần cứng và chúng tôi có thể nói tạm biệt để làm tròn! (cũng không hoàn toàn nhưng ít nhất có sự hiểu biết tốt hơn về nó) Tôi biết điều đó trước đây, nhưng tôi đã thúc đẩy tôi đọc nó! Cảm ơn bạn đã nhập số – mark

Các vấn đề liên quan