# Equivalent Quantile Function of Two Distributions Stemming from Left-Hand Threshold to Setup Conditional Probability Computations EXPERIMENTAL! This function computes the nonexceedance probability of a given quantile from a linear weighted combination of two quantile functions---a mixed distribution---when the data have been processed through the `x2xlo` function setting up left-hand thresholding and conditional probability compuation. The `par2qua2lo` function is a partial generalization of the `par2qua2` function (see there for the basic mathematics). The Examples section has an exhaustive demonstration. The resulting weighted- or mixed-quantile function is not rigorously checked for monotonic increase with $F$, which is a required property of quantile functions. However, a first-order difference on the mixed quantiles with the probabilities is computed and a warning issued if not monotonic increasing. ```r par2qua2lo(f, para1, para2, xlo1, xlo2, wfunc=NULL, weight=NULL, addouts=FALSE, inf.as.na=TRUE, ...) ``` ## Arguments - `f`: Nonexceedance probability ($0 \le F \le 1$). - `para1`: The first distribution parameters from `lmom2par` or `vec2par`. - `para2`: The second distribution parameters from `x2xlo`. - `xlo1`: The first distribution parameters from `x2xlo`. - `xlo2`: The second distribution parameters from `lmom2par` or similar. - `wfunc`: A function taking the argument `f` and computing a weight for the `para2` curve for which the complement of the computed weight is used for the weight on `para1`. - `weight`: An optional weighting argument to use in lieu of `F`. If `NULL` then weights are a function of `length(xlo1$xin)` and `length(xlo2$xin)` for the first and second distribution respectively, if `weight` has length 1, then weight on first distribution is the complement of the weight, and the weight on second distribution is `weight[1]`, and if `weight` had length 2, then `weight[1]` is the weight on the first distribution, and `weight[2]` is the weight on the second distribution. - `addouts`: In the computation of weight factors when the `xlo1$xin` and `xlo2$xin` are used by other argument settings, the `addouts` arguments triggers the inclusion of the lengths of the `xlo1$xout` and `xlo2$xout` (see source code). - `inf.as.na`: A logical controlling whether quantiles for each distribution that are non-finite are to be converted to `NA`s. If they are converter to `NA`s, then when the application of the weight or weights are made then that those indices of `NA` quantiles become a zero and the weight for the other quantile will become unity. It is suggested to review the source code. - `...`: Additional arguments to pass if needed. ## Returns The mixed quantile values for likely a subset of the provided `f` from the two distributions depending on the internals of `xlo1` and `xlo2` require the quantiles to actually start. This requires this function to return an `data.frame` that was only optional for `par2qua2`: - **f**: Nonexceedance probabilities. - **quamix**: The mixed quantiles. - **delta_curve1**: The computation `quamix` minus curve for `para1`. - **delta_curve2**: The computation `quamix` minus curve for `para2`. Alternatively, the returned value could be a weighting function for subsequent calls as `wfunc` to `par2qua2lo` (see Examples ). This alternative operation is triggered by setting `wfunc` to an arbitrary character string, and internally the contents of `xlo1` and `xlo2`, which themselves have to be called as named arguments, are recombined. This means that the `xin` and `xout` are recombined, into their respective samples. Each data point is then categorized with probability zero for the `xlo1` values and probability unity for the `xlo2` values. A logistic regression is fit using logit-link function for a binomial family using a generalized linear model. The binomial (0 or 1) is regressed as a function of the plotting positions of a sample composed of `xlo1` and `xlo2`. The coefficients of the regression are extracted, and a function created to predict the probability of event ‘xlo2’ . The `attributes` of the computed value inside the function store the coefficients, the regression model, and potentially useful for graphical review, a `data.frame` of the data used for the regression. This sounds more complicated than it really is (see source code and Examples ). ## Author(s) W.H. Asquith ## See Also `par2qua`, `par2cdf2`, `par2qua2`, `x2xlo` ## Examples ```r ## Not run: XloSNOW <- list( # data from "snow events" from prior call to x2xlo() xin=c(4670, 3210, 4400, 4380, 4350, 3380, 2950, 2880, 4100), ppin=c(0.9444444, 0.6111111, 0.8888889, 0.8333333, 0.7777778, 0.6666667, 0.5555556, 0.5000000, 0.7222222), xout=c(1750, 1610, 1750, 1460, 1950, 1000, 1110, 2600), ppout=c(0.27777778, 0.22222222, 0.33333333, 0.16666667, 0.38888889, 0.05555556, 0.11111111, 0.44444444), pp=0.4444444, thres=2600, nin=9, nout=8, n=17, source="x2xlo") # RAIN data from prior call to x2xlo() are XloRAIN <- list( # data from "rain events" from prior call to x2xlo() xin=c(5240, 6800, 5990, 4600, 5200, 6000, 4500, 4450, 4480, 4600, 3290, 6700, 10600, 7230, 9200, 6540, 13500, 4250, 5070, 6640, 6510, 3610, 6370, 5530, 4600, 6570, 6030, 7890, 8410), ppin=c(0.41935484, 0.77419355, 0.48387097, 0.25806452, 0.38709677, 0.51612903, 0.22580645, 0.16129032, 0.19354839, 0.29032258, 0.06451613, 0.74193548, 0.93548387, 0.80645161, 0.90322581, 0.64516129, 0.96774194, 0.12903226, 0.35483871, 0.70967742, 0.61290323, 0.09677419, 0.58064516, 0.45161290, 0.32258065, 0.67741935, 0.54838710, 0.83870968, 0.87096774), xout=c(1600), ppout=c(0.03225806), pp=0.03225806, thres=2599, nin=29, nout=1, n=30, source="x2xlo") QSNOW <- c(XloSNOW$xin, XloSNOW$xout ) # collect all of the snow QRAIN <- c(XloRAIN$xin, XloRAIN$xout ) # collect all of the rain PSNOW <- c(XloSNOW$ppin, XloSNOW$ppout) # probabilities collected PRAIN <- c(XloRAIN$ppin, XloRAIN$ppout) # probabilities collected # Logistic regression to blend the proportion of snow versus rain events as # ***also*** a function of nonexceedance probability wfunc <- par2qua2lo(xlo1=XloSNOW, xlo2=XloRAIN, wfunc="wfunc") # weight function # Plotting the data and the logistic regression. This shows how to gain access # to the attributes, in order to get the data, so that we can visualize the # probability mixing between the two samples. If the two samples are not a # function of probability, then each systematically would have a regression- # predicted weight of 50/50. For the RAIN and SNOW, the SNOW is likely to # produce the smaller events and RAIN the larger. opts <- par(las=1) # Note the 0.5 in the next line is arbitrary, we simply bin <- attr(wfunc(0.5), "data") # have to use wfunc() to get its attributes. FF <- seq(0,1,by=0.01); HH <- wfunc(FF); n <- length(FF) plot(bin$f, bin$prob, tcl=0.5, col=2*bin$prob+2, xlab="NONEXCEEDANCE PROBABILITY", ylab="RAIN-CAUSED EVENT RELATIVE TO SNOW") lines(c(-0.04,1.04), rep(0.5,2), col=8, lwd=0.8) # origin line at 50/50 chance text(0, 0.5, "50/50 chance line", pos=4, cex=0.8) segments(FF[1:(n-1)], HH[1:(n-1)], x1=FF[2:n], y1=HH[2:n], lwd=1+4*abs(FF-0.5), col=rgb(1-FF,0,FF)) # line grades from one color to other text(1, 0.1, "Events caused by snow", col=2, cex=0.8, pos=2) text(0, 0.9, "Events caused by rain", col=4, cex=0.8, pos=4) par(opts) # Suppose that the Pearson type III is thought applicable to the SNOW # and the AEP4 for the RAIN, now estimate respective parameters. parSNOW <- lmr2par(log10(XloSNOW$xin), type="nor" ) parRAIN <- lmr2par(log10(XloRAIN$xin), type="wak") # Two distributions are chosen to show the user than we are not constrained to one. Qall <- c(QSNOW, QRAIN) # combine into a "whole" sample XloALL <- x2xlo(Qall, leftout=2600, a=0) # apply the low-outlier threshold parALL <- lmr2par(log10(XloALL$xin), type="nor") # estimate Wakeby # Wakey has five parameters and is very flexible. FF <- nonexceeds() # useful nonexceedance probabilities col <- c(rep(0,length(QSNOW)), rep(2,length(QRAIN))) # for coloring plot(0, 0, col=2+col, ylim=c(1000,20000), xlim=qnorm(range(FF)), log="y", xlab="STANDARD NORMAL VARIATE", ylab="QUANTILE", type="n") lines(par()$usr[1:2], rep(2600, 2), col=6, lty=2, lwd=0.5) # draw threshold points(qnorm(pp(Qall, sort=FALSE)), Qall, col=2+col, lwd=0.98) # all record points(qnorm(PSNOW), QSNOW, pch=16, col=2) # snow events points(qnorm(PRAIN), QRAIN, pch=16, col=4) # rain events lines( qnorm(f2f( FF, xlo=XloSNOW)), # show fitted curve for snow events 10^par2qua(f2flo(FF, xlo=XloSNOW ), parSNOW), col=2) lines( qnorm(f2f( FF, xlo=XloRAIN)), # show fitted curve for rain events 10^par2qua(f2flo(FF, xlo=XloRAIN ), parRAIN), col=4) lines( qnorm(f2f( FF, xlo=XloALL )), # show fitted curve for all events combined 10^par2qua(f2flo(FF, xlo=XloALL ), parALL ), col=1, lty=3) PQ <- par2qua2lo( FF, parSNOW, parRAIN, XloSNOW, XloRAIN, wfunc=wfunc) lines(qnorm(PQ$f), 10^PQ$quamix, lwd=2) # draw the mixture legend(-3,20000, c("Rain curve", "Snow curve", "All combined (all open circles)", "MIXED CURVE by par2qua2lo()"), bty="n", lwd=c(1,1,1,2), lty=c(1,1,3,1), col=c(4,2,1,1)) text(-3, 15000, "A low-outlier threshold of 2,600 is used throughout.", col=6, pos=4) text(-3, 2600, "2,600", cex=0.8, col=6, pos=4) mtext("Mixed population frequency computation of snow and rainfall streamflow")# ## End(Not run) ## Not run: nsim <- 50000; FF <- runif(nsim); WF <- wfunc(FF) rB <- rbinom(nsim, 1, WF) RF <- FF[rB == 1]; SF <- FF[rB == 0] RAIN <- 10^qlmomco(f2flo(runif(length(RF)), xlo=XloRAIN), parRAIN) SNOW <- 10^qlmomco(f2flo(runif(length(SF)), xlo=XloRAIN), parSNOW) RAIN[RAIN < XloRAIN$thres] <- XloRAIN$thres SNOW[SNOW < XloSNOW$thres] <- XloSNOW$thres RAIN <- c(RAIN,rep(XloRAIN$thres, length(RF)-length(RAIN))) SNOW <- c(SNOW,rep(XloSNOW$thres, length(SF)-length(SNOW))) ALL <- c(RAIN,SNOW) lines(qnorm(pp(ALL)), sort(ALL), cex=0.6, lwd=0.8, col=3) RF <- FF[rB == 1]; SF <- FF[rB == 0] RAIN <- 10^qlmomco(RF, parRAIN) SNOW <- 10^qlmomco(SF, parSNOW) RAIN[RAIN < XloRAIN$thres] <- XloRAIN$thres SNOW[SNOW < XloSNOW$thres] <- XloSNOW$thres RAIN <- c(RAIN,rep(XloRAIN$thres, length(RF)-length(RAIN))) SNOW <- c(SNOW,rep(XloSNOW$thres, length(SF)-length(SNOW))) ALL <- c(RAIN,SNOW) lines(qnorm(pp(ALL)), sort(ALL), cex=0.6, lwd=0.8, col=3) RF <- FF[rB == 1]; SF <- FF[rB == 0] RAIN <- 10^qlmomco(f2flo(RF, xlo=XloRAIN), parRAIN) SNOW <- 10^qlmomco(f2flo(SF, xlo=XloRAIN), parSNOW) RAIN[RAIN < XloRAIN$thres] <- XloRAIN$thres SNOW[SNOW < XloSNOW$thres] <- XloSNOW$thres RAIN <- c(RAIN,rep(XloRAIN$thres, length(RF)-length(RAIN))) SNOW <- c(SNOW,rep(XloSNOW$thres, length(SF)-length(SNOW))) ALL <- c(RAIN,SNOW) lines(qnorm(pp(ALL)), sort(ALL), cex=0.6, lwd=0.8, col=3) # ## End(Not run) ```

par2qua2lo function