# Optimal strategy for wandering robots

Consider four agents sitting on four different number lines, all starting at the origin. In this game your only move is to choose one of the four agents to move. When you do that the agent you have chosen moves left or right by 1 with prob 1/2. You must choose which agent to move at each step.

If the agent is at the origin it always moves right when you choose it.

The goal is to get one agent to 10.

The cost is the sum of the square of the total number of moves (left or right) taken by each individual agent until the first agent gets to 10. You want to minimise the expected cost

Once you have fixed a strategy you should run your code ten thousand times (each time up to the point where the first agent gets to 10) and report the average cost.

# Input and output

There is no input to this challenge. Your code need only output the mean cost over 10000 100,000 runs as specified. However please do include an explanation of your strategy for choosing which agent to move when.

# Why isn't this trivial?

It is tempting to think we should always move the agent that is closest to 10. However, if agent 1 has moved 10 steps but is at position 1 and agent 2 hasn't moved at all, it costs $$\11^2-10^2\$$ to move agent 1 and $$\1^1 = 1\$$ to move agent 2. In this case, it looks like you should move the agent that is further away from 10.

In short, the cost of moving a particular agent effectively increases each time you move it. This is a result of the scoring system. As another example, if you move agent 1 twice for your first two moves and it ends up back at the origin, it is now cheaper to move any other agent.

If two answers have the same strategy then the first answer wins. If two different strategies give similar mean costs then they should be rerun with 100,000 a million (or even more) iterations to see which is better.

A very simple strategy gives a mean cost of 6369. So the first hurdle is to do better than that.

# Winners

This question was first brilliantly solved optimally by Anders Kaseorg using a bespoke hand crafted method. usul subsequently explained how a variant of a Gittins index can in principle solve it optimally and user1502040 then showed how to do this in practice simply (sped up in the comments by Anders Kaseorg).

Before this a number of excellent almost optimal solutions were given.

• What information can the "strategy" use?
– att
Aug 8 at 18:42
• @att The positions of all the agents at all times in the past and how many moves they have made and the scoring function. What else are you thinking of? You can train a model on as many examples are you choose. Aug 8 at 18:51
• @AndersKaseorg are you sure? I mean if I compute it by simulation it's pretty consistent. When you say it diverges do you mean for a fixed strategy it has no mean (like the Cauchy distribution)? Aug 8 at 19:16
• Aha, I missed the rule “If the agent is at the origin it always moves right when you choose it”. (Probably because it contradicts the earlier rule “the agent you have chosen moves left or right by 1 with prob 1/2”.) That changes things. Aug 8 at 19:20
• @AndersKaseorg Equivalently, we want to get an agent to ±10 with ±1 steps
– att
Aug 8 at 19:22

# Ruby, modified minimum product, average cost ~4950

f=->a,b{(0..3).min_by{|q|(9-a[q])**3*b[q]}}


Try it online!

Always select the agent that minimizes (9-position)^3*(number of moves):

• On the first round, every agent makes a move.
• By using 9-a instead of 10-a, any agent in position 9 will try to move on the next iteration.
• Position seems to be more important than the cost projection, so after trying out different values, I found out that the third power is the best way to represent it.

The easiest strategy (round-robin) has a cost of about 7000.

If I let it run in a continuous loop, after a while the average stabilizes around 4950.

