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.
 A: 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)

