A robustness release for very heavy-tailed populations, found with a 20-seed repeat of the simulation study.
method = "cobyla" on data: one of the starting points now comes from a
quick, coarse exact search (dynamic programming over about 400 quantile cut
points, a few milliseconds). On infinite-variance samples such as a Pareto
with shape 2, whose range can be hundreds of times the interquartile range,
COBYLA previously missed the narrow lower strata at larger h (up to 170%
above the optimum); it is now within 1% of the exact optimum in all our tests.
method = "cobyla": random restarts now perturb the starting probabilities
rather than the boundary values, and the minimum separation between
boundaries is a tiny fraction of the range instead of range/1000.
method = "global" on data now searches on the rank scale (boundaries are
parameterised by their position in the sorted data), so DIRECT-L covers the
lower tail as finely as the rest of the data.
A failed COBYLA solve on extreme data could stop strata.data() with an
obscure findInterval() error; it is now handled.
When the DP recursion cannot find valid boundaries for the fitted
distribution (possible on very heavy-tailed data), strata.data() now stops
with a clear message suggesting method = "exact" or method = "cobyla".
method = "dp": after the DP search on the fitted model, a short local
refinement (BOBYQA) of the boundaries is tried. It is kept only when it
lowers the model cost and gives valid strata, so it can never make the
DP solution worse. This removes most of the error caused by the grid.
method = "dp": two fixes in the original DP traceback. Grid positions
were converted to table indices by truncation, so a value such as
3.9999999 grid steps was read as step 3; they are now rounded. And the
second, finer pass only fills a window of the table around the first
answer; when its answer moved to the edge of that window, earlier stages
read empty cells and the solver stopped. It now keeps the first-pass
answer in that case (the final refinement improves it anyway).
method = "exact": when very large data are thinned to a candidate set,
the smallest and largest values are now always kept.
dpareto() and dtriang() on the search
path), so these families were never chosen and DP used a wrong model on
such data. Both are now fitted by maximum likelihood with their own code:
the Pareto (Lomax) shape and scale by Nelder-Mead and BFGS on the log
scale, and the triangle by profiling the likelihood over the mode with
the bounds set just outside the data range.strata.data() reports right-triangular data as "rtriangle" correctly,
and fitted parameters are moved back to the data scale correctly for the
Pareto and triangular families.method = "cobyla" and method = "global" on data: the objective rebuilt
the outer limits of the range as initval * maxval and finval * maxval,
which can round a hair inside the data, so the smallest or largest unit
was silently left out while optimising. On data with extreme values in
both tails (for example a Cauchy sample) the solvers then optimised the
wrong cost and could end up far from the optimum (up to 41% in a 20-seed
study). Units are now assigned with the interior boundaries only, exactly
as in the final allocation, and both solvers reach the exact optimum on
these samples. The same fix is applied in the interactive boundary plot.
method = "cobyla" on data now keeps the coarse exact start as a
candidate answer, so its result can never be worse than that start, and
the coarse search keeps more cut points in both tails.
method = "cobyla" and method = "global" on data now place each boundary
halfway between the two data values on either side of it. Before, on data
with many tied values, rounding could move the tied units at a boundary
into the next stratum, so the reported WhShTot could be slightly higher
than the cost the solver had found (0.2% on the math data at h = 3).
strata.distr() now stops with a clear message when the range
[initval, initval + dist] lies outside the support of a uniform or
triangular distribution, instead of failing inside the allocation step.
strata.distr(): the stratum population sizes Nh now add up to N
exactly (largest remainder). Before, each N * Wh was rounded on its own,
so the total could be off by one.
summary() output inside R Markdown and Quarto documents now renders as
formatted HTML. Before, the heading lines appeared as raw HTML text.
h = 4; COBYLA and GLOBAL equal to the
exact optimum on tied data; the new range check in strata.distr();
COBYLA equal to the exact optimum on a two-sided heavy-tailed sample.strata.data(..., method = "exact") runs Dynamic Programming on the data
themselves. It returns the partition with the smallest empirical Neyman cost
over all partitions into h contiguous strata with at least two units each,
never splitting tied values. It is exact for data with up to 20,000 distinct
cut points (about 15 to 20 seconds for 14,000 units at h = 6); beyond that
it searches 20,000 cut points spread evenly in rank, which can be changed
with options(stratifyR.exact_max_cand = ). The result records this in
exact_full. The "dp" solver remains the default: it is globally optimal
on a fine grid for the fitted distribution.method = "exact" places each boundary exactly on a data value. The
allocation step rebuilt it as (b / maxval) * maxval, which on some
platforms (for example Apple silicon) can round one bit below b and move
the boundary unit into the next stratum, so the reported cost could be a
little above the true optimum (0.25% for a Cauchy sample at h = 3). The
allocation step now uses a cut point halfway to the next data value.data.alloc() treated each stratum as closed at both ends, so a value equal
to a boundary was counted in two strata. Every unit is now assigned to
exactly one stratum, (b[h-1], b[h]].WhShTot is now reported to six decimals instead of three.n and never exceed N_h. Rounding each n_h on its
own could give a total of n + 1 or n - 1 (for example 63 + 103 + 135 =
301 for the anaemia Iron data at h = 3, n = 300).h boundaries in OSB, the last being the top of
the range.method = "cobyla" and method = "global" no longer force every stratum to
be at least (range)/(3h) wide. That rule excluded the narrow strata that
skewed and heavy-tailed populations need, so on Pareto, log-normal and Cauchy
data the two solvers returned designs with up to several times the optimal
Neyman cost. Boundaries now only need to be range/1000 apart, as in the
problem formulation.method = "global" no longer accepts a result whose objective is the
infeasibility penalty (a design with empty strata). If the global phase finds
no feasible point, it falls back to multi-start COBYLA.method = "global" no longer reports converged = FALSE when the BOBYQA
refinement confirms the DIRECT-L point without improving it.N argument to strata.data(), which has none. Previously this made the
H curve drop points and the generated code fail when a population size was set.method = "exact".compare_designs(), cost-aware allocation, plots and
large frames, with a worked example on household expenditure.DESCRIPTION now describes the four solvers correctly (GLOBAL refines with
BOBYQA, not COBYLA).fit_distribution() is documented as internal; it was marked for export
but never exported.cobyla.optim.R updated to match the 2.0-2 behaviour.strata.data() corrected: "dp" is the default, and the
"exact" method is documented.strata.data() and strata.distr() gain a method argument supporting
three optimisation strategies:
"dp" (default) — the original exhaustive Dynamic Programming grid search,
guaranteed globally optimal on the discrete grid."cobyla" — multi-start gradient-free local optimiser (COBYLA) via
nloptr; faster than DP for large h, suitable for exploratory work."global" — two-phase global optimiser combining DIRECT-L (deterministic
space-filling search) followed by COBYLA local refinement; closest to DP
accuracy among the two new methods.n_starts, max_iter, tol, and verbose in both functions
to control COBYLA/GLOBAL behaviour.DESCRIPTION updated: nloptr added to Imports.compare_designs() function and print method. For the same total
sample size it compares simple random sampling, proportional allocation and
Neyman allocation, with design effects, the equivalent SRS sample size and
the percentage cost saving.plot() method for "strata" objects with three types: "2d"
(histogram or density with boundaries and the fitted density formula),
"3d" and "interactive" (both need plotly).stratifyRApp() launches a self-contained shiny application:
upload a CSV or Excel file or pick a built-in dataset, compute OSB and OSS
with any solver, choose the number of strata, compare designs, explore the
boundaries, download results and generate the equivalent R code.cost = TRUE with per-stratum costs ch).summary().distr.root() and data.root()
were using the unregularized upper incomplete gamma function UGamma(s, z)
(via zipfR) where the regularized form Q(s, z) was required by the
mathematical formulation. Fixed by replacing all calls with
pgamma(z, shape = s, lower.tail = FALSE), which correctly implements
Q(s, z) = Γ(s, z) / Γ(s) as defined in Abramowitz & Stegun (1972),
Chapter 6, equation 6.5.2. Results for Weibull and Gamma populations are now
consistent with the published theory..refit_real_scale(): simultaneous MLE of all three
triangular parameters (min, max, mode) produced degenerate estimates
(min ≈ max ≈ mode). Fixed by fixing min and max to sample bounds and
estimating only the mode via MLE, matching the approach already used in
get.dist().get.dist(): triangular distribution now only selected over a tail model
(Gamma, Weibull, Log-normal, Exponential) when its AIC advantage exceeds 10
units, preventing spurious triangular wins on mildly right-skewed data.Q(r, y) (regularized) vs Γ(r, y) (unregularized)
in the Weibull and Gamma sections, with explicit R implementation note.mean/sd parameters).fix.arg
approach.NeedsCompilation field, removed stray
\url{} in Introduction.strata.data() / strata.distr() call signatures from v1.0-5 work
without modification (method defaults to "dp")."strata" S3 class and all slot names are unchanged.Initial CRAN-submission-ready version. Implemented DP solver, 10-family
distribution fitting via AIC, Neyman allocation with cost support, and
coloured summary.strata() output.