piLaboratory/sads

qbs does not allow non-integer total abundance

piklprado opened this issue · 4 comments

The Broken-stick is a continuous distribution but the quantile function returns NaN if we use non integer N:

> y <- seq(0.1,0.9, by=0.1)
> qbs(y, N=1000, S=100)
[1]  1.066649  2.260756  3.608093  5.152799  6.978123  9.221124 12.091714
[8] 16.131053 22.990529
> qbs(y, N=1000.1, S=100)
[1] NaN NaN NaN NaN NaN NaN NaN NaN NaN

We have an old PR to make BS a continuous distribution (#44), but the check for integer N was kept. Had we only missed that or was there any good reason to keep this condition in the function?

Fixed by 038b8bf

> y <- seq(0.1,0.9, by=0.1)
> qbs(y, N=1000, S=100)
[1]  1.066649  2.260756  3.608093  5.152799  6.978123  9.221124 12.091714
[8] 16.131053 22.990529
> qbs(y, N=1000.1, S=100)
[1]  1.066759  2.260985  3.608448  5.153331  6.978788  9.222071 12.092972
[8] 16.132725 22.992716

Still, a small discrepancy persists between quantile and cdf:

> pbs(qbs(y, N=100, S=10), N=100.5, S=10)
[1] 0.09998271 0.19984691 0.29900369 0.39903102 0.49871400 0.59845752 0.69842370
[8] 0.79848259 0.89883420

Something to be concerned about?

I noticed the problem with 'qbs' trying to plot a Broken-stick model fitted to continuous data. Returned an error because the 2nd plot is needs radpred which depends on the corresponding quantile function.

I have always thought that N and S should be constrained to be integers.

On the issue of pbs(qbs()), of course there are some approximation errors. These errors happen at all the functions in which we use the appproxfun to approximate quantiles (bs, gs, mand, mzsm...)

In this case, however, the errors are below 0.5%, independent of the N being integer or not:

> (y <- seq(0.1, 0.9, length.out=10))
 [1] 0.1000000 0.1888889 0.2777778 0.3666667 0.4555556
 [6] 0.5444444 0.6333333 0.7222222 0.8111111 0.9000000
> pbs(qbs(y, N=100, S=10), N = 100, S= 10)
 [1] 0.1004594 0.1895345 0.2784676 0.3667881 0.4561131
 [6] 0.5448992 0.6337408 0.7224840 0.8111949 0.9001463
> (pbs(qbs(y, N=100, S=10), N = 100, S= 10) - y)/y
 [1] 0.0045942513 0.0034178555 0.0024832532 0.0003311258
 [5] 0.0012238945 0.0008352793 0.0006433581 0.0003625235
 [9] 0.0001033395 0.0001625313
> max((pbs(qbs(y, N=100, S=10), N = 100, S= 10) - y)/y)
[1] 0.004594251
> max((pbs(qbs(y, N=100.5, S=10), N = 100.5, S= 10) - y)/y)
[1] 0.004676224

Your changes have fixed the problem of qbs not accepting non-integer N; does the same apply for qrbs?

# in qrbs.R
if (N <= 0 | S <= 0 | !is.wholenumber(N) | !is.wholenumber(S))  return(rep(NaN, length(p)))

I think it does not, because the rad for Broken-Stick is a discrete distribution.