# plotmoments.S #Written by Vadim Teverovsky # These functions use the absolute moments of a series to estimate H. # plotmoment subtracts the mean and calculates the moments. # plotmoment2 computes the moments without first subtracting the mean. If a # series has non-zero mean, the second method may give incorrect estimate, and # non-linear plot. plotmoment_function(data, nvar = 50, minnpts = 5, power1 = 0.7, power2 = 2.5, slopes = 1, plotflag = T, moment = 2) { # This is a function which gives a plot of the moments(absolute) for # various aggregation levels. It calls the function Smoments.S, which, # in turn, calls Cmoments.c. The estimate of the slope, beta, gives # an estimate for the parameter, H, via : H = (beta+moment)/moment. # # Variables: # data -- The data being used. # nvar -- The number of aggregation levels. # minnpts -- The minimum number of points to be used to # estimate the variance at any aggregation level # power1,power2-- Cut-offs for points to be used in estimating # beta. # slopes -- Flag. If 1, draw the estimated line whose # slope is beta. # plotflag -- Flag. If T, plot is performed. 12/27/94 # moment -- What absolute moment to calculate. cat("\n") # n is the length of the data vector, variances is the vector of length nvar+1 # which is the output of Smoments. It is the vector of moments at the # different levels of aggregation. n <- length(data) variances <- Smoments(data, nvar, minnpts, moment) # 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/minnpts))/nvar # lv and lm are going to be the y and x values respectively of the points in # the plot. lv <- log10(variances) lm <- rep(0, nvar + 1) # power1 and power2 determine whether the points will be included in the # estimate of beta. for (i in 1:(nvar + 1)) lm[i] <- log10(floor(10^((i - 1) * increment))) # Above line changed 2/28/95, to make the plotting consistent with # calculations. lmadj <- lm[(lm >= power1) & (lm <= power2)] lvadj <- lv[(lm >= power1) & (lm <= power2)] lmadjc <- lm[(lm < power1) | (lm > power2)] lvadjc <- lv[(lm < power1) | (lm > power2)] # The line is fitted to the appropriate points. fit <- lsfit(lmadj, lvadj) intercept <- fit$coef[1] beta <- fit$coef[2] # The plot is done here. # Changed the normalization from dividing by lv[1] to subtracting it.(Logs!) # 5/23/1994. Plotflag added 12/27/94 if (plotflag == T) { plot(c(0, lm), c(0, lv - lv[1]), xlab = "log10(m)", ylab = "log10(moments)", type = "n") points(lmadj, lvadj - lv[1], pch = 16, mkh = 0.06) points(lmadjc, lvadjc - lv[1], pch = 3, mkh = 0.06) # abline(0, - moment/2, lty = 2) if (slopes == 1) { lines(lmadj, intercept - lv[1] + (beta) * lmadj) } } # Returning results 12/27/94. cat(" beta = ", beta, "\n\n") cat("H= ", (beta + moment)/moment, "\n\n") result <- (beta + moment)/moment names(result) <- NULL return(result, variances) } Smoments function(data, nvar, minpts, moment) { n <- length(data) if(!is.loaded(C.symbol("Cmoments"))) dyn.load("Cmoments.o") .C("Cmoments", data = as.double(data), n = n, nvar = as.integer(nvar), minpts = as.integer(minpts), moment = as.double(moment), v = double(nvar + 1))$v } plotmoment2 _ function(data, nvar = 50, minnpts = 5, power1 = 0.7, power2 = 2.5, slopes = 1, plotflag = T, moment = 2) { # This is a function which gives a plot of the variances for various # aggregation levels. It calls the function Smoments2.S, which, in # turn, calls Cmoments2.c. The estimate of the slope, beta, gives # an estimate for the parameter, H, via : H = (beta+2)/2. # Variables: # # data -- The data being used. # nvar -- The number of aggregation levels. # minnpts -- The minimum number of points to be used to # estimate the variance at any aggregation level # power1,power2-- Cut-offs for points to be used in estimating # beta. # slopes -- Flag. If 1, draw the estimated line whose # slope is beta. # plotflag -- Flag. If T, plot is performed. 12/27/94 # cat("\n") # n is the length of the data vector, variances is the vector of length nvar+1 # which is the output of Svariances. It is the vector of variances at the # different levels of aggregation. n <- length(data) variances <- Smoments2(data, nvar, minnpts, moment) # 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/minnpts))/nvar # lv and lm are going to be the y and x values respectively of the points in # the plot. lv <- log10(variances) lm <- rep(0, nvar + 1) # power1 and power2 determine whether the points will be included in the # estimate of beta. for (i in 1:nvar + 1) lm[i] <- log10(floor(10^((i - 1) * increment))) # Above line changed 2/28/95, to make the plotting consistent with # calculations. lmadj <- lm[(lm >= power1) & (lm <= power2)] lvadj <- lv[(lm >= power1) & (lm <= power2)] lmadjc <- lm[(lm < power1) | (lm > power2)] lvadjc <- lv[(lm < power1) | (lm > power2)] # The line is fitted to the appropriate points. fit <- lsfit(lmadj, lvadj) intercept <- fit$coef[1] beta <- fit$coef[2] # The plot is done here. # Changed the normalization from dividing by lv[1] to subtracting it.(Logs!) # 5/23/1994. Plotflag added 12/27/94 if (plotflag == T) { plot(c(0, lm), c(0, lv - lv[1]), xlab = "log10(m)", ylab = "log10(moments)", type = "n") points(lmadj, lvadj - lv[1], pch = 16, mkh = 0.06) points(lmadjc, lvadjc - lv[1], pch = 3, mkh = 0.06) # abline(0, - moment/2, lty = 2) if (slopes == 1) { lines(lmadj, intercept - lv[1] + (beta) * lmadj) } } # Returning results 12/27/94. cat(" beta = ", beta, "\n\n") cat("H= ", (beta + moment)/moment, "\n\n") result <- (beta + moment)/moment names(result) <- NULL return(result) } Smoments2 function(data, nvar, minpts, moment) { n <- length(data) if(!is.loaded(C.symbol("Cmoments2"))) dyn.load("Cmoments2.o") .C("Cmoments2", data = as.double(data), n = n, nvar = as.integer(nvar), minpts = as.integer(minpts), moment = as.double(moment), v = double(nvar + 1))$v } # eof plotmoments.S