Wednesday, July 6, 2016

bug: missing index of an assignment

A bug of missing index during assignment costed me a whole afternoon to fix. I tried various ways to trace this bug.


for( i in 3:4){
  tb = read.table( paste("../qinlab_rls/",files[i],sep=''), sep="\t")
  GompFlex = flexsurvreg(formula = Surv(tb[,1]) ~ 1, dist = 'gompertz')
  WeibFlex = flexsurvreg(formula = Surv(tb[,1]) ~ 1, dist = 'weibull')

  report$avgLS[i] =  mean(tb[,1])
  report$stdLS[i] =  sd(tb[,1])
  report$CV[i] = report$stdLS[i] / report$avgLS[i]

  report$GompGFlex[i] = GompFlex$res[1,1]
  report$GompRFlex[i] = GompFlex$res[2,1]
  report$GompLogLikFlex[i] = round(GompFlex$loglik, 1)
  report$GompAICFlex[i] = round(GompFlex$AIC)

  report$WeibShapeFlex[i] = WeibFlex$res[1,1]
  report$WeibRateFlex[i] = WeibFlex$res[2,1]
  report$WeibLogLikFlex[i] = round(WeibFlex$loglik, 1)  
  report$WeibAICFlex[i] = round(WeibFlex$AIC)

  Rhat = report$GompRFlex[i]; # 'i' was missing. a bug costed HQ a whole afternoon.
  Ghat = report$GompGFlex[i];
  nhat = 5; 
  t0= (nhat-1)/Ghat;
  fitBinom = optim ( c(Rhat, t0, nhat),  llh.binomialMortality.single.run,
                     lifespan=tb[,1], method="L-BFGS-B", 
                     lower=c(1E-10,0.05, 1), upper=c(10,100,20) );  
  
  report[i, c("R", "t0", "n")] = fitBinom$par[1:3]
}

No comments:

Post a Comment