From 0db7c423aa053d09ae0b11c6388b612a56b884b8 Mon Sep 17 00:00:00 2001 From: rakow Date: Thu, 28 Nov 2024 20:59:51 +0100 Subject: [PATCH] option to estimate error component --- .../matsim/prepare/choices/ComputePlanChoices.java | 2 ++ .../run/scoring/AdvancedScoringConfigGroup.java | 4 ++++ .../org/matsim/run/scoring/PseudoRandomScorer.java | 13 ++++++++++++- .../choicemodels/estimate_biogeme_plan_choice.py | 12 ++++++++++-- 4 files changed, 28 insertions(+), 3 deletions(-) diff --git a/src/main/java/org/matsim/prepare/choices/ComputePlanChoices.java b/src/main/java/org/matsim/prepare/choices/ComputePlanChoices.java index 19a5b59c..67f7bff1 100644 --- a/src/main/java/org/matsim/prepare/choices/ComputePlanChoices.java +++ b/src/main/java/org/matsim/prepare/choices/ComputePlanChoices.java @@ -205,6 +205,7 @@ public Integer call() throws Exception { header.add("income"); header.add("util_money"); header.add("choice"); + header.add("n_trips"); for (int i = 1; i <= topK; i++) { @@ -285,6 +286,7 @@ public void run(Person person) { // choice, always the first one row.add(1); + row.add(model.trips()); List candidates = ctx.generator.generate(model, modes, null); diff --git a/src/main/java/org/matsim/run/scoring/AdvancedScoringConfigGroup.java b/src/main/java/org/matsim/run/scoring/AdvancedScoringConfigGroup.java index 6f3200ca..188cc174 100644 --- a/src/main/java/org/matsim/run/scoring/AdvancedScoringConfigGroup.java +++ b/src/main/java/org/matsim/run/scoring/AdvancedScoringConfigGroup.java @@ -33,6 +33,10 @@ public final class AdvancedScoringConfigGroup extends ReflectiveConfigGroup { @Comment("Scale for pseudo random errors. 0 disables it completely.") public double pseudoRamdomScale = 0; + @Parameter + @Comment("Distribution of the random error terms.") + public VariationType pseudoRandomDistribution = VariationType.normal; + private final List scoringParameters = new ArrayList<>(); public AdvancedScoringConfigGroup() { diff --git a/src/main/java/org/matsim/run/scoring/PseudoRandomScorer.java b/src/main/java/org/matsim/run/scoring/PseudoRandomScorer.java index f15bbd07..bfa4dfec 100644 --- a/src/main/java/org/matsim/run/scoring/PseudoRandomScorer.java +++ b/src/main/java/org/matsim/run/scoring/PseudoRandomScorer.java @@ -5,6 +5,7 @@ import org.apache.commons.math3.util.FastMath; import org.apache.commons.rng.UniformRandomProvider; import org.apache.commons.rng.core.source64.XoRoShiRo128PlusPlus; +import org.apache.commons.rng.sampling.distribution.ZigguratSampler; import org.matsim.api.core.v01.Id; import org.matsim.api.core.v01.population.Person; import org.matsim.core.config.Config; @@ -23,12 +24,14 @@ public final class PseudoRandomScorer { private final PseudoRandomTripError tripScore; private final long seed; private final double scale; + private final AdvancedScoringConfigGroup.VariationType distribution; @Inject public PseudoRandomScorer(PseudoRandomTripError tripScore, Config config) { this.tripScore = tripScore; this.seed = config.global().getRandomSeed(); this.scale = ConfigUtils.addOrGetModule(config, AdvancedScoringConfigGroup.class).pseudoRamdomScale; + this.distribution = ConfigUtils.addOrGetModule(config, AdvancedScoringConfigGroup.class).pseudoRandomDistribution; } /** @@ -47,7 +50,11 @@ public double scoreTrip(Id personId, String routingMode, String prevActi rng.nextLong(); } - return sampleGumbel(rng, 0, scale); + return switch (distribution) { + case gumbel -> sampleGumbel(rng, 0, scale); + case normal -> sampleNormal(rng, 0, scale); + default -> throw new IllegalStateException("Unsupported distribution: " + distribution); + }; } @@ -71,4 +78,8 @@ private double sampleGumbel(UniformRandomProvider rng, double mu, double beta) { return mu - FastMath.log(-FastMath.log(v)) * beta; } + private double sampleNormal(UniformRandomProvider rng, double mu, double sigma) { + return mu + ZigguratSampler.NormalizedGaussian.of(rng).sample() * sigma; + } + } diff --git a/src/main/python/choicemodels/estimate_biogeme_plan_choice.py b/src/main/python/choicemodels/estimate_biogeme_plan_choice.py index a6f7950c..2011908d 100644 --- a/src/main/python/choicemodels/estimate_biogeme_plan_choice.py +++ b/src/main/python/choicemodels/estimate_biogeme_plan_choice.py @@ -17,7 +17,7 @@ if __name__ == "__main__": parser = argparse.ArgumentParser(description="Estimate choice model for daily trip usage") parser.add_argument("--input", help="Path to the input file", type=str, - default="../../../../plan-choices-subtour_70.csv") + default="../../../plan-choices-subtour_70.csv") parser.add_argument("--mxl-modes", help="Modes to use mixed logit for", nargs="+", type=str, default=["pt", "bike", "ride", "car"]) parser.add_argument("--no-mxl", help="Disable mixed logit", action="store_true") @@ -27,6 +27,7 @@ parser.add_argument("--exp-income", help="Exponent for income", type=float, default=1) parser.add_argument("--util-money", help="Utility of money", type=float, default=1) parser.add_argument("--est-util-money", help="Estimate utility of money", action="store_true") + parser.add_argument("--est-error-component", help="Add a normal error component to each trip choice", action="store_true") parser.add_argument("--est-price-perception-car", help="Estimate price perception", action="store_true") parser.add_argument("--est-price-perception-pt", help="Estimate price perception", action="store_true") parser.add_argument("--same-price-perception", help="Only estimate one fixed price perception factor", action="store_true") @@ -110,6 +111,13 @@ print("Estimating car asc, instead of daily utility") B_CAR = 0 + if args.est_error_component: + print("Estimating error component") + B_ERROR_S = Beta('B_ERROR_S', 0.5, None, None, ESTIMATE) + B_ERROR = B_ERROR_S * bioDraws('B_ERROR_RND', 'NORMAL') + else: + B_ERROR = 0 + print("Using MXL modes", args.mxl_modes) U = {} AV = {} @@ -130,7 +138,7 @@ u += v[f"plan_{i}_car_used"] * B_CAR - U[i] = u + U[i] = u + B_ERROR * v["n_trips"] AV[i] = v[f"plan_{i}_valid"] if args.no_mxl: