Komputasi fungsi hypergeometrik di R

8

Saya mengalami kesulitan luar biasa mengevaluasi dengan paket di R. Dalam kasus saya, nilai , , selalu bilangan real positif. Meski begitu, fungsi hypergeometrik sangat sensitif terhadap nilai-nilai mereka. Saya tidak mencari presisi ekstrim; Saya dapat menggunakan Excel untuk mendapatkan estimasi kasar dari hypergeometrik Guass yang baik untuk tujuan saya.2F1(a,b;c;z)hypergeoabc

Adakah saran untuk implementasi dalam R yang akan memberikan perhitungan hypergeometrik Gaussian yang cepat, bebas kesalahan, jika tidak super akurat dari bilangan real positif dengan berbagai nilai?

Sunting: sepertinya ada jauh lebih banyak kode untuk ini di MATLAB daripada R. Ada pemikiran mengapa?

Benrolls
sumber
Saya akan berpikir pendekatan yang baik adalah dengan menemukan istilah terbesar dalam jumlah hypergeometrik dan memotong sekitar istilah ini. dengan cara ini Anda dapat memperhitungkan faktor terbesar dan bekerja dengan istilah kurang dari 1 dalam jumlah. Anda juga bisa menggunakan metode laplace pada representasi integral (sesuai parameterised untuk menghindari masalah batas).
probabilityislogic

Jawaban:

14

Kecuali jika Anda perlu mengevaluasi fungsi hypergeometrik Gauss untuk nilai-nilai kompleks dari parameter atau variabel, lebih baik menggunakan gslpaket Robin Hankin .

Berdasarkan pengalaman saya, saya juga merekomendasikan untuk hanya mengevaluasi fungsi hypergeometrik Gauss untuk nilai variabel yang terletak pada , dan menggunakan rumus transformasi untuk nilai dalam .[0,1]],0]

library(gsl)
Gauss2F1 <- function(a,b,c,x){
    if(x>=0 & x<1){
        hyperg_2F1(a,b,c,x)
    }else{
            hyperg_2F1(c-a,b,c,1-1/(1-x))/(1-x)^b
        }
}

Memperbarui

Berikut ini adalah implementasi alternatif saya dengan paket gmp (setidaknya, untuk bersenang-senang)

Stéphane Laurent
sumber
1
Terima kasih atas solusi sederhana dan bersih untuk masalah yang tidak terlalu sederhana! Tidak ada banyak dokumentasi tentang fungsi hypergeometric gauss untuk R sehingga pendekatan ini sangat berharga.
Benrolls
@ Benrolls Saya telah mengalami fungsi hypergeometrik Gauss dengan paket hypergeo, paket gsl, metode lain yang diterapkan oleh saya sendiri, dan Mathematica juga. Saya menjamin bahwa solusi yang saya berikan di atas sangat performan, saya tidak pernah menemukan bug saat menggunakannya.
Stéphane Laurent
Formula Stéphane di atas bagus, tetapi saya perhatikan bahwa ia menghasilkan NaNs ketika negatif dengan nilai absolut yang besar. Penggunaan rumus transformasi lain mengarah ke: library (gsl)ca
pglpm
1

@ Stéphane Formula Laurent di atas sangat bagus. Aku telah memperhatikan bahwa kadang-kadang menghasilkan NaNs ketika a, b, cbesar dan znegatif - saya belum mampu pin kondisi yang tepat ke bawah. Dalam kasus ini kita dapat menggunakan transformasi hypergeometrik lain mulai dari ekspresi alternatif Stéphane. Ini mengarah ke formula alternatif ini:

library(gsl)
Gauss2F1b <- function(a,b,c,x){
    if(x>=0 & x<1){
        hyperg_2F1(a,b,c,x)
    }else{
        hyperg_2F1(a,c-b,c,1-1/(1-x))/(1-x)^a
    }
}

Sebagai contoh:

> Gauss2F1(80.2,50.1,61.3,-1)
[1] NaN
>
> Gauss2F1b(80.2,50.1,61.3,-1)
[1] 5.498597e-20
>
>
> Gauss2F1(80.2,50.1,61.3,-3)
[1] NaN
> Gauss2F1b(80.2,50.1,61.3,-3)
[1] 5.343807e-38
>
>
> Gauss2F1(80.2,50.1,61.3,-0.4)
[1] NaN
> Gauss2F1b(80.2,50.1,61.3,-0.4)
[1] 3.322785e-10

ketiganya setuju dengan Mathematica Hypergeometric2F1. Formula ini tampaknya berperilaku baik juga untuk lebih kecil a, b, c. Perhatikan bahwa ada kasus-kasus di mana formula ini memberi NaNdan Stéphane tidak . Terbaik untuk memeriksa kasus per kasus.

pglpm
sumber