• Thank you for this first answer! Now you are the one to beat :) Aug 8 at 19:23
• This is a very good score already Aug 8 at 19:34
• I tried fractional powers too, but anything between 2 and 3 seems to yield the same result.
– G B
Aug 9 at 10:07
• The simplest solution with very good score (only beaten by Andres' latest answer) Aug 10 at 3:22
• (My most recent runs with fifty million attempts each say 4947.5 ± 2.3 versus 4942 ± 2.3, with 95% confidence.) Aug 10 at 13:11

# Rust, cost = 4928.296960314256

Computed exactly (modulo floating point imprecision) with a non-random algorithm, but also measured as $$\4923.8 ± 2.6\$$ ($$\1σ\$$ confidence) with $$\10^7\$$ random runs of the actual game as a check.

### Strategy

Consider the following equivalent but wetter version of the game. There’s a tower of cells $$\(p, m)\$$ with eleven columns $$\0 ≤ p ≤ 10\$$ and infinitely many rows $$\m ≥ 0\$$ going downwards. We’ll call the first ten columns $$\p = 0, \dotsc, 9\$$ “live” and the last column $$\p = 10\$$ “dead”. Each live cell has a plugged drain in the floor, connected by equal-width pipes to its diagonal neighbors below. We start with $$\1\$$ liter of water in $$\(0, 0)\$$, with $$\4\$$ grains of sand floating in the water. Our goal is to control the plugs until one of the grains ends up in the dead column, while minimizing the expected sum of the squared row numbers of the grains at that time.

The reason to formulate the game like this is that our strategy does not need to look at the grains. Instead, we’re going to open the plugs in some fixed order, letting the grains float where they may, until 100% of the water flows into the dead column in the limit as all of the plugs are open. (We’re very good at opening plugs.)

I’ll explain in a moment how to decide which order we’ll use, but to motivate that, I first need to explain how to evaluate our expected score. Define “time $$\t\$$” to be the time when $$\t\$$ liters of water have flowed into the dead column. Then, because the grains of sand flow equivalently to water molecules, the distribution of time quadruples $$\(t_1, \dotsc, t_4)\$$ when the four grains enter the dead column is just the uniform distribution over $$\[0, 1]^4\$$. Therefore, the time $$\t_\min = \min \{t_1, \dotsc, t_4\}\$$ when the score should be evaluated has a beta distribution $$\t_\min \sim \mathrm{Beta}(1, 4)\$$, whose PDF is $$\4(1 - t)^3\$$.

Let $$\q_{(p, m)}(t)\$$ be the amount of water in cell $$\(p, m)\$$ at time $$\t\$$; let $$\c_{(p, m)}(t) = m^2q_{(p, m)}(t)\$$ be its weighted cost, and $$\c_{\textrm{live}}(t)\$$ (resp. $$\c_{\textrm{dead}}(t)\$$) the sum over all live (resp. dead) cells. Then we can write our expected score as

$$\begin{split} &\int_0^1 4(1 - t)^3 \left(3\frac{c_{\textrm{live}}(t)}{1 - t} + c'_{\textrm{dead}}(t)\right)\,dt \\ &= \int_0^1 12(1 - t)^2 c_{\textrm{live}}(t)\,dt + \int_0^1 4(1 - t)^3 c'_{\textrm{dead}}(t)\,dt \\ &= \Bigl[-4(1 - t)^3 c_{\textrm{live}}(t)\Bigr]_0^1 - \int_0^1 -4(1 - t)^3 c'_{\textrm{live}}(t)\,dt + \int_0^1 4(1 - t)^3 c'_{\textrm{dead}}(t)\,dt \\ &= \int_0^1 4(1 - t)^3 (c'_{\textrm{live}}(t) + c'_{\textrm{dead}}(t))\,dt. \end{split}$$

Therefore, we see that our goal at each step is to minimize the total growth rate of the weighted cost of all cells ($$\c'_{\textrm{live}}(t) + c'_{\textrm{dead}}(t)\$$), relative to the rate of flow into the dead column (i.e. growth of $$\t\$$ itself). $$\4(1 - t)^3\$$ is decreasing, so we prefer to have smaller ratios earlier. This makes our strategy obvious: at each step, we’ll open the plug that minimizes this ratio, wait until all water has flowed out of that cell, and leave it open while proceeding to the following cells.

Since $$\c'_{\textrm{live}}(t) + c'_{\textrm{dead}}(t)\$$ is piecewise constant, the integral evaluates to an (infinite) sum, one that converges much more quickly than the Monte Carlo simulation. That’s how I determined the precise answer above.

### Code

Build and run with cargo run --release.

Cargo.toml

[package]
name = "wandering"
version = "0.1.0"
edition = "2021"

[dependencies]
rand = "0.8.5"


src/main.rs

use rand::prelude::*;
use rand::rngs::StdRng;
use std::collections::HashMap;

#[derive(Default, Debug)]
struct Rate {
cost: f64,
live: HashMap<(u32, u32), f64>,
}

fn add(map: &mut HashMap<(u32, u32), f64>, state: (u32, u32), value: f64) {
if value != 0.0 {
*map.entry(state).or_insert(0.0) += value;
}
}

struct SortingHat {
distribution: HashMap<(u32, u32), f64>,
moves_frontier: [u32; 10],
rates: HashMap<(u32, u32), Rate>,
orders: HashMap<(u32, u32), usize>,
}

impl Default for SortingHat {
fn default() -> SortingHat {
let mut distribution = HashMap::new();
distribution.insert((0, 0), 1.0);
SortingHat {
distribution,
moves_frontier: [0_u32, 1, 2, 3, 4, 5, 6, 7, 8, 9],
rates: HashMap::new(),
orders: HashMap::new(),
}
}
}

impl SortingHat {
fn get_order(&mut self, query_position: u32, query_moves: u32) -> usize {
if let Some(&order) = self.orders.get(&(query_position, query_moves)) {
return order;
}

loop {
let (_, position, moves) = (0_u32..10)
.zip(&self.moves_frontier)
.filter(|&(position, &moves)| {
position + 1 == 10 || self.moves_frontier[position as usize + 1] > moves + 1
})
.map(|(position, &moves)| {
let mut cost = 0.0;
let mut dead = 0.0;
for new_position in [position.abs_diff(1), position + 1] {
cost += moves as f64 + 1. / 2.;
if new_position == 10 {
dead += 1. / 2.;
} else if let Some(new_rate) = self.rates.get(&(new_position, moves + 1)) {
cost += new_rate.cost / 2.;
}
}
(cost / dead, position, moves)
})
.min_by(|a, b| a.partial_cmp(&b).unwrap())
.unwrap();

self.orders.insert((position, moves), self.orders.len());

let mut rate = Rate::default();
for new_position in [position.abs_diff(1), position + 1] {
rate.cost += moves as f64 + 1. / 2.;
if new_position == 10 {
rate.dead += 1. / 2.;
} else if let Some(new_rate) = self.rates.get(&(new_position, moves + 1)) {
rate.cost += new_rate.cost / 2.;
for (&newer_state, live_rate) in &new_rate.live {
}
} else {
add(&mut rate.live, (new_position, moves + 1), 1. / 2.);
}
}

if let Some(mass) = self.distribution.remove(&(position, moves)) {
for (&new_state, live_rate) in &rate.live {
add(&mut self.distribution, new_state, live_rate * mass);
}
}

for old_rate in self.rates.values_mut() {
if let Some(live_rate) = old_rate.live.remove(&(position, moves)) {
old_rate.cost += rate.cost * live_rate;
for (&new_state, new_live_rate) in &rate.live {
*old_rate.live.entry(new_state).or_insert(0.0) += new_live_rate * live_rate;
}
}
}

self.rates.remove(&(position + 1, moves + 1));
self.rates.insert((position, moves), rate);
self.moves_frontier[position as usize] += 2;

if (position, moves) == (query_position, query_moves) {
return self.orders.len() - 1;
}
}
}

fn strategy(&mut self, positions: &[u32; 4], num_moves: &[u32; 4]) -> usize {
(0..4)
.min_by_key(|&agent| self.get_order(positions[agent], num_moves[agent]))
.unwrap()
}
}

fn test(mut rng: impl Rng, mut strategy: impl FnMut(&[u32; 4], &[u32; 4]) -> usize, runs: u32) {
let mut count = 0;
let mut mean = 0.0;
let mut m2 = 0.0;
let mut mean_stdev = f64::NAN;

for _ in 0..runs {
let mut positions = [0; 4];
let mut num_moves = [0; 4];
loop {
let agent = strategy(&positions, &num_moves);
if positions[agent] == 0 || rng.gen() {
positions[agent] += 1;
} else {
positions[agent] -= 1;
}
num_moves[agent] += 1;
if positions[agent] == 10 {
break;
}
}
let sample: f64 = num_moves.into_iter().map(|m| (m as f64).powi(2)).sum();

count += 1;
let delta = sample - mean;
mean += delta / count as f64;
m2 += delta * delta * (1. - 1. / count as f64);
mean_stdev = (m2 / (count as f64 * (count - 1) as f64)).sqrt();
if count % 100 == 0 {
eprint!("\r\x1b[K{count}: {mean} ± {mean_stdev}");
}
}
eprint!("\r\x1b[K");
println!("{count}: {mean} ± {mean_stdev}");
}

fn main() {
let rng = StdRng::seed_from_u64(0);
let mut sorting_hat = SortingHat::default();
test(
rng,
|positions, num_moves| sorting_hat.strategy(positions, num_moves),
10000000,
)
}

• That's a great improvement! What is your new strategy? Aug 9 at 15:10
• @Anders Kaseorg I believe a better estimate of your program's performance is 4927.3 ± 0.67. Here's a version of your code with variance-reduced evaluation: controlc.com/3a765461a Aug 9 at 15:48
• @graffe If you have a function v(s) which estimates the future cost of a state, you can use the sum over all timesteps of v after the coin flip minus the expected v before the coin flip as a control variate. The same idea has been used successfully in computer poker. Aug 9 at 18:06
• I also have a completely non-random analysis that converges to an exact answer of 4928.296960314256 ± 0.00000000001. (It’s pretty complicated but I hope to add an explanation soon.) Aug 9 at 20:26
• @graffe I’ve added an attempt to explain my strategy. (Since I designed this before hearing about Gittins theory, it’s possible I’ve reinvented machinery that might be explained in more detail elsewhere?) Aug 11 at 10:31

# Python 3 (PyPy), 4957.6 +/- 2.6 (10M trials)

def expco(positions, num_moves,depth=0) -> float:
if depth==0:
return min((10**2 - positions[i] ** 2)
* (5 * 10**2 - positions[i] ** 2 + 6 * num_moves[i] - 2)/3 for i in range(4))
return (10 not in positions)and min(
expco([abs(p-(i==j)) for j,p in enumerate(positions)],
[m+(i==j) for j,m in enumerate(num_moves)],depth-1)/2 +
expco([p+(i==j) for j,p in enumerate(positions)],
[m+(i==j) for j,m in enumerate(num_moves)],depth-1)/2 +
num_moves[i]*2+1 for i in range(4))
def strategy(positions, num_moves, depth=1) -> int:
return min(range(4), key = lambda i:
expco([abs(p-(i==j)) for j,p in enumerate(positions)],
[m+(i==j) for j,m in enumerate(num_moves)],depth)/2 +
expco([p+(i==j) for j,p in enumerate(positions)],
[m+(i==j) for j,m in enumerate(num_moves)],depth)/2 +
num_moves[i]*2+1)



Try it online!

This is a simple brute force search on top of @Anders Kaseorg's formula (and code).

The score given was obtained with search depth 1, i.e. for each agent compute the expected cost if we can switch agents once immediately after the first move and take the best.

Apologies for the messy code.

• If the mean is really below 5000 that is very impressive. Can you try it for a lot of iterations to check? Aug 9 at 7:10
• @graffe I think I have now enough samples to say we are well below 5000. At search depth 1. Aug 9 at 9:05
• That is very cool! Aug 9 at 9:13
• @loopywalt I believe your strategy achieves 4961 +/- 0.5, if you want a more accurate estimate. Aug 9 at 13:04
• Nice one @user1502040. I take it you have a fast implementation? How many samples is your estimate based on? And would it be possible to try larger depths? (My code is too slow for that.) Aug 9 at 13:19

# Python 3, cost ≈ 5034.3 ± 2.7

$$\1σ\$$ confidence, evaluated with $$\10^7\$$ runs.

def strategy(positions, num_moves):
return min(
range(4),
key=lambda i: (10**2 - positions[i] ** 2)
* (5 * 10**2 - positions[i] ** 2 + 6 * num_moves[i] - 2),
)


Try it online! ($$\10^5\$$ runs)

### How it works

If we were to pick one agent at position $$\p\$$ that’s previously moved $$\m\$$ times, and use it for the rest of the game, the expected additional cost would be exactly

$$f(p, m) = \frac13(10^2 - p^2)(5⋅10^2 - p^2 + 6m - 2).$$

(You can verify that

$$m^2 + f(p, m) = (m + 1)^2 + \frac12(f(\lvert p - 1\rvert, m + 1) + f(p + 1, m + 1))$$

and $$\f(10, m) = 0\$$.)

At each step, we pick the agent that minimizes this quantity. (We don’t restrict ourselves to using it for the rest of the game.)

• Thank you! I wonder if yours has the same expected cost as @G.B.'s. Aug 8 at 20:08
• @graffe I got 5122.2 ± 8.4 for G B’s strategy. Aug 8 at 20:10
• Something is wrong with the typesetting in the "You can verify" part. The "and" comes directly after (m+1)^2. This is on my phone. It might be ok on a desktop though so out of your control. Aug 8 at 20:26
• @graffe This should help (switched that equation to display mode). Aug 8 at 20:28
• @att From $(m + n)^2 = m^2 + 2mn + n^2$, I knew the answer would have the form $f(p, m) = 2ma_p + b_p$ where $a_p$ is the expected number of remaining moves and $b_p$ is the expected square of the number of remaining moves. I solved the resulting 20 linear equations in $a_0, b_0, \dotsc, a_9, b_9$, observed the quadratic pattern in $a_p$ and quartic pattern in $b_p$, and verified correctness as explained in the answer. Aug 8 at 20:59

# 5050? ± 26.64

def strategy(positions, num_moves):
E = [100, 99, 96, 91, 84, 75, 64, 51, 36, 19]
_, n = min(((E[p] + m) ** 2 - m ** 2, i) for p, m, i in zip(positions, num_moves, range(4)))
return n


Try it online!

Use minimal $$\\left(E(P_i)+M_i\right)^2-M_i^2\$$, where $$\E(X)\$$ is a pre-calculated array $$\left[E(0),\dots,E(9)\right]=\left[100, 99, 96, 91, 84, 75, 64, 51, 36, 19\right]$$

generated by following codes

syms E0 E1 E2 E3 E4 E5 E6 E7 E8 E9 EA
solve([E0-E1-1, E1-(E0+E2)/2-1, E2-(E1+E3)/2-1, E3-(E2+E4)/2-1, E4-(E3+E5)/2-1, E5-(E4+E6)/2-1, E6-(E5+E7)/2-1, E7-(E6+E8)/2-1, E8-(E7+E9)/2-1, E9-(E8+EA)/2-1, EA])


Try it online!

The driver code is used from other posts.

# Python 3, cost ≈ 4928.4 ± 1

I believe this is essentially optimal. Thank you to @usul for the inspiration.

import functools

f = lambda p, m: (1. / 3) * (100 - p * p) * (500 - p * p + 6 * m - 2)

@functools.lru_cache(None)
def gittins_value(p, m, l):
if m >= 400:
return f(p, m)
if p == 0:
v_next = gittins_value(p + 1, m + 1, l)
else:
v_next = 0.5 * gittins_value(p - 1, m + 1, l)
if p != 9:
v_next += 0.5 * gittins_value(p + 1, m + 1, l)
return min((2 * m + 1) + v_next, l)

@functools.lru_cache(None)
def gittins_index(p, m):
low = 2 * m + 1
high = f(p, m)
while True:
mid = 0.5 * (low + high)
v = gittins_value(p, m, mid)
if v < mid:
high = mid
else:
low = mid
if abs(high - low) <= 1e-6:
break
return 0.5 * (low + high)

def strategy(p, m):
indices = [gittins_index(p_i, m_i) for p_i, m_i in zip(p, m)]
return indices.index(min(indices))

• Your score is essentially identical to @AndersKaseorg's but your code is written in a language I can understand. Would you be able to give an explanation or references for your answer please? Aug 11 at 3:05
• Very nice - how did you get such an accurate estimate of the score / how many trials were you able to perform?
– usul
Aug 11 at 4:09
• Would you be able to add the testing code too? I think everyone would be interested to see it. Aug 11 at 8:39
• @graffe, my interpretation: p is location, m is number of steps, l is the outside option in the 'magic game' of my post. For each l, there is a value of playing the magic game with outside option l, which is gittins_value(p,m,l). The Gittins index is the value of l where you break even, which gittins_index(p,m) finds by binary search. Finally, f is a heuristic for when the number of steps gets too high.
– usul
Aug 11 at 11:18
• Here’s a faster and more accurate non-recursive implementation without the step limit, using a Newton’s method search. Aug 12 at 0:51

# R / Rcpp: 5006 ± 7.31, or 4932 ± 3.63

Score computed by taking the average of 20 million simulated games. Bounds are a 95% bootstrap confidence interval with 10,000 iterations.

Thanks to Anders Kaseorg for pointing out a typo in my code that caused the solution to accidentally be optimal. Reader, may all your future bugs behave thusly. The following description explains the 5006 solution with a note for how to get the 4932 solution.

## Strategy:

Let's consider a simplified version of the environment in order to construct a risk heuristic. We'd like to be able to treat n consecutive movements of an agent as independent, but that won't be true if the agent ever reaches position 0 or 10 during those movements.

Suppose instead that each of the number lines is infinite in both directions, and that the probability of moving right is always 0.5 (even at position 0). In this simplified case, if we move a single agent n consecutive times, we may treat the resulting k number of movements to the right as binomially distributed.

At each step, let n = 10 - x, where x is the agent's position on the number line. Since the agent's movements are independent in the simplified environment, we can calculate the probability that the agent has moved backwards given k moves to the right as

$$Prob = 1 - \sum_{k=\lceil \frac{n}{2} \rceil}^n \dbinom n k \cdot 0.5^k \cdot (1-0.5)^{n-k}$$

Note: If you incorrectly move the one-minus inside the sum, the solution becomes near optimal for reasons that are unclear.

Let the cost function be the increase in cost accrued by moving the agent n times:

$$Cost = (x + n)^2 - (x)^2$$

Then define the risk of moving an agent as

$$Risk = Prob \cdot Cost$$

and at each step move the agent with the lowest risk.

## Intuition:

The idea is to calculate the risk that over the minimum number of steps an agent would need to reach the goal, the agent ends up moving backwards instead of forwards, i.e. We try to maximize the chance that the chosen agent will move closer to the goal. Minimizing this risk heuristic appears to perform fairly well even though the method does not take into account the effect of the order of the movements, nor the effect of reaching the zero position.

#include <Rcpp.h>
using namespace Rcpp;

//* Initialize Environment
// [[Rcpp::export]]
NumericMatrix initEnv() {
return(NumericMatrix (4, 2));
}

//* Update Environment
// [[Rcpp::export]]
NumericMatrix updateEnv(NumericMatrix x, int a) {

//Update move counter
x(a, 1) = x(a, 1) + 1;

//Update position
if(x(a, 0) == 0){
x(a, 0) = 1;
} else {
RNGScope scope;
NumericVector r = Rcpp::runif(1, 0.0, 1.0);
if(r <= 0.5){
x(a, 0) = x(a, 0) - 1;
} else {
x(a, 0) = x(a, 0) + 1;
}
}

//Return
return(x);

}

//* Calculate Loss
// [[Rcpp::export]]
double envCost(NumericMatrix x) {
double lambda = 0;
for(int i = 0; i < 4; i++){
lambda = lambda + pow(x(i, 1), 2);
}
return(lambda);
}

//* Calculate Risk
// [[Rcpp::export]]
double agentRisk(NumericMatrix x, int a) {

//Parameters
int pos = x(a, 0);
int n = 10 - pos;

//Calculate risk correctly
double prob = 0;
double cost = pow(x(a, 1) + n, 2) - pow(x(a, 1), 2);
for(int k = ceil(n / 2); k <= n; k++){
prob = prob + (Rf_choose(n, k) * pow(0.5, k) * pow(0.5, n - k));
}
prob = 1 - prob;

//Calculate risk with typo
//double prob = 0;
//double cost = pow(x(a, 1) + n, 2) - pow(x(a, 1), 2);
//for(int k = ceil(n / 2); k <= n; k++){
//  prob = prob + (1 - Rf_choose(n, k) * pow(0.5, k) * pow(0.5, n - k));
//}

//Return
return(prob * cost);

}

//* Simulate Game
// [[Rcpp::export]]
List simGame() {

//Generate environment
NumericMatrix en = initEnv();

//Complete game
bool done = false;
while(!done){

//Parameters
int a = 0;
double lambda = agentRisk(en, 0);

//Calculate risks and choose action
for(int i = 1; i < 4; i++){
double lambda_new = agentRisk(en, i);
if(lambda_new < lambda){
a = i;
lambda = lambda_new;
}
}

//Update game
en = updateEnv(en, a);

//Check victory conditions
for(int i = 0; i < 4; i++){
if(en(i, 0) >= 10){
done = true;
break;
}
}

}

//Return
return(
List::create(
Named("state") = en,
Named("cost") = envCost(en)
)
);

}

• This is a surprising and yet very good solution! Aug 12 at 7:33
• You have a bug in prob = prob + (1 - Rf_choose(n, k) * pow(0.5, k) * pow(0.5, n - k)): the 1 - is inside the loop, while your mathematical description has it outside. Surprisingly, this bug seems to have accidentally improved your score significantly. Aug 12 at 7:35
• Do you have a typo where you define n = 10 - x and then have "Cost = (x + n)^2 - (x)^2"? This doesn't seem right as x + n = 10 if n = 10 - x. Aug 12 at 7:55
• @graffe Thank you! The cost is intentionally written to include the penalty the agent accrues even during the move that takes it to the goal (because it might move backwards), so the minimum number of moves n to the terminal state from position x should be 10 - x Aug 12 at 11:52
• @AndersKaseorg Now having thought about this another hour, I think you're correct. I'm not convinced I know why this method worked. Since this is my first post, I'm not sure how these situations are resolved. Should I delete my answer, or update the description to point out my mistake? A solution with no / an incorrect explanation doesn't seem to have much utility. Aug 12 at 13:53

# Rust, cost = ~5032 (1M iterations)

const depth: usize = 1200;
const iterations: usize = 1_000_000;

lazy_static! {
static ref EXPECTED_SCORES: [[usize; 10]; depth] = {
let mut expected_scores = [[1_000_000usize; 10]; depth];

for i in (0..depth-1).rev() {
for j in 0..10 {
if j == 0 {
expected_scores[i][j] = expected_scores[i + 1][j + 1].saturating_add(2 * i + 1);
} else if j == 9 {
expected_scores[i][j] =
(expected_scores[i + 1][j - 1] / 2).saturating_add(2 * i + 1);
} else {
expected_scores[i][j] = (expected_scores[i + 1][j + 1]/ 2)
.saturating_add(expected_scores[i + 1][j - 1] / 2)
.saturating_add(2 * i + 1);
}
}
}

expected_scores
};
}

#[derive(Copy, Clone, Eq, PartialEq, Debug)]
struct Agent {
position: u8,
actions: u32,
}

fn strategy(agents: [Agent; 4]) -> usize {
return (0..4)
.min_by_key(|i| EXPECTED_SCORES[(agents[*i].actions as usize).min(depth - 1)][agents[*i].position as usize])
.unwrap();
}


Basically pre-computes a table of expected scores then uses that.

Playground link Need to reduce iterations if you don't want it to time out.

• What are the expected scores it precomputes? Aug 8 at 20:22
• This computes an approximation (due to integer truncation and depth limit) to the same values that I computed exactly in closed form. Aug 8 at 20:40
• I found the score didn't really change at all when varying the depth limit by several powers of 2 in either direciton Aug 8 at 21:00

# Python 3 (PyPy), 5027.1 ± 8.3

def strategy(positions, num_moves) -> int:
return min(
range(4),
key=lambda i: (10 - positions[i] + 2*num_moves[i]) * (10 - positions[i]),
)


Try it online!

Copies the function/testing template from Anders Keorg's answer.

Chooses the agent with the minimal best-case cost $$\(10-\mathit{position}+\mathit{moves})^2-\mathit{moves}^2=(10-\mathit{position}+2\mathit{moves})(10-\mathit{position})\$$.

# C++: ~4953 (100 million iterations)

#include <bits/stdc++.h>

using namespace std;

int rand(int a, int b) { return uniform_int_distribution<int>(a, b)(rng);}

int main() {
const int NRIT = 100'000'000;
auto go = [&] {
int l = { 10,10,10,10 }, t = {}, sol = 0;
auto cost = [&](int nr) {return l[nr] * (4 * t[nr] * l[nr] + 6); };
while (1) {
pair<int, int> best = { INT_MAX, -1 };
for (int i = 0; i < 4; ++i)
best = min(best, { cost(i), i });
int opt = best.second;
sol += 2 * t[opt]++ + 1;
if (l[opt] == 10)
--l[opt];
else {
if (rand(0,1))
--l[opt];
else
++l[opt];
}
if (!l[opt])
return sol;
}
};
long long total = 0;
for (int i = 0; i < NRIT; ++i)
total += go();
cout << (double)(total) / NRIT;

}


I misread the question at first and thought it had a 1/2 chance to move backwards or stay still. This chooses the agent with the least expected cost to finish, according to my calculations and the (wrong) assumption it can't go backwards.

$$\ \text{cost} = (10-p)(4 \cdot (10-p) \cdot m+6)\$$

I want to add 10 000 iterations is definitely not enough, I could cherrypick 10k iterations to get 4790, which is far from the real average. Standardised testing also encounters the issue of different number generators, which can actually make a small difference (it did from my testing).

• Welcome to Code Golf, and nice answer! Aug 9 at 16:15
• Thank you. Yes 10000 turns out to be far too little Aug 9 at 16:59

# Python: ~4973

TL;DR By Gittins theory, there is an optimal solution of the following form, but it has some complicated replacement for approx_grades and this is a heuristic.

import random,numpy

if x <= 0:
return 2*n+1 + approx_grade(1, n+1)
else:
return 2**(10-x) * (9 - x) + (2**(11-x) -2)*(2*n + 1)

def sim(num_players):
xs = *num_players
ns = *num_players
while True:
i = numpy.argmin(grades)  # STRATEGY HERE
ns[i] += 1
if xs[i] == 0:
xs[i] = 1
elif random.random() <= 0.5:
xs[i] += 1
if xs[i] == 10:
return sum([n*n for n in ns])
else:
xs[i] -= 1

T = 100000
scores = [sim(4) for t in range(T)]
print(numpy.mean(scores))


This is a special case of the Gittins index theorem, and there is a very nice exposition of essentially this game in the 2003 paper On Playing Golf with Two Balls

Notice we can rephrase the cost as: when we move an agent that has taken n steps so far, we pay 2n+1 immediately, and the total cost at the end is the sum of all payments made along the way. Note 2n+1 = (n+1)^2 - n^2.

The key point is that for any given agent i and state of the agent (i.e. location x[i] and number of steps so far n[i]), there is an index or "grade" we can give the agent that perfectly captures how desirable it is to choose that agent next. The optimal strategy is to pick the agent with the best grade. One agent's grade does not depend on the others' states!

The grade is the following number. Imagine we had a fifth agent with the ability to magically teleport to the end for a cost g. And imagine we are playing the game with only agent i and the magic agent. If g is very small, we will just choose the magic agent. But if g is large, we will take a risk and play agent i, see what happens, and then decide whether to use the magic agent or play agent i again.

There is some magic number g(x[i],n[i]) where we are indifferent between the magic agent and i. That number is the grade.

Now, in theory the grade is polynomial-time computable, but it's pretty complex. So I used the following approximation: suppose the optimal strategy in the magic game is to keep trying agent i as long as they move right, but as soon as they move left, pay g to stop. Then the index for step 9 would satisfy

$$g(9,n) = 2n + 1 + \frac{1}{2}g(9,n) .$$

The left side is if we pay g(9,n) for the magic agent to halt immediately, and the right side is if we pay 2n+1 to advance agent i: with probability half we terminate and pay nothing, and with probability half we go left, at which point we pay g(9,n) to stop.

This solves to g(9,n) = 4n+2.

Similarly, at location 8, we should be indifferent between paying g(8,n) to stop, or paying 2n+1 to get a half chance to go left -- when we pay g(8,n) to stop -- and a half chance to go right -- when we repeat the game at step 9.

$$g(8,n) = 2n+1 + \frac{1}{2} g(8,n) + \frac{1}{2}\left(2(n+1)+1 + \frac{1}{2}g(8,n)\right) .$$

If we repeat this logic, we get an approximation to the grade which is

$$g(x,n) = 2^{10-x}(9-x) + (2^{11-x}-2)(2n+1) .$$

After dealing with the special case of x=0, I get this code. Apparently it's still pretty suboptimal as another answer gets ~4927.

Note: I used multiprocessing to run 8 million total reps in 36min on a laptop, and got an average of just under 4973. I also reimplemented in C++ with probably hundreds of millions of reps and it seemed to converge to about 4974.5. But the quality of the RNG and precision arithmetic matters! C's built in rand() gave basically nonsense.

• This is great. There was a paper on how to compute the gittins index in practice at arxiv.org/abs/1909.05075 Aug 10 at 16:11
• Does the gittins index optimality theorem still apply when the cost is not discounted? I thought there was usually a term like gamma^i in there with gamma < 1. Aug 10 at 16:14
• @gaffe Check out the paper, it's a really beautiful read. Sounds like finite expected termination time is sufficient. It's definitely not the usual Gittins, but close enough that they felt it should be called a variation!
– usul
Aug 10 at 16:17
• Welcome to Code Golf! Nice answer! Aug 10 at 16:20
• There must be! I was thinking about solving the magic game either analytically or with simulations, to find the right values of g. First when agent i is at location 9, then 8, then....
– usul
Aug 10 at 16:36