Блог

У меня возникли две связанные проблемы с функцией Deriv() из пакета Deriv (из CRAN).

Проблема (1): Я считаю, что код в правиле для дифференцирования dbinom() w.r.t «prob» неверен. Доказательства этого (и мое предлагаемое исправление) приведены в следующем коде. Я бы предпочел прикрепить текстовые файлы, содержащие код, но, насколько я вижу, это невозможно сделать.

 # # Script demo01.R. #  library(Deriv) # Plot dbinom as a function of probabiity. plot(function(p){dbinom(3,8,p)},from=0,to=1,xlab="parameter "prob"",  ylab="binomial probability",main="dbinom") abline(v=3/8,col="red") readline("Go? ")  # Plot the derivative of dbinom, with respect to prob, calculated by # Deriv(), as a function of probability. Ddb lt;- Deriv(dbinom,"prob") plot(function(p){Ddb(3,8,p)},from=0,to=1,xlab="parameter "prob"",  ylab="",main="derivative of dbinom") readline("Go? ")  # Replace what I believe to be incorrect code for the rule for # differentiating dbinom() with what I believe to be correct code. # This rule should, strictly speaking, be placed in a new environment # rather than over-writing the existing rule, but this seems to # break down when second derivatives are taken. drule[["dbinom"]] lt;- alist(x=NULL,size=NULL,prob={  .e1 lt;- 1 - prob  .e2 lt;- size - 1  if (x == 0)  -(x * .e1^(x - 1))  else if (x == size)  prob^.e2 * size  else size*(dbinom(x-1,.e2,prob) - dbinom(x,.e2,prob)) *  (if (log) dbinom(x, size, prob) else 1)  })  # Plot the derivative of dbinom, with respect to prob, calculated by # the corrected version of Deriv(), as a function of probability. Ddb lt;- Deriv(dbinom,"prob") plot(function(p){Ddb(3,8,p)},from=0,to=1,xlab="parameter "prob"",  ylab="",main="derivative of dbinom, corrected") abline(v=3/8,col="red") abline(h=0,col="blue")  

Вы заметите, что производная от dbinom() должно быть положительным для вероятности lt; 3/8 и отрицательным для вероятности gt; 3/8. Моя исправленная версия обладает этим свойством, в то время как производная, полученная из исправленной версии, везде отрицательна.

Может ли кто-нибудь подтвердить, что я прав насчет того, что в пакете Deriv есть ошибка? (т. Е. что я не совершаю какую-то глупую ошибку?)

Проблема (2). Я «перепроверил» вычисления, выполняемые Deriv (), применив эту функцию к «собственной» версии dbinom (), для которой не требуется никаких специальных правил. Я также применил (исправленную версию) функции Deriv() к dbinom(). Код, который я использовал, выглядит следующим образом:

 # # Script demo02. # library(Deriv)   # Replace what I believe to be incorrect code for the rule for # differentiating dbinom() with what I believe to be correct code. # This rule should, strictly speaking, be placed in a new environment # rather than over-writing the existing rule, but this seems to # break down when second derivatives are taken. drule[["dbinom"]] lt;- alist(x=NULL,size=NULL,prob={  .e1 lt;- 1 - prob  .e2 lt;- size - 1  if (x == 0)  -(x * .e1^(x - 1))  else if (x == size)  prob^.e2 * size  else size*(dbinom(x-1,.e2,prob) - dbinom(x,.e2,prob)) *  (if (log) dbinom(x, size, prob) else 1)  },log=NULL)  fooB1 lt;- function(x,prob,size) {  dbinom(x,size,prob) }  fooB2 lt;- function(x,prob,size) {  choose(size,x)*prob^x*(1-prob)^(size-x) }  dfooB1 lt;- Deriv(fooB1,"prob") dfooB2 lt;- Deriv(fooB2,"prob") d2fooB1 lt;- Deriv(fooB1,"prob",nderiv=2) d2fooB2 lt;- Deriv(fooB2,"prob",nderiv=2)  vB1 lt;- fooB1(x=3,prob=0.6,size=8) vB2 lt;- fooB2(x=3,prob=0.6,size=8) dB1 lt;- dfooB1(x=3,prob=0.6,size=8) dB2 lt;- dfooB2(x=3,prob=0.6,size=8) d2B1 lt;- d2fooB1(x=3,prob=0.6,size=8) d2B2 lt;- d2fooB2(x=3,prob=0.6,size=8)  

Если вы запустите этот код, вы увидите, что значения функций (vB1 и vB2) совпадают, оба имеют значение 0,123863. Аналогично, значения первой производной dB1 и dB2 совпадают: -0,9289728.

Однако вторые производные не согласны. Значение d2B1 равно 1,769472, тогда как значение d2B2 равно 2,064384. Я понятия не имею, какой из этих ответов (если таковой имеется) правильный.

Что-то (правило цепочки?) работает не так, как следовало бы.

Есть ли какие-либо действия, которые я могу предпринять, чтобы устранить это несоответствие?