Problems with applying Deriv() to dbinom(). about deriv

[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.