# Estimating polymomial coefficient in R

I used a Taylor series to expand log(1 - ax) so I could estimate the value of parameter 'a'.

The expansion becomes -ax - a^2*x^2/2 - a^3*x^3/3 . . .

Now I need to estimate the parameter 'a' using regression and for simplicity I am only using the first 3 terms in the expansion.

The equation to be estimated becomes y ~ ax + a^2*x^2/2 + a^3*x^3/3 [I have absorbed the negative sign on the left hand side of the equation]

I wanted to ask if there is a way to estimate the coefficients a , a^2 and a^3 in the above equation, keeping in mind that all the three coefficients are powers of each other.

Is there a package in R for this?

Please do note that the Taylor series expansion was necessary as there are several other terms in the original equation which I haven't mentioned here.

Edit:

The original equation I have is:

Y ~ (1 - aX)(B^b)(C^c)(D^d)

In the above equation I have to estimate a,b,c,d, where a is to be estimated as aconstant while b,c and d as smooth splines.

So I have taken log on both side, which makes it:

log(Y) ~ log(1 - aX) + blog(B) + clog(C) + d*log(D)

If there is a better way to approach the entire equation, do mention.

• It's no simplification--instead, it creates errors--to use the Taylor series instead of the function itself. I would like to suggest you will get answers that are much more useful for your problem if you would ask us about the problem you actually have, rather than pursuing this approach. – whuber Aug 21 '17 at 20:08
• Worth a read for the OP: meta.stackexchange.com/questions/66377/what-is-the-xy-problem – Matthew Drury Aug 21 '17 at 20:10
• even if there is a better way to approach the original problem, I would still like to know how to estimate coefficients when the coefficients are the power of each other. – show_stopper Aug 21 '17 at 20:13
• One option would be using non-linear least squares to estimate a,b,c,d without doing any uncessary approximations. This algorithm is available via the nls function in R. – Jarle Tufto Aug 21 '17 at 20:31
• thanks alot for all the help that has been provided. I haven't figured out the best way to solve the problem yet, but there are several directions which I have gotten which I will be pursuing. Firstly, as mentioned I will be looking at the nls function to optimize the function. Secondly, given that as suggested, log-normal distribution assumption might not the the best way, so I will be looking into how to estimate the parameters assuming a Poisson or NB distribution – show_stopper Aug 23 '17 at 17:31

In this particular problem, $Y$ is linear in $a$, holding $b,c,d$ constant, and $\log Y$ is linear in $b,c,d$ holding $a$ constant. So we could probably use coordinate descent and least squares. Although, it's entirely possible that nls would be faster.

(Didn't test code) Given vectors $Y,X,B,C,D$,

    #setup maximum iteration, tolerance, and initialize parameters
max.iter <- 50
tol <- 1e-4
i <- a <- b <- c <- d <- 1
log.y <- log(Y); log.b <- log(B)
log.c <- log(C); log.d <- log(D)
while(i <= max.iter){
#a-update holdig b,c,d constant
y.new <- Y/(B^b*C^c*D^d) - 1
my.qr <- qr(-X)
a <- qr.coef(my.qr, y.new)
#b,c,d update holding a constant
y.new <- log.y - log(1-a*X)
my.qr <- qr(cbind(log.b, log.c, log.d))
bcd <- qr.coef(my.qr, y.new)
b <- bcd[1]
c <- bcd[2]
d <- bcd[3]
#check convergence
param.new <- c(a, bcd)
if(norm(param.new-param.old) < tol) break
else param.old <- param.new
i <- i + 1
}
print(param.new)

• Given additional information that the error is log-normal, this solves the wrong optimization problem. – Jirapat Samranvedhya Aug 22 '17 at 13:23