## Copyright (C) 2010 Marius Hofert and Martin Maechler
##
## This program is free software; you can redistribute it and/or modify it under
## the terms of the GNU General Public License as published by the Free Software
## Foundation; either version 3 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
## FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
## details.
##
## You should have received a copy of the GNU General Public License along with
## this program; if not, see .
## CAUTION: sign.fafac is not correct [see below for a corrected version]
## Utility for polyG() in ./R/aux-acopula.R :
sign.fafac <- function(alpha, d) { ## also on (d ; k:= 1:d ; k1:= 0:(d-1))
stopifnot(length(alpha)==1, 0 < alpha, alpha < 1,
length(d) == 1, d == round(d), d >= 1)
k <- 1:d
k1 <- k - 1L
s <- unlist(lapply(alpha*k, function(z) prod(z-k1)))
ss <- sign(s)
sn <- (-1)^d * (2*(floor(alpha*k) %% 2) - 1)
## --------------------------------------
stopifnot(ss[s != 0] == sn[s != 0])
ss
}
sign.fafac(0.8, 4)
sign.fafac(0.8, 5)
sign.fafac(0.8, 6)
sign.fafac(0.8, 16)
sign.fafac(0.8, 17)
## determine signs of the falling factorials
## s <- unlist(lapply(alpha*k, function(z) prod(z-(0:(d-1)))))
## signs <- (-1)^(d-k) * sign(s) ## see ../misc/sign-polyG.R
## signs <- (-1)^(d-k) * (-1)^d * (2*(floor(alpha*k) %% 2) - 1)
## signs = (-1)^k * (2*(floor(alpha*k) %% 2) - 1)
## CAUTION: sign.fafac is not correct [see below for a corrected version]
## Ok: It is "proven"
##
## sign(s) == (-1)^d * (2*(floor(alpha*k) %% 2) - 1)
##
## at whenever sign(s) != 0 {and for the s == 0 case, sign does not matter}
## test 1
alpha <- 0.9
d <- 2
j <- 1:d
signFF(alpha, 1:d, )
sign(choose(alpha*j, d)*(-1)^(d-j))
(-1)^d * (2*(floor(alpha*j) %% 2) - 1) # => formula above is not correct!
## test 2
alpha <- 0.9
d <- 3
j <- 1:d
all(signFF(alpha, 1:d, d) == sign(choose(alpha*j, d)*(-1)^(d-j)))