biasbox <- function(nls.obj){
  theta <- summary(nls.obj)$coef[,1]
  sd.theta <- summary(nls.obj)$coef[,2]
  F <- attr(nls.obj$m$fitted(), "gradient")
  H <- attr(nls.obj$m$fitted(), "hessian")
  sig <- summary(nls.obj)$sigma
  n <- dim(F)[1]
  FlFi <- t(F)%*%F
  d <- -(sig^2/2)*sapply(1:n, function(x){
    sum(diag(solve(FlFi)%*%H[x, , ]))})
  bias <- as.vector(solve(FlFi)%*%t(F)%*%d)
  names(bias) <- names(coef(nls.obj))
  bias.sd <- 100*bias/sd.theta
  bias.th <- 100*bias/theta
  return(list("viés bruto"=bias,
              "%viés(theta)"=bias.th,
              "%viés(sd(theta))"=bias.sd))
}