plotratio <- function(data, nvar = 50, power1 = 0.7, power2 = 2.5, plotflag = T) { # Variables: # data -- The data being used. # nvar -- The number of aggregation levels. # power1,power2-- Cut-offs for points to be used in estimating # beta. # Calls functions Speng and Speng.av which themselves call Cpeng.o and # Cpeng_av.o. Written by Vadim Teverovsky. cat("\n") # n is the length of the data vector #Compute the MEDIAN Variance of Residuals n <- length(data) med <- Speng(data, nvar) # The increment to be used in plotting and estimating the slope. Depends on # the number of levels, and the minimum number of points needed at each # level. increment <- (log10(n/3))/nvar # lv and lm are going to be the y and x values respectively of the points # in the plot. lr <- log10(med) lm <- rep(0, nvar) # power1 and power2 determine whether the points will be included in the # estimate of beta. for(i in 1:nvar) lm[i] <- log10(3 * floor(10^(i * increment))) lmadj <- lm[(lm >= power1) & (lm <= power2)] lradj <- lr[(lm >= power1) & (lm <= power2)] lmadjc <- lm[(lm < power1) | (lm > power2)] lradjc <- lr[(lm < power1) | (lm > power2)] #Compute the AVERAGE Variance of Residuals ave <- Speng.av(data, nvar) lv <- log10(ave) lv <- lr - lv lvadj <- lv[(lm >= power1) & (lm <= power2)] lvadjc <- lv[(lm < power1) | (lm > power2)] # The line is fitted to the appropriate points. fit <- lsfit(lmadj, lvadj - lv[1])$coef if(plotflag == T) { plot(c(0, lm), c(0, lv - lv[1]), xlab = "log(aggregation level)", type = "n", cex = 1.3) points(lmadj, lvadj - lv[1], pch = 16, mkh = 0.08) points(lmadjc, lvadjc - lv[1], pch = 3, mkh = 0.06) abline(fit) } # estimate of alpha returned. return(2/(fit[2] + 1)) } plotratio.reg <- function(data, nvar = 50, power1 = 0.7, power2 = 2.5, plotflag = T) { # Variables: # data -- The data being used. # nvar -- The number of aggregation levels. # power1,power2-- Cut-offs for points to be used in estimating # beta. cat("\n") # n is the length of the data vector, variances is the vector of length nvar+1 # which is the output of Speng1. It is the vector of variances at the # different levels of aggregation. n <- length(data) med <- Speng(data, nvar) # The increment to be used in plotting and estimating the slope. Depends on # the number of levels, and the minimum number of points needed at each # level. increment <- (log10(n/3))/nvar # lv and lm are going to be the y and x values respectively of the points # in the plot. lr <- log10(med) lm <- rep(0, nvar) # power1 and power2 determine whether the points will be included in the # estimate of beta. for(i in 1:nvar) lm[i] <- log10(3 * floor(10^(i * increment))) lmadj <- lm[(lm >= power1) & (lm <= power2)] lradj <- lr[(lm >= power1) & (lm <= power2)] lmadjc <- lm[(lm < power1) | (lm > power2)] lradjc <- lr[(lm < power1) | (lm > power2)] ave <- Speng.av(data, nvar) lv <- log10(ave) lv <- lr - lv lvadj <- lv[(lm >= power1) & (lm <= power2)] lvadjc <- lv[(lm < power1) | (lm > power2)] # The line is fitted to the appropriate points. Here is the one difference # from plotratio!!!! fit <- ltsreg(lmadj, lvadj - lv[1])$coefficients if(plotflag == T) { plot(c(0, lm), c(0, lv - lv[1]), xlab = "log(aggregation level)", type = "n", cex = 1.3) points(lmadj, lvadj - lv[1], pch = 16, mkh = 0.08) points(lmadjc, lvadjc - lv[1], pch = 3, mkh = 0.06) abline(0, 0, lty = 2) abline(fit) } return(2/(fit[2] + 1)) } Speng <- function(data, nvar) { n <- length(data) if(!is.loaded(C.symbol("Cpeng"))) dyn.load("Cpeng.o") minpts <- 1 .C("Cpeng", data = as.double(data), n = as.integer(n), nvar = as.integer(nvar), minpts = as.integer(minpts), v = double(length = nvar + 1))$v[2:(nvar + 1)] } Speng.av <- function(data, nvar) { n <- length(data) if(!is.loaded(C.symbol("Cpeng_av"))) dyn.load("Cpeng_av.o") minpts <- 1 .C("Cpeng_av", data = as.double(data), n = as.integer(n), nvar = as.integer(nvar), minpts = as.integer(minpts), v = double(length = nvar + 1))$v[2:(nvar + 1)] }