diff --git a/R/wtd.stats.s b/R/wtd.stats.s index f84dea91..4558c47e 100644 --- a/R/wtd.stats.s +++ b/R/wtd.stats.s @@ -66,20 +66,21 @@ wtd.quantile <- function(x, weights=NULL, probs=c(0, .25, .5, .75, 1), weights <- weights[! i] } if(type == 'quantile') { - w <- wtd.table(x, weights, na.rm=na.rm, normwt=normwt, type='list') - x <- w$x - wts <- w$sum.of.weights - n <- sum(wts) - order <- 1 + (n - 1) * probs - low <- pmax(floor(order), 1) - high <- pmin(low + 1, n) - order <- order %% 1 - ## Find low and high order statistics - ## These are minimum values of x such that the cum. freqs >= c(low,high) - allq <- approx(cumsum(wts), x, xout=c(low,high), - method='constant', f=1, rule=2)$y - k <- length(probs) - quantiles <- (1 - order)*allq[1:k] + order*allq[-(1:k)] + sorted_xi = sort(x, index.return=TRUE) + sorted_x = sorted_xi$x + sorted_weights = weights[sorted_xi$ix] + weighted_s = c() + cum_w <- cumsum(sorted_weights) + for (i in c(1:length(sorted_weights))){ + if (i > 1){ + sk = (i-1) * sorted_weights[i] + (length(sorted_weights)-1) * cum_w[i-1]} + else{ + sk = 0 + } + weighted_s = append(weighted_s, sk) + } + quantiles <- approx(weighted_s, sorted_x, xout=probs*weighted_s[length(weighted_s)], + method='linear', f=1, rule=2)$y names(quantiles) <- nams return(quantiles) }