-
Notifications
You must be signed in to change notification settings - Fork 40
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Feature request: Owen's t function #99
Comments
@maximerischard @johnczito |
Oh, what do you know. Thanks for the heads up. If that gets merged and released, we can close this and you can open the PR in |
Looks like JuliaMath/SpecialFunctions.jl#242 flamed out due to licensing issues -- the code it was derived from was LGPL. I can work on this next, basing it on Fast and Accurate Calculation of Owen's T Function, which won't have a license issue (The Journal of Statistical Software uses the Creative Commons License). I'm working on getting the MVN cdf (qsmivnv.jl) over the finish line, and the Johnson Distributions for Distributions.jl right now. Next I can turn my attention to Owen's T. |
@blackeneth that would be awesome! |
@azev77 and others ... History Six algorithms were given, and which one used depended on the values of (h,a) in owent(h,a)
They developed code for these algorithms on a DEC VAX 750. The VAX 750 came out in 1980 and had a processor clock speed of 3.125 MHz. The reason for 6 algorithms was to speed up the function when possible, with T1 being faster than T2, T2 faster than T3, etc. TODAY We have faster and better computers today. If we skip T1 through T4, we can go straight to 48 point Gauss-Legendre quadrature for a <= 0.999999. For a > 0.999999, we can fall back to quadrature using QuadGK. This is sufficient for 15 digit accuracy for Float64. Gauss-Legendre quadrature takes about 3 microseconds, and QuadGK takes about 20 microseconds. Questions:
Code I call this function "owensteasy(h,a)":
You can test it with this code:
Here are the results I get from the above testing:
|
@johnczito @andreasnoack @mschauer |
Whether using generic quadrature schemes is "sufficient" depends on the use case. It is typically orders of magnitude slower than specialized expansions like continued fractions, so you will probably be losing compared to other languages. That being said, I agree that the optimal choice of algorithm these days is probably different than that old paper. |
Does it make sense to discuss this in the larger context of the package? Are our other implementations typically "as fast as reasonably can be" or is there a continuous work on adding correct but perhaps not optimal implementations and later performance improvements? |
@mschauer I agree. |
True, but it's also nice to avoid "embarrassingly slow compared to R or Python", and it would be good to check whether the quadrature approach falls into this category. |
As discussed in #114, we should not add a circular dependency on Distributions to StatsFuns. It's seems it is not actually needed and you could just call I am also worried that a dependency on QuadGK will make the package slower to load. In the past this has been a reason to not merge PRs (even though I think #106 should be fine 😛 - it only adds a lightweight dependency and IMO the timings are increased only to a small extent). |
@blackeneth can you compare your functions performance w implementations in other languages? |
I can replace the Distributions.jl call to Normal() with As for R & Python, I can't say I'm an expert at benchmarking those, but my initial checks have R at 4.74 μs and Python at 2.5 μs I will look at using "T6" for 0.999999 < a < 1.0, which could eliminate the QuadGK dependency. |
That's precisely what I was referring to since StatsFuns.jl/src/distrs/norm.jl Line 45 in 36e48bb
|
Using method "T6" for a > 0.999999 has similar accuracy to using QuadGK:
The performance is quite good for 0.999999 < a < 1.0 Some additional checks for type stability, then a few more tests, and it should be good for a PR. |
The function has been ready for some time. The barrier is git. My VS Code is somehow setup wrong, so it wants to release the wrong files. I rarely use git, so I forget how to use it after every use. But now I'm moving, so won't have time to figure it out and release it for quite a while. Attached is the code for the function, owent(h,a). Also attached is the code for the tests. If someone wants to take them and release them, please do so. |
Owen's T function currently does not seem to have a native julia implementation. Most existing implementations in Fortran or C are under GPL, which makes this a little tricky.
This would be useful for the Skew Normal distribution in Distributions.jl.
The text was updated successfully, but these errors were encountered: