On Mon, May 5, 2014 at 10:16 PM, David Winsemius <dwinsem...@comcast.net> wrote:
>
> On May 5, 2014, at 6:05 PM, Gabor Grothendieck wrote:
>
>> On Mon, May 5, 2014 at 9:43 AM, Niloofar.Javanrouh
>> <javanrou...@yahoo.com> wrote:
>>>
>>>
>>> hello,
>>> i want to differentiate of L with respect to b
>>> when:
>>>
>>> L= k*ln (k/(k+mu)) + sum(y) * ln (1-(k/mu+k)) #(negative binomial ln
>>> likelihood)
>>> and
>>> ln(mu/(mu+k)) = a+bx #link function
>>>
>>> how can i do it in R?
>>
>> Try this. First we solve for 'mu' in terms of the other variables
>> using the link equation:
>>
>>> library(Ryacas)
>>> k <- Sym("k")
>>> mu <- Sym("mu")
>>> y <- Sym("y")
>>> L <- Sym("L")
>>> a <- Sym("a")
>>> b <- Sym("b")
>>> x <- Sym("x")
>>> sumy <- Sym("sumy")
>>> Solve(log(mu/(mu+k)) == a+b*x, "mu")
>> expression(list(mu == k * exp(a + b * x)/(1 - exp(a + b * x))))
>>>
>>
>> Now in 'k*log(k/(k+mu)) + sumy * log(1-(k/mu+k))' substitute 'k *
>> exp(a + b * x)/(1 - exp(a + b * x)))' for 'mu' using 'Subst' and take
>> the derivative with respect to 'a' using 'deriv':
>>
>>> s <- Subst(k*log(k/(k+mu)) + sumy * log(1-(k/mu+k)), mu,
>> + k * exp(a + b * x)/(1 - exp(a + b * x)))
>>> deriv(s, a)
>> expression(sumy * (k * ((1 - exp(a + b * x)) * (k * exp(a + b *
>> x)) + k * exp(a + b * x)^2))/((1 - exp(a + b * x))^2 * (k *
>> exp(a + b * x)/(1 - exp(a + b * x)))^2 * (1 - (k * (1 - exp(a +
>> b * x))/(k * exp(a + b * x)) + k))) - k * ((k + k * exp(a +
>> b * x)/(1 - exp(a + b * x))) * (k * ((1 - exp(a + b * x)) *
>> (k * exp(a + b * x)) + k * exp(a + b * x)^2)))/((1 - exp(a +
>> b * x))^2 * ((k + k * exp(a + b * x)/(1 - exp(a + b * x)))^2 *
>> k)))
>>
>
> But can we maybe get the Taylor series approximation to first or second order?
>
> Taylor(d, a, 0, 2)
expression(sumy * (k * ((1 - exp(b * x)) * (k * exp(b * x)) +
k * exp(b * x)^2))/((1 - exp(b * x))^2 * (k * exp(b * x)/(1 -
exp(b * x)))^2 * (1 - (k * (1 - exp(b * x))/(k * exp(b *
x)) + k))) - k * ((k + k * exp(b * x)/(1 - exp(b * x))) *
(k * ((1 - exp(b * x)) * (k * exp(b * x)) + k * exp(b * x)^2)))/((1 -
exp(b * x))^2 * ((k + k * exp(b * x)/(1 - exp(b * x)))^2 *
k)) + (((1 - exp(b * x))^2 * (k * exp(b * x)/(1 - exp(b *
x)))^2 * (1 - (k * (1 - exp(b * x))/(k * exp(b * x)) + k)) *
(sumy * (k * ((1 - exp(b * x)) * (k * exp(b * x)) - k * exp(b *
x)^2 + k * (2 * exp(b * x)^2)))) - sumy * (k * ((1 -
exp(b * x)) * (k * exp(b * x)) + k * exp(b * x)^2)) * ((1 -
exp(b * x))^2 * (k * exp(b * x)/(1 - exp(b * x)))^2 * (k *
((1 - exp(b * x)) * (k * exp(b * x)) + k * exp(b * x)^2))/((1 -
exp(b * x))^2 * (k * exp(b * x)/(1 - exp(b * x)))^2) + ((1 -
exp(b * x))^2 * (k * exp(b * x) * (2 * ((1 - exp(b * x)) *
(k * exp(b * x)) + k * exp(b * x)^2)))/(1 - exp(b * x))^3 +
-2 * exp(b * x) * (1 - exp(b * x)) * (k * exp(b * x)/(1 -
exp(b * x)))^2) * (1 - (k * (1 - exp(b * x))/(k * exp(b *
x)) + k))))/((1 - exp(b * x))^2 * (k * exp(b * x)/(1 - exp(b *
x)))^2 * (1 - (k * (1 - exp(b * x))/(k * exp(b * x)) + k)))^2 -
((1 - exp(b * x))^2 * ((k + k * exp(b * x)/(1 - exp(b * x)))^2 *
k) * (k * ((k + k * exp(b * x)/(1 - exp(b * x))) * (k *
((1 - exp(b * x)) * (k * exp(b * x)) - k * exp(b * x)^2 +
k * (2 * exp(b * x)^2))) + k * ((1 - exp(b * x)) *
(k * exp(b * x)) + k * exp(b * x)^2)^2/(1 - exp(b * x))^2)) -
k * ((k + k * exp(b * x)/(1 - exp(b * x))) * (k * ((1 -
exp(b * x)) * (k * exp(b * x)) + k * exp(b * x)^2))) *
((1 - exp(b * x))^2 * (k * ((k + k * exp(b * x)/(1 -
exp(b * x))) * (2 * ((1 - exp(b * x)) * (k *
exp(b * x)) + k * exp(b * x)^2))))/(1 - exp(b *
x))^2 + -2 * exp(b * x) * (1 - exp(b * x)) *
((k + k * exp(b * x)/(1 - exp(b * x)))^2 * k)))/((1 -
exp(b * x))^2 * ((k + k * exp(b * x)/(1 - exp(b * x)))^2 *
k))^2) * a + a^2 * ((((1 - exp(b * x))^2 * (k * exp(b *
x)/(1 - exp(b * x)))^2 * (1 - (k * (1 - exp(b * x))/(k *
exp(b * x)) + k)))^2 * ((1 - exp(b * x))^2 * (k * exp(b *
x)/(1 - exp(b * x)))^2 * (1 - (k * (1 - exp(b * x))/(k *
exp(b * x)) + k)) * (sumy * (k * ((1 - exp(b * x)) * (k *
exp(b * x)) - k * exp(b * x)^2 - k * (2 * exp(b * x)^2) +
k * (4 * exp(b * x)^2)))) + ((1 - exp(b * x))^2 * (k * exp(b *
x)/(1 - exp(b * x)))^2 * (k * ((1 - exp(b * x)) * (k * exp(b *
x)) + k * exp(b * x)^2))/((1 - exp(b * x))^2 * (k * exp(b *
x)/(1 - exp(b * x)))^2) + ((1 - exp(b * x))^2 * (k * exp(b *
x) * (2 * ((1 - exp(b * x)) * (k * exp(b * x)) + k * exp(b *
x)^2)))/(1 - exp(b * x))^3 + -2 * exp(b * x) * (1 - exp(b *
x)) * (k * exp(b * x)/(1 - exp(b * x)))^2) * (1 - (k * (1 -
exp(b * x))/(k * exp(b * x)) + k))) * (sumy * (k * ((1 -
exp(b * x)) * (k * exp(b * x)) - k * exp(b * x)^2 + k * (2 *
exp(b * x)^2)))) - (sumy * (k * ((1 - exp(b * x)) * (k *
exp(b * x)) + k * exp(b * x)^2)) * (((1 - exp(b * x))^2 *
(k * exp(b * x)/(1 - exp(b * x)))^2 * ((1 - exp(b * x))^2 *
(k * exp(b * x)/(1 - exp(b * x)))^2 * (k * ((1 - exp(b *
x)) * (k * exp(b * x)) - k * exp(b * x)^2 + k * (2 * exp(b *
x)^2))) + ((1 - exp(b * x))^2 * (k * exp(b * x) * (2 * ((1 -
exp(b * x)) * (k * exp(b * x)) + k * exp(b * x)^2)))/(1 -
exp(b * x))^3 + -2 * exp(b * x) * (1 - exp(b * x)) * (k *
exp(b * x)/(1 - exp(b * x)))^2) * (k * ((1 - exp(b * x)) *
(k * exp(b * x)) + k * exp(b * x)^2))) - (1 - exp(b * x))^2 *
(k * exp(b * x)/(1 - exp(b * x)))^2 * (k * ((1 - exp(b *
x)) * (k * exp(b * x)) + k * exp(b * x)^2)) * ((1 - exp(b *
x))^2 * (k * exp(b * x) * (2 * ((1 - exp(b * x)) * (k * exp(b *
x)) + k * exp(b * x)^2)))/(1 - exp(b * x))^3 + -2 * exp(b *
x) * (1 - exp(b * x)) * (k * exp(b * x)/(1 - exp(b * x)))^2))/((1 -
exp(b * x))^2 * (k * exp(b * x)/(1 - exp(b * x)))^2)^2 +
(((1 - exp(b * x))^2 * (k * exp(b * x) * (2 * ((1 - exp(b *
x)) * (k * exp(b * x)) + k * exp(b * x)^2)))/(1 - exp(b *
x))^3 + -2 * exp(b * x) * (1 - exp(b * x)) * (k * exp(b *
x)/(1 - exp(b * x)))^2) * (k * ((1 - exp(b * x)) * (k *
exp(b * x)) + k * exp(b * x)^2))/((1 - exp(b * x))^2 *
(k * exp(b * x)/(1 - exp(b * x)))^2) + (((1 - exp(b *
x))^3 * ((1 - exp(b * x))^2 * (k * exp(b * x) * (2 *
((1 - exp(b * x)) * (k * exp(b * x)) - k * exp(b * x)^2 +
k * (2 * exp(b * x)^2))) + k * exp(b * x) * (2 *
((1 - exp(b * x)) * (k * exp(b * x)) + k * exp(b * x)^2))) +
-2 * exp(b * x) * (1 - exp(b * x)) * (k * exp(b * x) *
(2 * ((1 - exp(b * x)) * (k * exp(b * x)) + k * exp(b *
x)^2)))) - (1 - exp(b * x))^2 * (k * exp(b *
x) * (2 * ((1 - exp(b * x)) * (k * exp(b * x)) + k *
exp(b * x)^2))) * (-3 * exp(b * x) * (1 - exp(b * x))^2))/(1 -
exp(b * x))^6 + (-2 * exp(b * x) * (1 - exp(b * x)) *
(k * exp(b * x) * (2 * ((1 - exp(b * x)) * (k * exp(b *
x)) + k * exp(b * x)^2)))/(1 - exp(b * x))^3 + (-2 *
exp(b * x) * (1 - exp(b * x)) - -2 * exp(b * x)^2) *
(k * exp(b * x)/(1 - exp(b * x)))^2)) * (1 - (k * (1 -
exp(b * x))/(k * exp(b * x)) + k)))) + sumy * (k * ((1 -
exp(b * x)) * (k * exp(b * x)) - k * exp(b * x)^2 + k * (2 *
exp(b * x)^2))) * ((1 - exp(b * x))^2 * (k * exp(b * x)/(1 -
exp(b * x)))^2 * (k * ((1 - exp(b * x)) * (k * exp(b * x)) +
k * exp(b * x)^2))/((1 - exp(b * x))^2 * (k * exp(b * x)/(1 -
exp(b * x)))^2) + ((1 - exp(b * x))^2 * (k * exp(b * x) *
(2 * ((1 - exp(b * x)) * (k * exp(b * x)) + k * exp(b * x)^2)))/(1 -
exp(b * x))^3 + -2 * exp(b * x) * (1 - exp(b * x)) * (k *
exp(b * x)/(1 - exp(b * x)))^2) * (1 - (k * (1 - exp(b *
x))/(k * exp(b * x)) + k))))) - ((1 - exp(b * x))^2 * (k *
exp(b * x)/(1 - exp(b * x)))^2 * (1 - (k * (1 - exp(b * x))/(k *
exp(b * x)) + k)) * (sumy * (k * ((1 - exp(b * x)) * (k *
exp(b * x)) - k * exp(b * x)^2 + k * (2 * exp(b * x)^2)))) -
sumy * (k * ((1 - exp(b * x)) * (k * exp(b * x)) + k * exp(b *
x)^2)) * ((1 - exp(b * x))^2 * (k * exp(b * x)/(1 - exp(b *
x)))^2 * (k * ((1 - exp(b * x)) * (k * exp(b * x)) +
k * exp(b * x)^2))/((1 - exp(b * x))^2 * (k * exp(b *
x)/(1 - exp(b * x)))^2) + ((1 - exp(b * x))^2 * (k *
exp(b * x) * (2 * ((1 - exp(b * x)) * (k * exp(b * x)) +
k * exp(b * x)^2)))/(1 - exp(b * x))^3 + -2 * exp(b *
x) * (1 - exp(b * x)) * (k * exp(b * x)/(1 - exp(b *
x)))^2) * (1 - (k * (1 - exp(b * x))/(k * exp(b * x)) +
k)))) * (2 * ((1 - exp(b * x))^2 * (k * exp(b * x)/(1 -
exp(b * x)))^2 * (k * ((1 - exp(b * x)) * (k * exp(b * x)) +
k * exp(b * x)^2))/((1 - exp(b * x))^2 * (k * exp(b * x)/(1 -
exp(b * x)))^2) + ((1 - exp(b * x))^2 * (k * exp(b * x) *
(2 * ((1 - exp(b * x)) * (k * exp(b * x)) + k * exp(b * x)^2)))/(1 -
exp(b * x))^3 + -2 * exp(b * x) * (1 - exp(b * x)) * (k *
exp(b * x)/(1 - exp(b * x)))^2) * (1 - (k * (1 - exp(b *
x))/(k * exp(b * x)) + k))) * ((1 - exp(b * x))^2 * (k *
exp(b * x)/(1 - exp(b * x)))^2 * (1 - (k * (1 - exp(b * x))/(k *
exp(b * x)) + k)))))/((1 - exp(b * x))^2 * (k * exp(b * x)/(1 -
exp(b * x)))^2 * (1 - (k * (1 - exp(b * x))/(k * exp(b *
x)) + k)))^4 - (((1 - exp(b * x))^2 * ((k + k * exp(b * x)/(1 -
exp(b * x)))^2 * k))^2 * ((1 - exp(b * x))^2 * ((k + k *
exp(b * x)/(1 - exp(b * x)))^2 * k) * (k * ((k + k * exp(b *
x)/(1 - exp(b * x))) * (k * ((1 - exp(b * x)) * (k * exp(b *
x)) - k * exp(b * x)^2 - k * (2 * exp(b * x)^2) + k * (4 *
exp(b * x)^2))) + k * ((1 - exp(b * x)) * (k * exp(b * x)) -
k * exp(b * x)^2 + k * (2 * exp(b * x)^2)) * ((1 - exp(b *
x)) * (k * exp(b * x)) + k * exp(b * x)^2)/(1 - exp(b * x))^2 +
((1 - exp(b * x))^2 * (k * (2 * ((1 - exp(b * x)) * (k *
exp(b * x)) - k * exp(b * x)^2 + k * (2 * exp(b * x)^2)) *
((1 - exp(b * x)) * (k * exp(b * x)) + k * exp(b * x)^2))) -
k * ((1 - exp(b * x)) * (k * exp(b * x)) + k * exp(b *
x)^2)^2 * (-2 * exp(b * x) * (1 - exp(b * x))))/(1 -
exp(b * x))^4)) + ((1 - exp(b * x))^2 * (k * ((k + k *
exp(b * x)/(1 - exp(b * x))) * (2 * ((1 - exp(b * x)) * (k *
exp(b * x)) + k * exp(b * x)^2))))/(1 - exp(b * x))^2 + -2 *
exp(b * x) * (1 - exp(b * x)) * ((k + k * exp(b * x)/(1 -
exp(b * x)))^2 * k)) * (k * ((k + k * exp(b * x)/(1 - exp(b *
x))) * (k * ((1 - exp(b * x)) * (k * exp(b * x)) - k * exp(b *
x)^2 + k * (2 * exp(b * x)^2))) + k * ((1 - exp(b * x)) *
(k * exp(b * x)) + k * exp(b * x)^2)^2/(1 - exp(b * x))^2)) -
(k * ((k + k * exp(b * x)/(1 - exp(b * x))) * (k * ((1 -
exp(b * x)) * (k * exp(b * x)) + k * exp(b * x)^2))) *
(((1 - exp(b * x))^2 * ((1 - exp(b * x))^2 * (k * ((k +
k * exp(b * x)/(1 - exp(b * x))) * (2 * ((1 - exp(b *
x)) * (k * exp(b * x)) - k * exp(b * x)^2 + k * (2 *
exp(b * x)^2))) + 2 * ((1 - exp(b * x)) * (k * exp(b *
x)) + k * exp(b * x)^2)^2/(1 - exp(b * x))^2)) +
-2 * exp(b * x) * (1 - exp(b * x)) * (k * ((k + k *
exp(b * x)/(1 - exp(b * x))) * (2 * ((1 - exp(b *
x)) * (k * exp(b * x)) + k * exp(b * x)^2))))) -
(1 - exp(b * x))^2 * (k * ((k + k * exp(b * x)/(1 -
exp(b * x))) * (2 * ((1 - exp(b * x)) * (k *
exp(b * x)) + k * exp(b * x)^2)))) * (-2 * exp(b *
x) * (1 - exp(b * x))))/(1 - exp(b * x))^4 +
(-2 * exp(b * x) * (1 - exp(b * x)) * (k * ((k +
k * exp(b * x)/(1 - exp(b * x))) * (2 * ((1 -
exp(b * x)) * (k * exp(b * x)) + k * exp(b *
x)^2))))/(1 - exp(b * x))^2 + (-2 * exp(b * x) *
(1 - exp(b * x)) - -2 * exp(b * x)^2) * ((k +
k * exp(b * x)/(1 - exp(b * x)))^2 * k))) + k *
((k + k * exp(b * x)/(1 - exp(b * x))) * (k * ((1 - exp(b *
x)) * (k * exp(b * x)) - k * exp(b * x)^2 + k * (2 *
exp(b * x)^2))) + k * ((1 - exp(b * x)) * (k * exp(b *
x)) + k * exp(b * x)^2)^2/(1 - exp(b * x))^2) * ((1 -
exp(b * x))^2 * (k * ((k + k * exp(b * x)/(1 - exp(b *
x))) * (2 * ((1 - exp(b * x)) * (k * exp(b * x)) + k *
exp(b * x)^2))))/(1 - exp(b * x))^2 + -2 * exp(b * x) *
(1 - exp(b * x)) * ((k + k * exp(b * x)/(1 - exp(b *
x)))^2 * k)))) - ((1 - exp(b * x))^2 * ((k + k * exp(b *
x)/(1 - exp(b * x)))^2 * k) * (k * ((k + k * exp(b * x)/(1 -
exp(b * x))) * (k * ((1 - exp(b * x)) * (k * exp(b * x)) -
k * exp(b * x)^2 + k * (2 * exp(b * x)^2))) + k * ((1 - exp(b *
x)) * (k * exp(b * x)) + k * exp(b * x)^2)^2/(1 - exp(b *
x))^2)) - k * ((k + k * exp(b * x)/(1 - exp(b * x))) * (k *
((1 - exp(b * x)) * (k * exp(b * x)) + k * exp(b * x)^2))) *
((1 - exp(b * x))^2 * (k * ((k + k * exp(b * x)/(1 - exp(b *
x))) * (2 * ((1 - exp(b * x)) * (k * exp(b * x)) + k *
exp(b * x)^2))))/(1 - exp(b * x))^2 + -2 * exp(b * x) *
(1 - exp(b * x)) * ((k + k * exp(b * x)/(1 - exp(b *
x)))^2 * k))) * (2 * ((1 - exp(b * x))^2 * (k * ((k +
k * exp(b * x)/(1 - exp(b * x))) * (2 * ((1 - exp(b * x)) *
(k * exp(b * x)) + k * exp(b * x)^2))))/(1 - exp(b * x))^2 +
-2 * exp(b * x) * (1 - exp(b * x)) * ((k + k * exp(b * x)/(1 -
exp(b * x)))^2 * k)) * ((1 - exp(b * x))^2 * ((k + k *
exp(b * x)/(1 - exp(b * x)))^2 * k))))/((1 - exp(b * x))^2 *
((k + k * exp(b * x)/(1 - exp(b * x)))^2 * k))^4)/2)
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
