# HW 3 boxplot6.R JLCY, 19 Dec 2013

############################################################ 

## Compute statistics from the Historical DAta. obsstats = stats(LFmonflow)

{

    write(t(arcor), file = "arcorrs1", ncol = 12)
    write(t(armean), file = "armeans", ncol = 12)
    write(t(arstdev), file = "arstdevs", ncol = 12)
    write(t(arskw), file = "arskews", ncol = 12)
    write(t(armax), file = "armax", ncol = 12)
    write(t(armin), file = "armin", ncol = 12)

    ############################################################ 

    # obs stat

    obsmean = apply(LFmonflow, 2, mean)
    obsvar = apply(LFmonflow, 2, var)
    obsskew = apply(LFmonflow, 2, skew)
    # obslag1 -- see below
    obsmax = apply(LFmonflow, 2, max)
    obsmin = apply(LFmonflow, 2, min)

    ############################################################ 

    ## Boxplot of statistics
    par(mfrow = c(2, 3))
    boxplot(armean, main = "Monthly Mean")
    lines(obsmean, col = "red")

    boxplot(arvar, main = "Monthly Var")
    lines(obsvar, col = "red")

    boxplot(arskw, main = "Monthly Skew")
    lines(obsskew, col = "red")

    # boxplot(arcor, main='Monthly Lag-1')

    xmeans = matrix(scan("arcorrs1"), ncol = 12, byrow = T)
    n = length(xmeans[, 1])
    # the first row is the means of the original data
    xmeans1 = xmeans[2:n, ]
    xs = 1:12
    zz = boxplot(split(xmeans1, col(xmeans1)), plot = F, cex = 1)
    zz$names = rep("", length(zz$names))
    z1 = bxp(zz, ylim = range(xmeans), xlab = "", ylab = "", cex = 1)
    lines(z1, xmeans[1, ], pch = 16, col = "red")
    title(main = "Lag-1 correlation")

    obslag1 = xmeans[1, ]

    boxplot(armax, main = "Monthly Max")
    lines(obsmax, col = "red")

    boxplot(armin, main = "Monthly Min")
    lines(obsmin, col = "red")

    ############################################################ 

    # plot May
    par(mfrow = c(2, 3))

    boxplot(armean[, 5], main = "May Mean")
    par(new = T)
    plot(obsmean[5], col = "red")

    boxplot(arvar[, 5], main = "May Var")
    par(new = T)
    plot(obsvar[5], col = "red")

    boxplot(arskw[, 5], main = "May Skew")
    par(new = T)
    plot(obsskew[5], col = "red")

    boxplot(arcor[, 5], main = "May Lag-1")
    par(new = T)
    plot(obslag1[5], col = "red")

    boxplot(armax[, 5], main = "May Max")
    par(new = T)
    plot(obsmax[5], col = "red")

    boxplot(armin[, 5], main = "May Min")
    par(new = T)
    plot(obsmin[5], col = "red")

    ############################################################ 

    ## boxplot of annual mean, variance, skew and lag-1 correlation

    ymean = apply(armean, 1, sum)
    yvar = apply(arvar, 1, sum)
    yskw = apply(arskw, 1, sum)
    ycor = apply(arcor, 1, sum)
    ymax = apply(armax, 1, sum)
    ymin = apply(armin, 1, sum)

    # obs avg
    yobsmean = mean(obsmean)
    yobsvar = mean(obsvar)
    yobsskew = mean(obsskew)
    yobslag1 = mean(obslag1)
    yobsmax = mean(obsmax)
    yobsmin = mean(obsmin)

    par(mfrow = c(2, 3))
    boxplot(ymean, main = "Ann Mean")
    par(new = T)
    plot(yobsmean, col = "red")

    boxplot(yvar, main = "Ann Var")
    par(new = T)
    plot(yobsvar, col = "red")

    boxplot(yskw, main = "Ann Skew")
    par(new = T)
    plot(yobsskew, col = "red")

    boxplot(ycor, main = "Ann Lag-1")
    par(new = T)
    plot(yobslag1, col = "red")

    boxplot(ymax, main = "Ann Max")
    par(new = T)
    plot(yobsmax, col = "red")

    boxplot(ymin, main = "Ann Min")
    par(new = T)
    plot(yobsmin, col = "red")

}

## Error: 找不到物件 'arcor'