statistics/rt_data.r

   1
   2
   3 source("functions.r");
   4
   5 # system("parse_rt_data.py 3 > rt_data.csv");
   6 #t <- read.csv('rt_data.csv', sep=',', header=TRUE)
   7 t <- read.csv('rt_data_2004-2010.csv', sep=',', header=TRUE)
   8
   9 par(mfrow=c(2,1))
  10
  11 h<-hist(log(log(t$replies)), breaks=50)
  12 lines(h$breaks[which(h$counts!=0)], h$counts[which(h$counts!=0)])
  13 h<-hist(log(log(log(t$replies))), breaks=50)
  14 lines(h$breaks[which(h$counts!=0)], h$counts[which(h$counts!=0)])
  15
  16
  17 par(mfrow=c(1,1))
  18
  19 t2 <- t[which(t$complete == 1),]
  20 d <- (t2$lastreply - t2$start)/(60*60)
  21
  22 #start_image("rt_hist_ttc_1000.png")
  23 #hist(d[which(d<1000)], xlab="hours from creation to last reply", breaks=30)
  24 #end_image()
  25 #
  26 #start_image("rt_hist_ttc_200.png")
  27 #hist(d[which(d<200)], xlab="hours from creation to last reply", breaks=30)
  28 #end_image()
  29 #
  30 #start_image("rt_hist_ttc_50.png")
  31 #hist(d[which(d<50)], xlab="hours from creation to last reply", breaks=30)
  32 #end_image()
  33 #
  34 #start_image("rt_hist_ttc_10.png")
  35 #hist(d[which(d<10)], xlab="hours from creation to last reply", breaks=30)
  36 #end_image()
  37 #
  38 #d2 <- (t2$lastreply - t2$start)
  39 #h<-hist(log(d2), plot=F, breaks=50)
  40 #lines(h$breaks[which(h$counts!=0)], h$counts[which(h$counts!=0)])
  41
  42 # this doesn't work as I would like.  I think the bins aren't as I expect
  43 #h <- hist(d, plot=F, breaks=c(seq(0,max(d)+1, .1)))
  44 #plot(h$counts, log="x", pch=20, col="blue",
  45 #       main="Log-normal distribution",
  46 #       xlab="Value", ylab="Frequency")
  47
  48 #plot(log(d2))
  49 #plot(ecdf(d2))
  50
  51 tstamp_45 <-unclass(as.POSIXct("2005-01-01", origin="1970-01-01"))[1]
  52 tstamp_56 <-unclass(as.POSIXct("2006-01-01", origin="1970-01-01"))[1]
  53 tstamp_67 <-unclass(as.POSIXct("2007-01-01", origin="1970-01-01"))[1]
  54 tstamp_78 <-unclass(as.POSIXct("2008-01-01", origin="1970-01-01"))[1]
  55 tstamp_89 <-unclass(as.POSIXct("2009-01-01", origin="1970-01-01"))[1]
  56 tstamp_90 <-unclass(as.POSIXct("2010-01-01", origin="1970-01-01"))[1]
  57
  58
  59 t_4 <- t2[which( t2$start <  tstamp_45 ),]
  60 t_5 <- t2[which( t2$start >= tstamp_45 & t2$start < tstamp_56 ),]
  61 t_6 <- t2[which( t2$start >= tstamp_56 & t2$start < tstamp_67 ),]
  62 t_7 <- t2[which( t2$start >= tstamp_67 & t2$start < tstamp_78 ),]
  63 t_8 <- t2[which( t2$start >= tstamp_78 & t2$start < tstamp_89 ),]
  64 t_9 <- t2[which( t2$start >= tstamp_89 & t2$start < tstamp_90 ),]
  65 t_10 <- t2[which( t2$start >= tstamp_90 ),]
  66
  67 par(mfrow=c(4,1))
  68 plot_rt_hist(t_4)
  69 plot_rt_hist(t_5)
  70 plot_rt_hist(t_6)
  71 plot_rt_hist(t_7)
  72 plot_rt_hist(t_8)
  73 plot_rt_hist(t_9)
  74 par(mfrow=c(1,1))
  75
  76 start_image("rt_support_seasonal.png")
  77 par(mfrow=c(6,1))
  78 par(mai=c(.3,.3,.3,.3))
  79
  80 # start dates on Sunday to align all weeks with weekend boundaries.
  81 year_hist(t_4, "2004", "2003/12/28", "2005/1/7", 85)
  82 year_hist(t_5, "2005", "2005/1/2", "2006/1/7", 85)
  83 year_hist(t_6, "2006", "2006/1/1", "2007/1/7", 85)
  84 year_hist(t_7, "2007", "2006/12/31", "2008/1/7", 85)
  85 year_hist(t_8, "2008", "2007/12/30", "2009/1/7", 85)
  86 year_hist(t_9, "2009", "2008/12/28", "2010/1/30", 85)
  87 end_image()
  88
  89 h4<-year_hist(t_4, "2004", "2003/12/28", "2005/2/7", 0, type='month', fmt="%b")
  90 h5<-year_hist(t_5, "2005", "2005/1/2", "2006/2/7", 0, type='month', fmt="%b")
  91 h6<-year_hist(t_6, "2006", "2006/1/1", "2007/2/7", 0, type='month', fmt="%b")
  92 h7<-year_hist(t_7, "2007", "2006/12/31", "2008/2/7", 0, type='month', fmt="%b")
  93 h8<-year_hist(t_8, "2008", "2007/12/30", "2009/2/7", 0, type='month', fmt="%b")
  94 h9<-year_hist(t_9, "2009", "2008/12/28", "2010/1/30", 0, type='month', fmt="%b")
  95
  96 hall<-year_hist(t2, "200x", "2004/1/1", "2010/3/28", 0, type='month', fmt="%b")
  97
  98 threshold <- function (hall, d, from, to, type, fmt="%b")
  99 {
 100     dates <-seq(as.Date(from), as.Date(to), type)
 101     months <- format(dates, fmt)
 102     hbreaks<-unclass(as.POSIXct(dates))
 103
 104     x<-seq(1,length(hall$breaks))
 105     a_x<-x[which(hall$counts>d)]
 106     a_y<-hall$counts[which(hall$counts>d)]
 107     b_x<-x[which(hall$counts<d)]
 108     b_y<-hall$counts[which(hall$counts<d)]
 109
 110     plot(a_x, a_y, type='p', col='red', ylim=c(0,260), xlim=c(0,81), axes=F)
 111     points(b_x, b_y, type='p', col='blue', ylim=c(0,260), xlim=c(0,81))
 112     axis(1, labels=months, at=x)
 113     axis(2)
 114     abline(v=seq(13,length(months),12))
 115 }
 116
 117 years <- 7
 118 b<- seq(1,years*12,12)
 119 yy<-NULL
 120 for (i in seq(1,years) )
 121 {
 122     if ( i+1 > length(b) ) {
 123         yy<- rbind(yy,hall$counts[b[i]:length(hall$counts)])
 124     } else {
 125         yy<- rbind(yy,hall$counts[b[i]:b[i+1]-1])
 126     }
 127 }
 128 yy[7,3:12]<-0   # no data for beyond feb.
 129 y2<-NULL ; for ( i in seq(1,12) ) { y2<-c(y2,sum(yy[,i])) }
 130
 131 start_image('rt_aggregate_months.png', width=600, height=300)
 132 barplot(y2, space=.1, width=.9, col=c('blue','red', 'red', 'red', 'red',
 133     'blue', 'blue', 'red', 'red', 'red', 'blue', 'blue'),
 134     xlab="Months", ylab="Sum of Tickets over 6 years")
 135 axis(1, labels=c('Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug',
 136     'Sep', 'Oct', 'Nov', 'Dec'), at=c(0,1,2,3,4,5,6,7,8,9,10,11)+.5)
 137 end_image()
 138
 139 cc<-NULL ;
 140 for (i in 1:length(yy))
 141 {
 142     if ( t(yy)[i] < 80 )
 143     {
 144         cc<- c(cc, 'blue')
 145     } else {
 146         cc<- c(cc, 'red')
 147     }
 148 }
 149 barplot(yy, col=cc)
 150
 151 # skip 2007
 152 start_image('rt_aggregate_months_no2007.png', width=600, height=300)
 153 y3<-NULL ; for ( i in seq(1,12) ) { y3<-c(y3,sum(yy[1:3,i], yy[5:7,i])) }
 154 barplot(y3, , space=.1, width=.9, col=c('blue','blue', 'red', 'red', 'red',
 155     'blue', 'blue', 'red', 'red', 'red', 'blue', 'blue'),
 156     xlab="Months", ylab="Sum of Tickets over 6 years")
 157 axis(1, labels=c('Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug',
 158     'Sep', 'Oct', 'Nov', 'Dec'), at=c(0,1,2,3,4,5,6,7,8,9,10,11)+.5)
 159 end_image()
 160
 161
 162
 163
 164 par(mai=c(0.7,0.7,0.7,0.7))
 165 par(mfrow=c(1,1))
 166
 167
 168 tstamp <-unclass(as.POSIXct("2008-01-01", origin="1960-01-01"))
 169 t_67 <- t2[which( t2$start <  tstamp[1] ),]
 170 t_89 <- t2[which( t2$start >= tstamp[1] ),]
 171
 172
 173 # install.packages('sn')
 174 require(sn)
 175 par(mfrow=c(6,1))
 176 par(mai=c(0.3,0.3,0.3,0.3))
 177
 178 #start_image("rt_hist_ttc_1000.png")
 179 time_hist <- function (t, lessthan, year, log=T, breaks=30, xlim=c(-4,10), ylim=c(0,150))
 180 {
 181     d <- (t$lastreply - t$start)/(60*60)
 182     main = sprintf("Histogram of d<%s for %s", lessthan, year);
 183     if ( log )
 184     {
 185         d <- log(d[which(d<lessthan)])
 186         avg <- round(mean(d), 2)
 187         main = sprintf("Histogram of d<%s for %s : %s", lessthan, year, avg);
 188         #h<-sn.mle(d, xlab="Hours to Complete Ticket", breaks=30, main=main, xlim=xlim, ylim=ylim)
 189         h<-sn.mle(y=d)
 190     } else {
 191         h<-hist(d[which(d<lessthan)], xlab="Hours to Complete Ticket", breaks=30, main=main, xlim=xlim, ylim=ylim)
 192     }
 193     return (h);
 194 }
 195
 196 median_time_to_resolve_window <- function (t, from, to, window, fmt="%b")
 197 {
 198     # find 'type' range of days
 199     dates <-seq(as.Date(from), as.Date(to), 'week')
 200     months <- format(dates, fmt)
 201     hbreaks<-unclass(as.POSIXct(dates))
 202
 203     xx<-NULL;
 204     yy<-NULL;
 205     yy_sd_high<-NULL;
 206     yy_sd_low<-NULL;
 207     date_index <- NULL;
 208     q_list <- NULL;
 209
 210     x<-seq(-20,20,0.01)
 211     for ( i in seq(1,length(hbreaks)-window-1) )
 212     {
 213         print (sprintf("round %s of %s", i, length(hbreaks)-window-1))
 214         # get range from t
 215         t_sub <- t[which(t$start > hbreaks[i] & t$start<= hbreaks[i+window]),]
 216         if ( length(t_sub$start) <= 1 )  { next }
 217         # take log, then sn.mle -> h
 218         d <- (t_sub$lastreply - t_sub$start)/(60*60)    # hours
 219         d <- log(d)                                     # log(hours)
 220             # sn.mle
 221         print (sprintf("length: %s", length(d)))
 222         q<-quantile(d)
 223         print(q)
 224
 225         date_index <- c(date_index, round(i+window/2))
 226
 227         xx<- c(xx, hbreaks[round(i+window/2)])
 228         q_list <- rbind(q_list, q)
 229
 230     }
 231     m<- months[date_index]
 232     return (cbind(xx,q_list, m))
 233 }
 234 mean_time_to_resolve_window <- function (t, from, to, window, fmt="%b")
 235 {
 236     # find 'type' range of days
 237     dates <-seq(as.Date(from), as.Date(to), 'week')
 238     months <- format(dates, fmt)
 239     hbreaks<-unclass(as.POSIXct(dates))
 240
 241     xx<-NULL;
 242     yy<-NULL;
 243     yy_sd_high<-NULL;
 244     yy_sd_low<-NULL;
 245     date_list <- NULL;
 246
 247     x<-seq(-20,20,0.01)
 248     for ( i in seq(1,length(hbreaks)-window-1) )
 249     {
 250         print (sprintf("round %s of %s", i, length(hbreaks)-window-1))
 251         # get range from t
 252         t_sub <- t[which(t$start > hbreaks[i] & t$start<= hbreaks[i+window]),]
 253         if ( length(t_sub$start) <= 1 )  { next }
 254         # take log, then sn.mle -> h
 255         d <- (t_sub$lastreply - t_sub$start)/(60*60)    # hours
 256         d <- log(d)                                     # log(hours)
 257             # sn.mle
 258         print (sprintf("length: %s", length(d)))
 259         avg<-mean(d)
 260         s<-sd(d)
 261         r<-shapiro.test(d) #, mean=avg, sd=s)
 262         if ( r$statistic < 0.9 ){
 263             print (r);
 264         }
 265
 266         m<-dnorm(x, mean=avg, sd=s)
 267         print(avg)
 268         # find max of y
 269         y_peak <- x[which(m==max(m))]
 270         print(y_peak)
 271         # plot point date, max(y)
 272         xx<- c(xx, hbreaks[round(i+window/2)])
 273         yy<- c(yy, y_peak)
 274         yy_sd_high<- y_peak + s
 275         yy_sd_low <- y_peak - s
 276         date_list <- c(date_list, dates[i])
 277         # plot whisker2(x0,y0,y0_hi,y0_lo)
 278     }
 279     l<-length(months)-window-1
 280     m<- months[1:l]
 281     return (rbind(xx,yy,yy_sd_high, yy_sd_low, m))
 282 }
 283 require(sn)
 284 sknorm_time_to_resolve_window <- function (t, from, to, window, fmt="%b")
 285 {
 286     # find 'type' range of days
 287     dates <-seq(as.Date(from), as.Date(to), 'week')
 288     months <- format(dates, fmt)
 289     hbreaks<-unclass(as.POSIXct(dates))
 290
 291     xx<-NULL;
 292     yy<-NULL;
 293     yy_sd_high<-NULL;
 294     yy_sd_low<-NULL;
 295     date_list <- NULL;
 296
 297     x<-seq(-20,20,0.01)
 298     for ( i in seq(1,length(hbreaks)-window-1) )
 299     {
 300         print (sprintf("round %s of %s", i, length(hbreaks)-window-1))
 301         # get range from t
 302         t_sub <- t[which(t$start > hbreaks[i] & t$start<= hbreaks[i+window]),]
 303         if ( length(t_sub$start) <= 1 )  { next }
 304         # take log, then sn.mle -> h
 305         d <- (t_sub$lastreply - t_sub$start)/(60*60)    # hours
 306         d <- log(d)                                     # log(hours)
 307             # sn.mle
 308         print (sprintf("length: %s", length(d)))
 309         h<-sn.em(y=d)
 310         if ( abs(h$cp['skewness']) > 0.95 )
 311         {
 312             print(h)
 313             next    # just skip it
 314         }
 315
 316         # find dsn() using h parameters  -> y
 317         m<-dsn(x, dp=cp.to.dp(h$cp))
 318         # find max of y
 319         y_peak <- x[which(m==max(m))]
 320         # plot point date, max(y)
 321         xx<- c(xx, hbreaks[round(i+window/2)])
 322         yy<- c(yy, y_peak)
 323         yy_sd_high<- y_peak + h$cp['s.d.']
 324         yy_sd_low <- y_peak - h$cp['s.d.']
 325         date_list <- c(date_list, dates[i])
 326         # plot whisker2(x0,y0,y0_hi,y0_lo)
 327     }
 328     l<-length(months)-window-1
 329     m<- months[1:l]
 330     return (rbind(xx,yy,yy_sd_high, yy_sd_low, m))
 331 }
 332
 333 # NOTE: Try something simpler, like median of the log of ttr.
 334 #       it's going to be a lot of work to explain lsn distributions.  something
 335 #       more obvious would be a lot easier.
 336
 337 par(mfrow=c(4,1))
 338 par(mai=c(.3,0.3,0.3,0.3))
 339 for ( s in c(7)) #,14,21) )
 340 {
 341     d<- median_time_to_resolve_window(t2, "2004/1/1", "2010/2/28", s, "%b%y")
 342     plot(d[,1], exp(as.numeric(d[,5]))/24, type='l', lty=1, xlab="",
 343             axes=F, ylim=c(0.01, 15), ylab="Days to Resolve", col='orange')
 344     lines(d[,1], exp(as.numeric(d[,4]))/24, lty=1, col='red')
 345     lines(d[,1], exp(as.numeric(d[,3]))/24, lty=1, col='black')
 346     axis(1, labels=d[,7], at=d[,1])
 347     axis(2, las=1)
 348     m<-round(max(exp(as.numeric(d[,4]))/24), 2)
 349     axis(2, labels=m, at=m, las=1)
 350     abline(h=m, lty=2, col='grey40')
 351 }
 352
 353 # monitor
 354     d2<- median_time_to_resolve_window(m2, "2007/02/1", "2010/2/28", s, "%b%y")
 355     plot(d[,1], exp(as.numeric(d[,2]))/24, type='l', lty=1, xlab="",
 356             axes=F, ylim=c(0.01, 165), ylab="Days to Resolve", col='white')
 357     lines(d2[,1], exp(as.numeric(d2[,5]))/24, lty=1, col='red')
 358     lines(d2[,1], exp(as.numeric(d2[,4]))/24, lty=1, col='red')
 359     lines(d2[,1], exp(as.numeric(d2[,3]))/24, lty=1, col='black')
 360     axis(1, labels=d[,7], at=d[,1])
 361     axis(2, las=1)
 362     m<-round(max(exp(as.numeric(d2[,4]))/24), 2)
 363     axis(2, labels=m, at=m, las=1)
 364     abline(h=m, lty=2, col='grey40')
 365
 366
 367
 368 mean_time_to_resolve <- function (t, from, to, type, fmt="%b")
 369 {
 370     # find 'type' range of days
 371     dates <-seq(as.Date(from), as.Date(to), type)
 372     months <- format(dates, fmt)
 373     hbreaks<-unclass(as.POSIXct(dates))
 374
 375     xx<-NULL;
 376     yy<-NULL;
 377     yy_sd_high<-NULL;
 378     yy_sd_low<-NULL;
 379     date_list <- NULL;
 380
 381     for ( i in seq(1,length(hbreaks)-1) )
 382     {
 383         # get range from t
 384         t_sub <- t[which(t$start > hbreaks[i] & t$start<= hbreaks[i+1]),]
 385         if ( length(t_sub$start) == 0 )  { next }
 386         # take log, then sn.mle -> h
 387         d <- (t_sub$lastreply - t_sub$start)/(60*60)    # hours
 388         d <- log(d)                                     # log(hours)
 389             # sn.mle
 390         h<-sn.em(y=d)
 391         if ( abs(h$cp['skewness']) > 0.95 )
 392         {
 393             print(h)
 394         }
 395
 396         # find dsn() using h parameters  -> y
 397         x<-seq(-8,10,0.01)
 398         m<-dsn(x, dp=cp.to.dp(h$cp))
 399         # find max of y
 400         y_peak <- x[which(m==max(m))]
 401         # plot point date, max(y)
 402         xx<- c(xx, hbreaks[i])
 403         yy<- c(yy, y_peak)
 404         yy_sd_high<- y_peak + h$cp['s.d.']
 405         yy_sd_low <- y_peak - h$cp['s.d.']
 406         date_list <- c(date_list, dates[i])
 407         # plot whisker2(x0,y0,y0_hi,y0_lo)
 408     }
 409     m<- months[1:length(months)-1]
 410     return (rbind(xx,yy,yy_sd_high, yy_sd_low, m))
 411 }
 412
 413
 414 par(mfrow=c(5,1))
 415 par(mai=c(.3,0.3,0.3,0.3))
 416 for ( s in c("10 days", "2 weeks", "3 weeks", "month", "2 months"))
 417 #for ( s in c("month") )
 418 {
 419     d<- mean_time_to_resolve(t2, "2004/1/1", "2010/2/28", s, "%b%y")
 420     plot(d[1,], exp(as.numeric(d[2,]))/24, type='l', axes=F)
 421     points(d[1,], exp(as.numeric(d[2,]))/24, pch=23)
 422     axis(1, labels=d[5,], at=d[1,])
 423     axis(2)
 424 }
 425
 426
 427
 428 tstamp <-unclass(as.POSIXct("2007-05-01", origin="1960-01-01"))
 429 t_7a <- t_7[t_rep <- read.csv('rt_replies.csv', sep=',', header=TRUE)
 430 t2_rep <- t_rep[which(t_rep$complete == 1),]
 431 t2_rep <- t_rep[which(t_rep$diff != 0),]
 432
 433 which(t_7$start < tstamp),]
 434 t_7b <- t_7[which(t_7$start >= tstamp),]
 435
 436 #end_image()
 437 h4<-time_hist(t_4, 10000, "2004")
 438 h5<-time_hist(t_5, 10000, "2005")
 439 h6<-time_hist(t_6, 10000, "2006")
 440 #h7<-time_hist(t_7, 10000, "2007")
 441 h7a<-time_hist(t_7a, 10000, "2007")
 442 h7b<-time_hist(t_7b, 10000, "2007")
 443 h8<-time_hist(t_8, 10000, "2008")
 444 h9<-time_hist(t_9, 10000, "2009")
 445
 446 tstamp <-unclass(as.POSIXct("2009-09-01", origin="1960-01-01"))
 447 m_9a <- m_9[which(m_9$start < tstamp),]
 448 m_9b <- m_9[which(m_9$start >= tstamp),]
 449
 450 split_by_time <- function (t, datestr)
 451 {
 452     tstamp <-unclass(as.POSIXct(datestr, origin="1960-01-01"))
 453     a <- t[which(t$start < tstamp),]
 454     b <- t[which(t$start >= tstamp),]
 455     v<- list('before'=a, 'after'=b)
 456     return (v);
 457 }
 458
 459 mh7 <- time_hist(m_7, 10000, '2007')
 460
 461 sm_8 <- split_by_time(m_8, "2008-07-01")
 462 #mh8a <- time_hist(rbind(m_7, m_8$before, m_8$after), 10000, '2008')
 463 #mh8a <- time_hist(rbind(m_7[which(log((m_7$lastreply-m_7$start)/(60*60))>2),]), 10000, '2008')
 464 # m_7 is junk data
 465
 466 mh_8 <- time_hist(sm_8$before, 10000, '2008')
 467
 468 sm_9 <- split_by_time(m_9, "2009-09-01")
 469
 470 mh_89 <- time_hist(rbind(sm_8$after, sm_9$before), 10000, '2009')
 471 mh_9 <- time_hist(sm_9$after, 10000, '2009')
 472
 473
 474 x<-seq(-8,10,0.01)
 475 #x<- exp(x)/24
 476
 477 #my7<-dsn(x, dp=cp.to.dp(mh7$cp))
 478 my8<-dsn(x, dp=cp.to.dp(mh_8$cp))
 479 my89<-dsn(x, dp=cp.to.dp(mh_89$cp))
 480 my9<-dsn(x, dp=cp.to.dp(mh_9$cp))
 481
 482 y4<-dsn(x, dp=cp.to.dp(h4$cp))
 483 y5<-dsn(x, dp=cp.to.dp(h5$cp))
 484 y6<-dsn(x, dp=cp.to.dp(h6$cp))
 485 y7a<-dsn(x, dp=cp.to.dp(h7a$cp))
 486 y7b<-dsn(x, dp=cp.to.dp(h7b$cp))
 487 y8<-dsn(x, dp=cp.to.dp(h8$cp))
 488 y9<-dsn(x, dp=cp.to.dp(h9$cp))
 489
 490 start_image("rt_time_to_resolve.png")
 491 par(mfrow=c(1,1))
 492 par(mai=c(1.0,0.7,0.7,0.7))
 493 # monitor
 494 plot(x, my9, col='blue', type='l', axes=F, xlab="Days to Resolve", ylab="Density")
 495 axis(1, labels=c(0.0001, 0.01, 0.1, 1, 5, 20, 100), at=c(0.0001, 0.01, 0.1, 1, 5, 20, 100))
 496 axis(2)
 497 lines(x, my8, col='dodgerblue')
 498 lines(x, my7, col='turquoise')
 499 abline(v=x[which(my8==max(my8))])
 500 abline(v=x[which(my9==max(my9))])
 501
 502 # heavy
 503 lines(x, y7a, col='green3')
 504 lines(x, y4, col='green4')
 505 lines(x, y5, col='greenyellow')
 506
 507 abline(v=x[which(y4==max(y4))])
 508 abline(v=x[which(y5==max(y5))])
 509 abline(v=x[which(y7a==max(y7a))])
 510
 511 # light
 512 lines(x, y7b, col='orange', type='l')
 513 lines(x, y6, col='orange3')
 514 lines(x, y8, col='firebrick2')
 515 lines(x, y9, col='firebrick4')
 516
 517 abline(v=x[which(y7b==max(y7b))])
 518 abline(v=x[which(y6==max(y6))])
 519 abline(v=x[which(y8==max(y8))])
 520 abline(v=x[which(y9==max(y9))])
 521
 522 end_image()
 523
 524 whisker <- function (x0,y0,sd, length=0.05)
 525 {
 526     arrows(x0, y0, x0, y0+sd, code=2, angle=90, length=length)
 527     arrows(x0, y0, x0, y0-sd, code=2, angle=90, length=length)
 528 }
 529
 530 whisker2 <- function (x0,y0, y0_high, y0_low, col="black", length=0.05)
 531 {
 532     arrows(x0, y0, x0, y0_high, code=2, angle=90, length=length, col=col)
 533     arrows(x0, y0, x0, y0_low, code=2, angle=90, length=length, col=col)
 534 }
 535
 536 # NOTE: ** monthly averages might make a more compelling case than annual averages.
 537 start_image("rt_aggregate_times.png")
 538 par(mfrow=c(1,1))
 539 par(mai=c(1,1,1,1))
 540 par(mar=c(5,4,4,4))
 541
 542 s_list <- c(1519, 1596, 1112, 1591, 1019, 815)
 543 m_list <- c(0,0,0,    119,  229,  251)
 544 x_tick_list <- c(1,   2.5,  4, 5.5, 7, 8.5)
 545 x_tt_resolve_list <- c(1,   2.5,  4, 5.2,5.8, 7, 8.5)
 546 y_tt_resolve_list <- c( x[which(y4==max(y4))],
 547                             x[which(y5==max(y5))],
 548                             x[which(y6==max(y6))],
 549                             x[which(y7a==max(y7a))],
 550                             x[which(y7b==max(y7b))],
 551                             x[which(y8==max(y8))],
 552                             x[which(y9==max(y9))])
 553
 554
 555 y_mean_list <- c( h4$cp['mean'],
 556                 h5$cp['mean'],
 557                 h6$cp['mean'],
 558                 h7a$cp['mean'],
 559                 h7b$cp['mean'],
 560                 h8$cp['mean'],
 561                 h9$cp['mean'])
 562
 563 y_sd_list <- c( h4$cp['s.d.'],
 564                 h5$cp['s.d.'],
 565                 h6$cp['s.d.'],
 566                 h7a$cp['s.d.'],
 567                 h7b$cp['s.d.'],
 568                 h8$cp['s.d.'],
 569                 h9$cp['s.d.'])
 570
 571 days_tt_resolve <- exp(y_tt_resolve_list)/24
 572 days_tt_resolve_low <- exp(y_tt_resolve_list-y_sd_list)/24
 573 days_tt_resolve_high <- exp(y_tt_resolve_list+y_sd_list)/24
 574
 575
 576 my_mean_list <- c( mh_8$cp['mean'],
 577                 mh_89$cp['mean'],
 578                 mh_9$cp['mean'])
 579
 580 my_sd_list <- c( mh_8$cp['s.d.'],
 581                 mh_89$cp['s.d.'],
 582                 mh_9$cp['s.d.'])
 583
 584 mx_tt_resolve_list <- c(7, 8, 8.5)
 585 my_tt_resolve_list <- c(x[which(my8==max(my8))],
 586                         x[which(my89==max(my89))],
 587                         x[which(my9==max(my9))] )
 588
 589 mdays_tt_resolve <- exp(my_tt_resolve_list)/24
 590 mdays_tt_resolve_low <- exp(my_tt_resolve_list-my_sd_list)/24
 591 mdays_tt_resolve_high <- exp(my_tt_resolve_list+my_sd_list)/24
 592
 593
 594 days_y_sd_list <- exp(y_sd_list)/24
 595 mdays_y_sd_list <- exp(my_sd_list)/24
 596
 597 days_y_sd_list <- exp(y_sd_list)/24
 598 mdays_tt_resolve <- exp(my_tt_resolve_list)/24
 599
 600 plot(x_tt_resolve_list, days_tt_resolve, type='p', pch=c(22), axes=FALSE,
 601         log='y', ylim=c(.01,350), xlab="Year", ylab='')
 602 #points(x_tt_resolve_list, days_tt_resolve, pch=c(22))
 603
 604 lines(c(x_tt_resolve_list[1:2], x_tt_resolve_list[4]), c(days_tt_resolve[1:2], days_tt_resolve[4]), col='red')
 605 lines(c(x_tt_resolve_list[3], x_tt_resolve_list[5:7]), c(days_tt_resolve[3], days_tt_resolve[5:7]), col='green')
 606 #lines(mx_tt_resolve_list, mdays_tt_resolve)
 607 #points(mx_tt_resolve_list, mdays_tt_resolve, pch=c(24))
 608
 609 lines(mx_tt_resolve_list, mdays_tt_resolve, col='blue')
 610 points(mx_tt_resolve_list, mdays_tt_resolve, pch=c(24))
 611
 612 ticks<-c(0,0.01, 0.1, 0.5,1,2,4,7,14,21, 28, 60, 120)
 613
 614 axis(1, labels=c('2004', '2005', '2006', '2007', '2008', '2009'), at=x_tick_list)
 615 axis(2, las=1, labels=ticks, at=ticks)
 616 mtext("Days to Resolve Message", 2, line=3)
 617 #axis(2, labels=ticks, at=ticks)
 618 #for (i in 1:length(days_y_sd_list) ) {
 619 #    whisker(x_tt_resolve_list[i], days_tt_resolve[i], days_y_sd_list[i])
 620 #}
 621 #for (i in 1:length(mdays_y_sd_list) ) {
 622 #    whisker(mx_tt_resolve_list[i], mdays_tt_resolve[i], mdays_y_sd_list[i])
 623 #}
 624 for (i in c(1,2,4) ) {
 625     whisker2(x_tt_resolve_list[i], days_tt_resolve[i],
 626             days_tt_resolve_high[i], days_tt_resolve_low[i], col='red')
 627 }
 628 for (i in c(3,5,6,7) ) {
 629     whisker2(x_tt_resolve_list[i], days_tt_resolve[i],
 630             days_tt_resolve_high[i], days_tt_resolve_low[i], col='green')
 631 }
 632 for (i in 1:length(mdays_y_sd_list) ) {
 633     whisker2(mx_tt_resolve_list[i], mdays_tt_resolve[i],
 634             mdays_tt_resolve_high[i], mdays_tt_resolve_low[i], col='blue')
 635 }
 636
 637 abline(h=120,col='grey80', lty=2)
 638 abline(h=21,col='grey80', lty=2)
 639 abline(h=7,col='grey80', lty=2)
 640 abline(h=2,col='grey80', lty=2)
 641 abline(h=0.5,col='grey80', lty=2)
 642 abline(h=0.1,col='grey80', lty=2)
 643
 644 legend(1, .05,
 645         cex=0.7,
 646         legend=c("Unstable Periods", "Stable Periods", "MyOps Messages"),
 647         pch=c(22, 22, 24),
 648         col=c('red', 'green', 'blue'),
 649         lty=c(1, 1,1), merge=T)
 650 end_image()
 651 # install.packages('UsingR')
 652 require(UsingR)
 653
 654 m<-min(t_4$start)
 655 d<-data.frame(
 656     '2004'=t_4$start-m,
 657     '2005'=t_5$start-m,
 658     '2006'=t_6$start-m)
 659 simple.violinplot(d)
 660
 661 par(mfrow=c(3,3))
 662 par(mai=c(.3,.3,.3,.3))
 663 sp <- function (t)
 664 {
 665     d <- (t$lastreply-t$start)/(60*60*24)
 666     simple.violinplot(log(d))
 667 }
 668 sp(t_4)
 669 sp(t_5)
 670 sp(t_6)
 671 sp(t_7)
 672 sp(t_8)
 673 sp(t_9)
 674 sp(m_8)
 675 sp(m_89)
 676 sp(m_9)
 677
 678
 679 t3 <- add_year (t2)
 680 m3 <- add_year (m2)
 681
 682 par(mfrow=c(1,2))
 683 par(mai=c(.5,.5,.5,.5))
 684 t4<-t3[which((t3$lastreply-t3$start)/(60*60*24) < 20),]
 685 t4<-t3
 686 simple.violinplot(log((lastreply-start)/(60*60*24)) ~ year, data=t4)
 687
 688 m3[which((m3$lastreply-m3$start)< 0),]
 689 m4<-m3[which((m3$lastreply-m3$start)/(60*60*24) < 100),]
 690 simple.violinplot(log((lastreply-start)/(60*60*24)) ~ year, data=m4, log='y')
 691
 692 meanof <- function (t, year)
 693 {
 694     tx <- t[which(t$year == year),]
 695     r<-sn.em(y=log((tx$lastreply-tx$start)/(60*60*24)))
 696     return (r)
 697 }
 698
 699 t_sd <- p
 700 t_p <- c( meanof(t3,2004)$cp['mean'],
 701         meanof(t3,2005)$cp['mean'],
 702         meanof(t3,2006)$cp['mean'],
 703         meanof(t3,2007)$cp['mean'],
 704         meanof(t3,2008)$cp['mean'],
 705         meanof(t3,2009)$cp['mean'],
 706         meanof(t3,2010)$cp['mean'])
 707 points(t_p)
 708 for (i in 1:length(t_sd) ) {
 709     whisker(i, t_p[i], exp(t_sd[i]))
 710 }
 711
 712
 713
 714
 715
 716
 717 #for (i in 1:length(y_tt_resolve_list) ) {
 718 #    whisker(x_tt_resolve_list[i], scale_by*y_tt_resolve_list[i], scale_by*2)
 719 #}
 720 #for (i in 1:length(my_tt_resolve_list) ) {
 721 #    whisker(mx_tt_resolve_list[i], scale_by*my_tt_resolve_list[i], scale_by*2)
 722 #}
 723
 724 #
 725 #end_image()
 726 #par(mfrow=c(2,1))
 727 #plot_rt_hist(t_67)
 728 #plot_rt_hist(t_89)
 729 par(mfrow=c(1,1))
 730
 731
 732 # system("./parse_rt_replies.py 3> rt_replies.csv")
 733 t_rep <- read.csv('rt_replies.csv', sep=',', header=TRUE)
 734 t2_rep <- t_rep[which(t_rep$complete == 1),]
 735 t2_rep <- t_rep[which(t_rep$diff != 0),]
 736
 737 mean_diff_time <- function (t, from, to, type, fmt="%b")
 738 {
 739     # find 'type' range of days
 740     dates <-seq(as.Date(from), as.Date(to), type)
 741     months <- format(dates, fmt)
 742     hbreaks<-unclass(as.POSIXct(dates))
 743
 744     xx<-NULL;
 745     yy<-NULL;
 746     yy_sd_high<-NULL;
 747     yy_sd_low<-NULL;
 748     date_list <- NULL;
 749
 750     for ( i in seq(1,length(hbreaks)-1) )
 751     {
 752         # get range from t
 753         t_sub <- t[which(t$prev > hbreaks[i] & t$prev <= hbreaks[i+1]),]
 754         if ( length(t_sub$start) == 0 )  { next }
 755         # take log, then sn.mle -> h
 756         d <- (abs(t_sub$diff)/(60*60))
 757         d <- log(d)                                     # log(hours)
 758             # sn.mle
 759         h<-sn.em(y=d)
 760         if ( abs(h$cp['skewness']) > 0.95 )
 761         {
 762             print(h)
 763         }
 764
 765         # find dsn() using h parameters  -> y
 766         x<-seq(-8,10,0.01)
 767         m<-dsn(x, dp=cp.to.dp(h$cp))
 768         # find max of y
 769         y_peak <- x[which(m==max(m))]
 770         # plot point date, max(y)
 771         xx<- c(xx, hbreaks[i])
 772         yy<- c(yy, y_peak)
 773         yy_sd_high<- y_peak + h$cp['s.d.']
 774         yy_sd_low <- y_peak - h$cp['s.d.']
 775         date_list <- c(date_list, dates[i])
 776         # plot whisker2(x0,y0,y0_hi,y0_lo)
 777     }
 778     m<- months[1:length(months)-1]
 779     return (rbind(xx,yy,yy_sd_high, yy_sd_low, m))
 780 }
 781
 782 par(mfrow=c(5,1))
 783 par(mai=c(.3,0.3,0.3,0.3))
 784 for ( s in c("2 weeks", "3 weeks", "month", "2 months"))
 785 #for ( s in c("month") )
 786 {
 787     d<- mean_diff_time(t2_rep, "2004/1/1", "2010/2/28", s, "%b%y")
 788     plot(d[1,], exp(as.numeric(d[2,]))/24, type='l', axes=F)
 789     points(d[1,], exp(as.numeric(d[2,]))/24, pch=23)
 790     axis(1, labels=d[5,], at=d[1,])
 791     axis(2)
 792 }