Problems with applying Deriv() to dbinom().

[rolfTurner]

(1) I believe the code in the rule for differentiating dbinom(), found
in the drule environment, is incorrect. The attached script "demo01.txt"
provides evidence for my belief. (The derivative produced by Deriv()
is always negative, whereas it should be positive for prob < y/size
and negative for prob > y/size.)

The script demo01.txt also contains code for a proposed correction to the
rule.

Note in addition that "demo01.txt" raises the issue that if a rule is placed
in a new environment, rather than being placed in "drule", then an
error is thrown when one attempts to calculate a second derivative.
(Calculation of the first derivative seems to be OK.) This phenomenon
is illustrated in the attached script "demo02.txt".

(2) Incorrect results seem to be produced when second derivatives are
calculated when the derivative depends on a rule in "drule". This is
illustrated in the attached script "scr01.txt" This script also produces
correct values using calculations of first and second derivative done
"by hand". Note that the function value f1 agrees with the correct
function value f0, as does the first derivative df1 with the correct
value df0. However the second derivative value d2f1 is 1.769472 which
differs from the correct value, d2f0 = 2.064384. When derivatives are
calculated without applying a rule from "drule", as in the attached
script "scr02.txt", d2f1 agrees with the correct value.

Note that the problem is illustrated using a "local" version of dbinom()
called ldb() which has been simplified so as to clarify the issue. The
rule that I placed in "drule" for differentiating ldb() is also made
to be very simplistic so as to illustrate the problem clearly. This
rule is obviously not "robust" and would result in problems with
borderline cases.

demo01.txt
demo02.txt
scr01.txt
scr02.txt

[sgsokol]

Thanks for reporting. I confirm both problems that I could reproduce on my side.
I will see how it could pass automatic tests which compare with numerical differences.

[sgsokol]

dbinom's rule is fixed in v4.1.4. Could you confirm that it works correctly on your side, please?

devtools::install_github("sgsokol/Deriv")
library(Deriv)
packageVersion("Deriv")
[1] ‘4.1.4’

[rolfTurner]

I sourced my dem01 script and the graph of the derivative of dbinom() now looks correct.

I also executed

library(Deriv)
packageVersion("Deriv")
[1] ‘4.1.4’
D2 <- Deriv(dbinom,"prob",nderiv=2)
D2(3,8,0.6)

and got 2.064384 which is the correct answer.

However I still get an error when I try to define a differentiation rule in a new environment rather than in "drule". E.g.:
tempEnv <- new.env()
tempEnv[["dbinom"]] <- alist(x=NULL,size=NULL,prob=
(if (x == 0) -size * (1 - prob)^(size - 1) else if (x == size) size *
prob^(size - 1) else dbinom(x, size, prob) * (x - prob *
size)/(prob - prob * prob))/(if (log) dbinom(x, size, prob,
log = FALSE) else 1))
xxx <- Deriv(dbinom,"prob",drule=tempEnv,nderiv=2)
Error in FUN(X[[i]], ...) : Could not retrieve body of '^()'
First derivatives seem to work OK..
Am I misunderstanding something about the syntax for using a new environmnt to contain differentiation rules?

I also remain puzzled by the fact that things still don't work with my somewhat shakey "local" version of dbinom(). If I execute

ldb <- function(x,size,prob) {
    choose(size,x)*prob^x*(1-prob)^(size-x)
}
drule[["ldb"]] <- alist(prob={
                             .e1 <- 1 - prob
                             .e2 <- size - 1
                             size*(ldb(x-1,.e2,prob) - ldb(x,.e2,prob))
                             })
lD2 <- Deriv(ldb,"prob",nderiv=2)

then I get lD2(3,8,0.6) equal to 1.769472, the "old" wrong answer.
Am I doing something silly here? I can't see what that might be.

[sgsokol]

