f(x)=(a*b*alpha*beta*theta*x^(-beta-1)*exp(-alpha*x^(-beta))*exp(theta*exp(-alpha*x^(-beta)))*((exp(theta)-1)^(-1))) *(((1-exp(theta*exp(-alpha*x^(-beta))))/(1-exp(theta)))^(a-1))*((1 - ((1-exp(theta*exp(-alpha*x^(-beta))))/(1-exp(theta)))^a)^(b-1))



 # xvec my data
xvec <- c(15,20,38 ,42, 61 ,76, 86, 98, 121, 146,
          149, 157, 175, 176, 180, 180, 198, 220, 224, 251,
          264, 282, 321, 325, 653)
n=length(xvec)


fn <- function(theta) {
  a=theta[1]
  b=theta[2]
  c=theta[3]
  d=theta[4]
  e=theta[5]
  -(n*log(a)+n*log(b)+n*log(c)+n*log(d)+n*log(e)-n*a*b*log(exp(e)-1)-(d+1)*sum(log(xvec))-c*sum(1/xvec^(d))+e*sum(exp(-c*(1/xvec^(d)))))
  +(a-1)*sum(e*exp(-c*(1/xvec^(d))))+(b-1)*sum(log((exp(e)-1)^(a)-exp(e*exp(-c*(1/xvec^(d))))))

}

m=nlm(fn, theta <- c(10,20,3,4,6), hessian=TRUE)

$minimum
[1] -4471.137

$estimate
[1]  6.9621176 18.4829665  3.0000015  3.9999876  0.7483349

$gradient
[1]  3.582204e+06 -2.621225e+02  8.706031e+00 -7.659562e+01  2.914793e+08

$hessian
              [,1]          [,2]          [,3]          [,4]          [,5]
[1,] -1.166002e+09  6.840246e+04 -8.306624e+03  74236.917635 -1.158880e+10
[2,]  6.840246e+04 -2.662299e-07  4.979732e-01     -4.380081  1.134472e+06
[3,] -8.306624e+03  4.979732e-01 -1.796654e+00    -10.434046 -8.428344e+04
[4,]  7.423692e+04 -4.380081e+00 -1.043405e+01    100.060620  7.423287e+05
[5,] -1.158880e+10  1.134472e+06 -8.428344e+04 742328.727945 -1.682691e+11

$code
[1] 2

$iterations
[1] 29