[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Overhaul of statistical distribution functions
From: |
Michael D Godfrey |
Subject: |
Re: Overhaul of statistical distribution functions |
Date: |
Wed, 21 Sep 2011 11:17:54 -0700 |
User-agent: |
Mozilla/5.0 (X11; Linux x86_64; rv:6.0.2) Gecko/20110906 Thunderbird/6.0.2 |
On 09/21/2011 10:35 AM, Rik wrote:
Michael,
I cleaned up a lot, but wasn't always trying to re-implement
algorithms which looked reasonably solid. The Loader paper looks good
though, and particularly for large N as you said. Are there any
copyright issues with filing a bug report on the Octave bug tracker
and attaching the Loader paper?
I am sure not. The paper is freely available and the algorithm from it
is widely used, including, for
example, in R which is also free software. The algorithm is pretty
simple, but n as a vector may
need an extra line or two.
I just don't have time to implement it myself.
In this regard, there are still several .m functions which use
iterative algorithms which are going to be slow in Octave. These are:
binoinv
nbininv
poissinv
betainv
gaminv
hygecdf
hygeinv
hygernd
The first three are problematic because they use a 'while (1)'
construction and the stopping condition is not always met. For example,
poissinv(1-eps,4)
creates an infinite loop that requires a Ctrl-C to kill. At the very
least, it would be good to change these to for loops so that they
would stop eventually.
I would rather rewrite those entirely to avoid the iterative nature,
but to do so I need a function to compute the inverse of the
incomplete Beta function. Matlab has such a function called
betaincinv
(http://www.mathworks.com/help/techdoc/ref/betaincinv.html). Do you
have reference papers or code for such a function? A little bit of
searching found this
http://people.sc.fsu.edu/~jburkardt/f_src/asa109/asa109.html which has
implementations in C++, Fortran, or Matlab. Given the .edu address it
may be possible to contact the author and get agreement to use the
code. I'm out of time to pursue that angle though which is why I was
going to file an issue report about those functions.
--Rik
The AS (from the Journal Applied Statistics) algorithms are published
under the Royal Statistical Society. Their
copyright statement reads:
The Royal Statistical Society holds the copyright to these routines,
but has given its permission for their distribution provided that no
fee is charged.
Everyone that I know accepts this as permission to use for
non-commercial purposes. The algorithms
appear in many free software distributions. I have never seen any
specific authorization from the RSS attached
to any of the algorithms. Also, these algorithms are generally quite
reliable and are very widely used.
Quite a few have been translated to *.m files and therefore should
usable with minor updates. The best
source for these is: http://lib.stat.cmu.edu/apstat/ There is no need
to ask permission to use within the
non-commercial restriction. Since the RSS claims the copyright it would
be to them that any requests should
be directed.
In any case, it would make sense to look at the Statlib algorithms
before doing much more work. The
more of these algorithms are adopted in Octave the better.
About the binopdf problem: do you know for what values of n the Mac OS
version fails? If it is
above n=10, then the solution in any case is to implement the Loader
algorithm. The current code
should never be used for n > 10.
Anyhow, thanks for getting this this far.
Michael
- Re: tolerance in binopdf.m, (continued)
- Re: tolerance in binopdf.m, Ben Abbott, 2011/09/21
- Re: tolerance in binopdf.m, Marco atzeri, 2011/09/21
- Re: tolerance in binopdf.m, Ben Abbott, 2011/09/21
- Re: tolerance in binopdf.m, Jordi GutiƩrrez Hermoso, 2011/09/21
- Re: tolerance in binopdf.m, Dr. Alexander Klein, 2011/09/21
Re: Overhaul of statistical distribution functions, Michael D Godfrey, 2011/09/21
Re: Overhaul of statistical distribution functions, Michael D Godfrey, 2011/09/21