Theprob pakiet numerycznie ocenia funkcje charakterystyczne dla bazowych rozkładów R. Dla prawie wszystkich dystrybucji istnieją istniejące formuły. Jednak w kilku przypadkach nie jest znane rozwiązanie zamknięte. Przypadek w punkcie: rozkład Weibulla (ale patrz poniżej).

Dla funkcji charakterystycznej Weibulla zasadniczo obliczam dwie całki i łączę je ze sobą:

fr <- function(x) cos(t * x) * dweibull(x, shape, scale)
fi <- function(x) sin(t * x) * dweibull(x, shape, scale)
Rp <- integrate(fr, lower = 0, upper = Inf)$value
Ip <- integrate(fi, lower = 0, upper = Inf)$value
Rp + (0+1i) * Ip

Tak, jest niezgrabny, ale działa zaskakująco dobrze! ...ahem, większość czasu. Użytkownik zgłosił niedawno, że następujące przerwy:

cfweibull(56, shape = 0.5, scale = 1)

Error in integrate(fr, lower = 0, upper = Inf) : 
  the integral is probably divergent

Teraz wiemy, że całka nie jest rozbieżna, więc musi być problemem numerycznym. Z jakimś skrzypkiem mogłem wykonać następujące czynności:

fr <- function(x) cos(56 * x) * dweibull(x, 0.5, 1)

integrate(fr, lower = 0.00001, upper = Inf, subdivisions=1e7)$value
[1] 0.08024055

To jest OK, ale nie jest to w porządku, a to wymaga sporej ilości bzdur, które nie skalują się dobrze. Zbadałem to pod kątem lepszego rozwiązania. Znalazłem niedawno opublikowaną „formę zamkniętą” dla funkcji charakterystycznej za pomocąscale > 1 (Spójrz tutaj), ale to dotyczyUogólniona zlewająca się funkcja hipergeometryczna Wrighta który nie jest jeszcze zaimplementowany w R (jeszcze). Szukałem w archiwachintegrate alternatywy, a jest tam mnóstwo rzeczy, które nie wydają się zbyt dobrze zorganizowane.

W ramach tego poszukiwania przyszło mi do głowy przetłumaczyć obszar integracji na skończony przedział przez odwrotną styczną, avoila! Sprawdź to:

cfweibull3 <- function (t, shape, scale = 1){
  if (shape <= 0 || scale <= 0) 
    stop("shape and scale must be positive")
  fr <- function(x) cos(t * tan(x)) * dweibull(tan(x), shape, scale)/(cos(x))^2
  fi <- function(x) sin(t * tan(x)) * dweibull(tan(x), shape, scale)/(cos(x))^2
  Rp <- integrate(fr, lower = 0, upper = pi/2, stop.on.error = FALSE)$value
  Ip <- integrate(fi, lower = 0, upper = pi/2, stop.on.error = FALSE)$value
  Rp + (0+1i) * Ip
}

> cfweibull3(56, shape=0.5, scale = 1)
[1] 0.08297194+0.07528834i

Pytania:Czy możesz zrobić lepiej niż to?Czy jest coś w numerycznych procedurach integracyjnych, które ludzie, którzy są ekspertami w takich sprawach, mogą rzucić nieco światła na to, co się tutaj dzieje? Mam podejrzenie, że dla dużycht cosinus szybko się zmienia, co powoduje problemy ...?Czy istnieją istniejące procedury / pakiety R, które są lepiej dostosowane do tego typu problemów i czy ktoś mógłby wskazać mi dobrze umiejscowioną pozycję (na górze), aby rozpocząć wspinaczkę?Komentarze:Tak, to zła praktykat jako argument funkcji.Obliczyłem dokładną odpowiedź nashape > 1 używając opublikowanego wyniku z Maple ibrute-force-integrate-by-the-definition-with-R kopnął tyłek Maple'a. To znaczy, otrzymuję tę samą odpowiedź (aż do precyzji numerycznej) w małym ułamku sekundy i jeszcze mniejszym ułamku ceny.Edytować:

Zamierzałem zapisać dokładne całki, których szukam, ale wygląda na to, że ta konkretna witryna nie obsługuje MathJAX, więc zamiast tego dam linki. Szukam liczbowej ocenyfunkcja charakterystyczna zRozkład Weibulla za rozsądne nakładyt (cokolwiek to znaczy). Wartość jest liczbą złożoną, ale możemy podzielić ją na części rzeczywiste i wyimaginowane i to właśnie dzwoniłemRp iIp powyżej.

Ostatnia uwaga: Wikipedia ma wyszczególnioną formułę (seria nieskończona) dla Weibulla c.f. i ta formuła pasuje do tej, którą udowodniono w artykule, o którym mowa powyżej,jednak, ta seria została tylko udowodnionashape > 1. Walizka0 < shape < 1 jest nadal otwartym problemem; szczegóły na ten temat.