One thing after another. I did not yet get a look at the problem on a new environment. But it is on my "todo" list.
As to your "silly" problem of ldb(), it is not you but Deriv which is badly functioning. There is kind of collision between ".e2" that you use in drule[["ldb"]] and ".e2" automatically generated and used by Deriv. If you rewrite the rule as

drule[["ldb"]] <- alist(prob=size*(ldb(x-1,size-1,prob) - ldb(x,size-1,prob)))

things work as expected.

[rolfTurner]

Ah!!! The light dawns. I really was being silly though --- it's now clear in
retrospect. (Hindsight is 20-20, and all that.). I put together the code for my rule for "ldb" by looking at Deriv(dbinom,"prob"). What I should have looked at of course was drule[["dbinom"]]. Had I done so (and noticed the obvious difference!) I would not have made my silly mistake.

Perhaps the help file for Deriv should warn the user NOT to use variables named ".e1", ";e2", ... etc. when constructing new rules.

Good luck with fixing the problem with putting rules in new environments.

[rolfTurner]

Here is another annoyance/item to add to your already overfull to-do list.

I have found that when the "x" argument, provided to the derivative of dbinom() ,is a vector, disconcerting warning messages result. E.g.:

d1 <- Deriv(dbinom,"prob")
x <- 3:5
d1(x,8,0.6)
[1] -0.9289728 -0.7741440 0.2322432
Warning messages:
1: In if (x == 0) -(size * .e1^.e2) else if (x == size) prob^.e2 * :
the condition has length > 1 and only the first element will be used
2: In if (x == size) prob^.e2 * size else dbinom(x, size, prob) * (x - :
the condition has length > 1 and only the first element will be used

I think that this can be fixed by changing the rule for dbinom() to something like the following:

prob=(ifelse(x %in% 0, -size * (1 - prob)^(size - 1),
ifelse(x %in% size, size * prob^(size - 1),
dbinom(x, size, prob) * (x - prob * size)/(prob - prob * prob)))/
(if (log) dbinom(x, size, prob,log = FALSE) else 1))

You think?

[sgsokol]

Even if Deriv() pretends to work with scalar expressions, you are right, it's better to make it as much compatible as possible with vector expressions. That will require to fix many other similar expressions.

[sgsokol]

Perhaps the help file for Deriv should warn the user NOT to use variables named ".e1", ";e2", ... etc. when constructing new rules.

The situation was revealed to be worthier than that. Even if you use tmp instead of .e2 (or whatever similar), it won't work. The problem is the incorrect reuse of a variable in generated subexpressions.

[sgsokol]

current version fixes both new.env() and compound expressions in drule. Does it work well on your side ?

[rolfTurner]

The "new.env()" problem seems to be solved:

tempEnv <- new.env()
tempEnv[["dpois"]] <- alist(lambda= dpois(x,lambda)*(x/lambda - 1))
xxx <- Deriv(dpois,"lambda",drule=tempEnv,nderiv=2)

No error thrown, and the answer looks correct to me.

However there seems to be a new (?) problem with providing a vector argument to the output of Deriv():

xxx <- Deriv(dbinom,"prob")
xxx(3,8,0.6) # -0.9289728 OK
xxx(3,8,0.5) # -0.875 OK
xxx(3,8,c(0.5,0.6)) # -0.875, i.e. only the first entry of "prob" is used.

Sorry 'bout that.

[sgsokol]

Right. It's because ifelse() always returns a result of length of the test. Maybe it's better to revert to if() and let the result be dependent on scalars x and size, and vector prob?

[rolfTurner]

Arrghhh! I see the problem. Thanks for working out what was going wrong.

In the application that I have in mind, I really need x to be allowed to be a vector quantity. I thought just now that I had worked out a rule that would accommodate this desideratum (not using ifelse()). However I cannot get my idea to work. Also, this idea involves using lists and unlist() so even if I could get it working for first derivatives, I suspect that it would break for second and higher derivatives.

So, yes. Go back to the if / else construction and eschew ifelse(). I think that I will be able to handle things in my application by using sapply() on vector valued x.

Sorry for introducing this red herring.