option to estimate error component

matsim-scenarios · Nov 28, 2024 · 0db7c42 · 0db7c42
1 parent 12d5415
commit 0db7c42
Show file tree

Hide file tree

Showing 4 changed files with 28 additions and 3 deletions.
diff --git 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<PlanCandidate> 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
@@ -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> 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
@@ -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<Person> 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
@@ -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: