The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.
This vignette is for when you already know which function to call
(see vignette("RASbound-intro")) and want to understand its
result. It shows how to:
fixed_cells;x0 <- matrix(c(2, 4, 9,
6, 8, 1,
3, 5, 7), nrow = 3, byrow = TRUE)
row_targets <- c(25, 20, 25)
col_targets <- c(20, 25, 25)
upper <- matrix(c(Inf, Inf, 6,
Inf, Inf, Inf,
Inf, Inf, Inf), nrow = 3, byrow = TRUE) # cap cell [1, 3] at 6
fit <- bras_freeze_hybrid(x0, row_targets, col_targets, upper = upper)
#> Warning: `x0` had 1 cell(s) outside their own [lower, upper] range before any
#> adjustment (e.g. (1,3)); imputed to the nearest bound before solving, as if
#> that value had been supplied as the seed. If this was unintentional -- a
#> mismatch between `x0` and the bounds you meant to set -- double-check your
#> inputs; the final result may otherwise reflect a seed value you did not intend.
fit
#> <bras_fit> method: freeze_hybrid
#> converged: TRUE | feasible: TRUE | outer passes: 2
#> cells fixed at a bound: 1fixed_cellsEach row is a cell that was pinned to a bound during the fit. The
column outer_pass says when. Cells fixed in pass 1 were the
clearest violations of the unconstrained solution. Cells fixed in later
passes became binding only after earlier fixations changed the problem
that was left.
This column is your audit trail. It shows that a cell reached its bound because of the equilibrium of the free problem at that point, not because of an arbitrary early choice.
If a cell binds when you did not expect it to, look at the data first. Usually the seed value, or the target of its row or column, is not consistent with the bounds you set.
Sometimes no matrix satisfies both the targets and the bounds. By
default (on_infeasible = "error") the function stops with a
message. The same information is also available in the result:
too_tight <- matrix(1, nrow = 3, ncol = 3)
fit_bad <- suppressWarnings( # x0 is also above `upper` here; that is not the point
bras_freeze_sim(x0, row_targets, col_targets, upper = too_tight,
on_infeasible = "warning")
)
fit_bad$feasible
#> [1] FALSE
fit_bad$diagnostics$infeasible_at
#> $rows
#> [1] 1 2 3
#>
#> $cols
#> [1] 1 2 3infeasible_at$rows and infeasible_at$cols
list the margins that could not be met. Start there instead of scanning
the whole matrix. There are two common causes.
The bounds are too tight for a row or column.
The sum of upper over its cells is smaller than its target,
or the sum of lower is larger. Check it directly:
A cell that has to contribute is seeded at zero.
The next section explains this case. The fix is different: change
x0, not the bounds.
feasible = TRUE but converged = FALSEThis is not the same as proven infeasibility. It means the outer loop
reached max_outer before it found a solution or ruled one
out:
fit_slow <- bras_freeze_sim(x0, row_targets, col_targets, upper = upper,
max_outer = 1, on_infeasible = "warning")
#> Warning: `x0` had 1 cell(s) outside their own [lower, upper] range before any
#> adjustment (e.g. (1,3)); imputed to the nearest bound before solving, as if
#> that value had been supplied as the seed. If this was unintentional -- a
#> mismatch between `x0` and the bounds you meant to set -- double-check your
#> inputs; the final result may otherwise reflect a seed value you did not intend.
#> Warning: Stopped after 1 outer pass(es) without fully resolving bound
#> violations. The returned matrix may still contain cells outside their bounds;
#> consider increasing `max_outer`.
fit_slow$feasible # TRUE: not proven infeasible
#> [1] TRUE
fit_slow$converged # FALSE: it just did not finish
#> [1] FALSEHere a larger max_outer may solve it. See
?bras_fit for the full distinction.
You do not need lower = upper = 0 for a cell that must
be zero. Every variant scales cells multiplicatively, so a cell with
x0[i, j] = 0 stays exactly 0:
x0_with_zero <- x0
x0_with_zero[1, 1] <- 0
fit_zero <- bras_freeze_sim(x0_with_zero, row_targets, col_targets)
fit_zero$matrix[1, 1] # exactly 0
#> [1] 0Two points to keep in mind.
A seed outside its own [lower, upper] range
is moved to the nearest bound, with a warning. For most cells
this only saves some work. For a cell seeded at 0 with a positive
lower, it is necessary, because no iteration could move the
cell away from zero:
bad_lower <- matrix(0, 3, 3); bad_lower[1, 1] <- 5
fit_imputed <- bras_freeze_sim(x0_with_zero, row_targets, col_targets,
lower = bad_lower)
#> Warning: `x0` had 1 cell(s) outside their own [lower, upper] range before any
#> adjustment (e.g. (1,1)); imputed to the nearest bound before solving, as if
#> that value had been supplied as the seed. For cell(s) below `lower` (including
#> any seeded at exactly 0 with a positive `lower`), this substitution is what
#> makes the cell solvable at all, not just a speed optimization. If this was
#> unintentional -- a mismatch between `x0` and the bounds you meant to set --
#> double-check your inputs; the final result may otherwise reflect a seed value
#> you did not intend.
fit_imputed$matrix[1, 1] # at least 5, not 0
#> [1] 7.015922Read the warning as a prompt to check your inputs. It usually means
that x0 and the bounds do not agree.
A row or column that depends on zero-seeded cells can be
infeasible even with the loosest bounds
(lower = 0, upper = Inf), because those cells
cannot contribute anything. The feasibility check of the
bras_freeze_*() functions knows this: it treats the
achievable range of a zero-seeded cell as {0}. Such a case
is reported like any other infeasibility.
If a cell is zero in your seed only because of a data problem (rounding, suppression, a missing value), remember that RAS can never make it positive. This holds for every form of RAS, bounded or not. If the cell should be able to become positive, seed it with a small positive value instead of 0.
All functions return the same bras_fit structure, so
comparing them is direct. This is useful as a sanity check, and it is
what a simulation study would automate:
fits <- list(
freeze_sim = bras_freeze_sim(x0, row_targets, col_targets, upper = upper),
freeze_seq = bras_freeze_seq(x0, row_targets, col_targets, upper = upper),
naive = bras_naive(x0, row_targets, col_targets, upper = upper),
soft = bras_soft(x0, row_targets, col_targets, upper = upper),
cap_iter = bras_cap_iter(x0, row_targets, col_targets, upper = upper)
)
#> Warning: `x0` had 1 cell(s) outside their own [lower, upper] range before any
#> adjustment (e.g. (1,3)); imputed to the nearest bound before solving, as if
#> that value had been supplied as the seed. If this was unintentional -- a
#> mismatch between `x0` and the bounds you meant to set -- double-check your
#> inputs; the final result may otherwise reflect a seed value you did not intend.
#> Warning: `x0` had 1 cell(s) outside their own [lower, upper] range before any
#> adjustment (e.g. (1,3)); imputed to the nearest bound before solving, as if
#> that value had been supplied as the seed. If this was unintentional -- a
#> mismatch between `x0` and the bounds you meant to set -- double-check your
#> inputs; the final result may otherwise reflect a seed value you did not intend.
#> Warning: `x0` had 1 cell(s) outside their own [lower, upper] range before any
#> adjustment (e.g. (1,3)); imputed to the nearest bound before solving, as if
#> that value had been supplied as the seed. If this was unintentional -- a
#> mismatch between `x0` and the bounds you meant to set -- double-check your
#> inputs; the final result may otherwise reflect a seed value you did not intend.
#> Warning: `x0` had 1 cell(s) outside their own [lower, upper] range before any
#> adjustment (e.g. (1,3)); imputed to the nearest bound before solving, as if
#> that value had been supplied as the seed. If this was unintentional -- a
#> mismatch between `x0` and the bounds you meant to set -- double-check your
#> inputs; the final result may otherwise reflect a seed value you did not intend.
#> Warning: `x0` had 1 cell(s) outside their own [lower, upper] range before any
#> adjustment (e.g. (1,3)); imputed to the nearest bound before solving, as if
#> that value had been supplied as the seed. If this was unintentional -- a
#> mismatch between `x0` and the bounds you meant to set -- double-check your
#> inputs; the final result may otherwise reflect a seed value you did not intend.fits$calibrate_raking <- bras_calibrate_raking(x0, row_targets, col_targets, upper = upper)
#> Warning: `x0` had 1 cell(s) outside their own [lower, upper] range before any
#> adjustment (e.g. (1,3)); imputed to the nearest bound before solving, as if
#> that value had been supplied as the seed. If this was unintentional -- a
#> mismatch between `x0` and the bounds you meant to set -- double-check your
#> inputs; the final result may otherwise reflect a seed value you did not intend.
fits$calibrate_linear <- bras_calibrate_linear(x0, row_targets, col_targets, upper = upper)
#> Warning: `x0` had 1 cell(s) outside their own [lower, upper] range before any
#> adjustment (e.g. (1,3)); imputed to the nearest bound before solving, as if
#> that value had been supplied as the seed. If this was unintentional -- a
#> mismatch between `x0` and the bounds you meant to set -- double-check your
#> inputs; the final result may otherwise reflect a seed value you did not intend.
# bras_calibrate_logit() is left out on purpose. The section
# "When bras_calibrate_logit() fails" below uses this same problem.data.frame(
method = names(fits),
converged = vapply(fits, function(f) f$converged, logical(1)),
max_row_err = vapply(fits, function(f) max(f$diagnostics$row_margin_error), numeric(1)),
max_col_err = vapply(fits, function(f) max(f$diagnostics$col_margin_error), numeric(1)),
dist_from_seed = vapply(fits, function(f) sum(abs(f$matrix - x0)), numeric(1))
)
#> method converged max_row_err max_col_err
#> freeze_sim freeze_sim TRUE 7.531138e-09 3.552714e-15
#> freeze_seq freeze_seq TRUE 7.531138e-09 3.552714e-15
#> naive naive TRUE 6.700361e-09 4.934506e-09
#> soft soft TRUE 6.912888e-09 3.552714e-15
#> cap_iter cap_iter TRUE 7.531138e-09 0.000000e+00
#> calibrate_raking calibrate_raking TRUE 6.883027e-11 3.698375e-11
#> calibrate_linear calibrate_linear TRUE 7.105427e-15 2.131628e-14
#> dist_from_seed
#> freeze_sim 31.35153
#> freeze_seq 31.35153
#> naive 31.35153
#> soft 31.00000
#> cap_iter 31.35153
#> calibrate_raking 31.35153
#> calibrate_linear 32.03954dist_from_seed is the L1 distance to the seed. It is a
quick measure of how much each method changed the data. It is not a
substitute for the divergence measures of a proper simulation study.
Do not read this column as a ranking of all methods. The freeze
family, bras_naive(), bras_soft(),
bras_cap_iter() and bras_calibrate_raking()
all change the data by multiplicative scaling, like classic RAS.
bras_calibrate_linear() minimizes a least-squares distance.
A smaller dist_from_seed for it does not mean it is closer
to the right answer. It means that distance function chose a different
matrix. Which distance to minimize is a modelling choice. Use
dist_from_seed only to compare methods of the same
family.
bras_calibrate_linear() and
bras_calibrate_raking() give different matricesBoth wrappers meet the same margins, but they do not return the same
matrix. This is a reproducible case, with a fixed seed. It is the case
used in the package’s own cross-validation test
(tests/testthat/test-crossvalidation-and-properties.R):
set.seed(20260919)
x0_cmp <- matrix(round(runif(16, 1, 20)), nrow = 4)
row_targets_cmp <- round(runif(4, 80, 150))
props <- runif(4); props <- props / sum(props)
col_targets_cmp <- round(sum(row_targets_cmp) * props)
col_targets_cmp[4] <- sum(row_targets_cmp) - sum(col_targets_cmp[1:3])
fit_ras_cmp <- ras(x0_cmp, row_targets_cmp, col_targets_cmp, tol = 1e-10, max_iter = 500)
fit_raking_cmp <- bras_calibrate_raking(x0_cmp, row_targets_cmp, col_targets_cmp)
fit_linear_cmp <- bras_calibrate_linear(x0_cmp, row_targets_cmp, col_targets_cmp)
# raking is very close to our own engine (same distance function):
max(abs(fit_ras_cmp$matrix - fit_raking_cmp$matrix))
#> [1] 5.457679e-11
# linear is not, although it meets the same margins to the same tolerance:
max(abs(fit_ras_cmp$matrix - fit_linear_cmp$matrix))
#> [1] 12.12134
rbind(row_ras = rowSums(fit_ras_cmp$matrix), row_linear = rowSums(fit_linear_cmp$matrix))
#> [,1] [,2] [,3] [,4]
#> row_ras 92 121 139 135
#> row_linear 92 121 139 135bras_calibrate_raking() agrees with ras()
up to numerical tolerance. This is expected, because both minimize the
same multiplicative distance. bras_calibrate_linear() meets
the same margins, but some of its cells differ a lot from the
ras() result. This is not numerical noise: its
least-squares distance prefers another matrix.
Both matrices are correct solutions of the constraints. They differ
in which notion of “closest to the seed” selects among all valid
matrices. Which one you need depends on what the adjustment is supposed
to represent. RASbound cannot decide that for you.
bras_calibrate_logit() failsbras_calibrate_logit() was left out of the comparison
above on purpose. On this problem it fails:
# Self-contained: the running example again, so this chunk gives the same
# result in whatever order the document is run.
x0 <- matrix(c(2, 4, 9,
6, 8, 1,
3, 5, 7), nrow = 3, byrow = TRUE)
row_targets <- c(25, 20, 25)
col_targets <- c(20, 25, 25)
upper <- matrix(c(Inf, Inf, 6,
Inf, Inf, Inf,
Inf, Inf, Inf), nrow = 3, byrow = TRUE)
bras_calibrate_logit(x0, row_targets, col_targets, upper = upper)
#> Warning: `x0` had 1 cell(s) outside their own [lower, upper] range before any
#> adjustment (e.g. (1,3)); imputed to the nearest bound before solving, as if
#> that value had been supplied as the seed. If this was unintentional -- a
#> mismatch between `x0` and the bounds you meant to set -- double-check your
#> inputs; the final result may otherwise reflect a seed value you did not intend.
#> Error:
#> ! survey::calibrate(calfun = "logit") failed: Valores infinitos o faltantes en 'x'. This calibration function requires finite bounds and can be numerically fragile when most cells are only loosely bounded (see the "Logit calibration cannot represent an unbounded cell" section of ?bras_calibrate_logit). Try increasing `unbounded_mult`, switching to bras_calibrate_linear() instead -- it does not require finite bounds and does not need this substitution -- or supplying a tighter, meaningful `upper` yourself.The reason is the following. The bound upper[1, 3] = 6
is tight, and the seed x0[1, 3] = 9 is above it. The seed
is therefore moved to exactly 6 before the calibration starts. The logit
calibration needs finite bounds, so every Inf is replaced
by unbounded_mult * max(row_targets, col_targets). By
default that is 20 * 25 = 500. One cell sits on a real
bound of 6 next to cells with a bound of 500, about 80 times larger.
With such different scales, the Newton-Raphson solve of the logit
distance can diverge, which gives the error above.
In practice, the automatic replacement of Inf works well
when the bounds are on a similar scale, or when only a few cells really
matter. It does not replace real finite bounds when you know that some
cells need a tight limit. If you meet this failure on your own data, use
bras_calibrate_raking() or
bras_calibrate_linear(). Both accept Inf and
need no replacement. Or use one of the bras_freeze_*()
functions.
bras_freeze_hybrid()’s per-pass fallback is not
enoughbras_freeze_hybrid() fixes all violating cells at once.
If that would make a single pass infeasible, it falls back to sequential
fixation for that pass only. This per-pass fallback does not make it as
robust as bras_freeze_seq() in every case, as this real
input shows:
x0 <- matrix(c(29, 16, 22,
14, 4, 20,
27, 20, 23,
29, 23, 13,
26, 20, 12), nrow = 5, byrow = TRUE)
row_targets <- c(91, 45, 67, 71, 70)
col_targets <- c(133, 89, 122)
upper <- matrix(c(40, 47, 27,
17, 19, 13,
44, 16, 17,
29, 19, 41,
25, 16, 33), nrow = 5, byrow = TRUE)
fit_sim <- suppressWarnings(
bras_freeze_sim(x0, row_targets, col_targets, upper = upper, on_infeasible = "warning")
)
fit_seq <- suppressWarnings(
bras_freeze_seq(x0, row_targets, col_targets, upper = upper, on_infeasible = "warning")
)
fit_hybrid <- bras_freeze_hybrid(x0, row_targets, col_targets, upper = upper)
#> Warning: `x0` had 6 cell(s) outside their own [lower, upper] range before any
#> adjustment (e.g. (5,1), (3,2), (4,2), (5,2), (2,3), and 1 more); imputed to the
#> nearest bound before solving, as if that value had been supplied as the seed.
#> If this was unintentional -- a mismatch between `x0` and the bounds you meant
#> to set -- double-check your inputs; the final result may otherwise reflect a
#> seed value you did not intend.
#> Warning: `x0` had 6 cell(s) outside their own [lower, upper] range before any
#> adjustment (e.g. (5,1), (3,2), (4,2), (5,2), (2,3), and 1 more); imputed to the
#> nearest bound before solving, as if that value had been supplied as the seed.
#> If this was unintentional -- a mismatch between `x0` and the bounds you meant
#> to set -- double-check your inputs; the final result may otherwise reflect a
#> seed value you did not intend.
data.frame(
method = c("sim", "seq", "hybrid"),
feasible = c(fit_sim$feasible, fit_seq$feasible, fit_hybrid$feasible),
converged = c(fit_sim$converged, fit_seq$converged, fit_hybrid$converged),
n_fixed = c(nrow(fit_sim$fixed_cells), nrow(fit_seq$fixed_cells), nrow(fit_hybrid$fixed_cells))
)
#> method feasible converged n_fixed
#> 1 sim FALSE FALSE 10
#> 2 seq TRUE TRUE 5
#> 3 hybrid TRUE TRUE 5sim ends infeasible after fixing 10 cells, while
seq solves the problem with only 5. hybrid
also solves it, with 5 fixed cells. The per-pass fallback alone is not
enough here. It looks at one pass at a time, and fixing a few more cells
than necessary in several different passes can add up to an infeasible
run, even if each pass looked fine on its own.
For that reason bras_freeze_hybrid() has a second,
global safety net. If the per-pass approach ends without a solution, it
retries the whole problem from the start with sequential fixation, as
bras_freeze_seq() does, before giving up. That is how
hybrid succeeds here, and fit$method shows
that this retry was used:
A stall is a third outcome, different from infeasibility and from
running out of max_outer. It is a known limitation of
alternating-projection methods such as RAS. After enough cells are
fixed, the free cells that remain in a row, and the columns they belong
to, can reach a true fixed point that does not satisfy the margins. The
margin error stops changing at all across iterations. It is not slow
convergence.
x0_stall <- matrix(c(4, 20, 5, 49,
10, 41, 44, 35,
29, 9, 24, 18,
12, 30, 16, 46,
46, 14, 6, 10), nrow = 5, byrow = TRUE)
row_targets_stall <- c(95, 120, 63, 67, 134)
col_targets_stall <- c(146, 215, 82, 36)
upper_stall <- matrix(c(33, 31, 13, 39,
35, 41, 13, 46,
17, 34, 48, 6,
23, 53, 51, 35,
53, 57, 50, 52), nrow = 5, byrow = TRUE)
fit_stall <- suppressWarnings(
bras_freeze_sim(x0_stall, row_targets_stall, col_targets_stall,
upper = upper_stall, on_infeasible = "warning")
)
fit_stall$feasible # TRUE: not proven infeasible
#> [1] TRUE
fit_stall$converged # FALSE: the margins are not matched
#> [1] FALSE
fit_stall$diagnostics$stalled # TRUE: a fixed point, not slow progress
#> [1] TRUEThe feasibility check done before solving cannot detect this. For each row and column separately, it checks whether the bounds of the remaining free cells can reach the residual target. That condition is necessary but not sufficient. A sufficient condition has to hold jointly for every subset of rows and columns (it is related to network-flow feasibility). That check is harder and is not implemented.
What to do. A larger max_iter or
max_outer will not help, and the warning says so. A stall
often means that no matrix satisfies the targets and the bounds
together. The last section of this vignette shows such a case on real
data, with a simple test that proves it. You can loosen the bounds of
the rows and columns involved. The largest values of
fit_stall$diagnostics$row_margin_error and
col_margin_error point to them. Or you can accept some
residual mismatch, if that is tolerable for your use.
A real 55 x 55 input-output-style table (rounded to integers, and
shipped in inst/extdata as the
infeasible_*.csv files) shows a stall on real data, and how
to confirm that the problem has no solution.
read_mat <- function(f) {
as.matrix(utils::read.csv(system.file("extdata", f, package = "RASbound"),
header = FALSE))
}
x0i <- unname(read_mat("infeasible_x0.csv"))
upperi <- unname(read_mat("infeasible_upper.csv"))
rti <- read_mat("infeasible_row_targets.csv")[, 1]
cti <- read_mat("infeasible_col_targets.csv")[, 1]
fit_i <- suppressWarnings(
bras_freeze_seq(x0i, rti, cti, upper = upperi,
tol = 1e-6, max_iter = 1000, max_outer = 500,
on_infeasible = "warning")
)
fit_i$converged
#> [1] FALSE
fit_i$diagnostics$stalled
#> [1] TRUEThe simple feasibility check compares the target of each row and column with the total capacity of its cells. It passes here:
cap <- ifelse(x0i > 0, upperi, 0) # effective capacity under RAS
any(rti > rowSums(cap))
#> [1] FALSE
any(cti > colSums(cap))
#> [1] FALSEThat check is necessary but not sufficient, because the cells of a
row also have to respect the column targets. A set of cells can put at
most min(col_target, capacity) into each column. For row
13:
i <- 13
rti[i] # what the row must receive
#> [1] 11737980
sum(cap[i, ]) # its raw capacity: it looks sufficient
#> [1] 13042200
sum(pmin(cti, cap[i, ])) # what it can really absorb
#> [1] 11270963Row 13 must receive more than it can absorb once the column totals are taken into account. So no matrix satisfies this problem, and the stall is the correct diagnosis. (Solving the transportation problem as a linear program or a maximum flow gives the same answer.) If every upper limit is relaxed by 25%, the problem has a solution:
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.