We are given a *multiset* $M$ of real numbers which initially is equal to $\{0,1\}$. In a sequential fashion, at each round $r\in\mathbb{N}$


* two *distinct instances* $x\_r$ and $y\_r$ of $M$'s numbers are selected uniformly at random from $M$ (which implies that they cannot be the *same instance* of any number contained in $M$, viz., $x\_r$ is selected and temporarily removed from $M$, thereafter $y\_r$ is selected from $M\setminus\{x\_r\}$ without removing it, and finally $x\_r$ is added back to $M$), and
* $z\_r=\frac{x\_r+y\_r}{2}$ is added to $M$.




---


**Question:** What is the probability $p\_{r,{\epsilon}}$ that we have $\left|z\_r-\frac{1}{2}\right|\le\epsilon$ for a given $\epsilon\in\left(0,\frac{1}{2}\right)$?




---


***Edit:*** *I show below the experimental results for different simulations of the random process run with an increasing total number of rounds, and for a single simulation of the random process keeping track of the evolution over time of the total average and the last element added - For the sake of convenience, by writing "rounds", here I am counting the initial insertion of both $0$ and $1$ (simultaneously) as the very first round:*


New simul. with 2^1 rounds - Avg: 0.5 | Last added: 0.5


New simul. with 2^2 rounds - Avg: 0.5 | Last added: 0.5


New simul. with 2^3 rounds - Avg: 0.472222 | Last added: 0.5


New simul. with 2^4 rounds - Avg: 0.430147 | Last added: 0.375


New simul. with 2^5 rounds - Avg: 0.413826 | Last added: 0.40625


New simul. with 2^6 rounds - Avg: 0.40012 | Last added: 0.40625


New simul. with 2^7 rounds - Avg: 0.38313 | Last added: 0.46875


New simul. with 2^8 rounds - Avg: 0.377516 | Last added: 0.378906


New simul. with 2^9 rounds - Avg: 0.366866 | Last added: 0.378906


New simul. with 2^10 rounds - Avg: 0.362607 | Last added: 0.342743


New simul. with 2^11 rounds - Avg: 0.360595 | Last added: 0.358353


New simul. with 2^12 rounds - Avg: 0.359471 | Last added: 0.343569


New simul. with 2^13 rounds - Avg: 0.364962 | Last added: 0.336161


New simul. with 2^14 rounds - Avg: 0.500135 | Last added: 0.497771


New simul. with 2^15 rounds - Avg: 0.49995 | Last added: 0.488623


New simul. with 2^16 rounds - Avg: 0.602851 | Last added: 0.590848


New simul. with 2^17 rounds - Avg: 0.376087 | Last added: 0.372888


New simul. with 2^18 rounds - Avg: 0.655107 | Last added: 0.62898


New simul. with 2^19 rounds - Avg: 0.182425 | Last added: 0.201142


New simul. with 2^20 rounds - Avg: 0.709139 | Last added: 0.713385


New simul. with 2^21 rounds - Avg: 0.219937 | Last added: 0.220374


New simul. with 2^22 rounds - Avg: 0.112707 | Last added: 0.112427




---




---


Same simul. r=2^1 - Avg: 0.5 | Last added: 0.5


Same simul. r=2^2 - Avg: 0.5 | Last added: 0.75


Same simul. r=2^3 - Avg: 0.545139 | Last added: 0.46875


Same simul. r=2^4 - Avg: 0.625 | Last added: 0.28125


Same simul. r=2^5 - Avg: 0.60393 | Last added: 0.59375


Same simul. r=2^6 - Avg: 0.568329 | Last added: 0.71875


Same simul. r=2^7 - Avg: 0.57769 | Last added: 0.697266


Same simul. r=2^8 - Avg: 0.573474 | Last added: 0.631714


Same simul. r=2^9 - Avg: 0.575036 | Last added: 0.576538


Same simul. r=2^10 - Avg: 0.577153 | Last added: 0.47583


Same simul. r=2^11 - Avg: 0.578355 | Last added: 0.617221


Same simul. r=2^12 - Avg: 0.576684 | Last added: 0.57461


Same simul. r=2^13 - Avg: 0.576757 | Last added: 0.581285


Same simul. r=2^14 - Avg: 0.577305 | Last added: 0.546254


Same simul. r=2^15 - Avg: 0.577683 | Last added: 0.592735


Same simul. r=2^16 - Avg: 0.577662 | Last added: 0.56319


Same simul. r=2^17 - Avg: 0.577692 | Last added: 0.576607


Same simul. r=2^18 - Avg: 0.577675 | Last added: 0.571428


Same simul. r=2^19 - Avg: 0.577657 | Last added: 0.572818


Same simul. r=2^20 - Avg: 0.577655 | Last added: 0.579482


Same simul. r=2^21 - Avg: 0.577652 | Last added: 0.575974


Same simul. r=2^22 - Avg: 0.577654 | Last added: 0.5777


Same simul. r=2^23 - Avg: 0.577659 | Last added: 0.585123


Same simul. r=2^24 - Avg: 0.577657 | Last added: 0.571693


Same simul. r=2^25 - Avg: 0.577659 | Last added: 0.579782


Same simul. r=2^26 - Avg: 0.577659 | Last added: 0.574194


Same simul. r=2^27 - Avg: 0.577659 | Last added: 0.579098



Convergence speed of a random dyadic rational generator
This is a cute problem... Let's normalise your multiset $M$ to a probability measure $\mu$ by setting $\mu = {1\over |M|} \sum\_{x \in M} \delta\_x$ (repeated elements are repeated in the sum). Write also $h \colon x \mapsto x/2$. If we then set $t = \log N$ for $N$ the number of steps, a very good model for the evolution of $\mu$ by the time $N$ is very large is given by
$$
\partial\_t \mu\_t = -\mu\_t + h^\*(\mu\_t \* \mu\_t) \tag{$\star$}
$$
where $\*$ denotes convolution. (This is the mean field behaviour, which is of course a bad estimate for small $N$. In particular, the distribution for the mean of $\mu$, which remains fixed under this evolution, would be dominated by the small $N$ behaviour and might depend on details like whether the two instances are drawn from $M$ with or without replacement.)


Now the first thing we notice is that the fixed points for $(\star)$ are precisely the Cauchy distributions, which strongly suggests that this is what Bill sees in the second simulation. On the other hand, we know of course that $\mu$ remains compactly supported, so that if it converges to a fixed point, it has to be a delta. This strongly suggests that for large $t$, $\mu\_t$ is well approximated by a Cauchy distribution $\mu\_t^\star$ at a fixed (but random) location $x\_0$ and with scale parameter $\gamma\_t \to 0$.


An interesting question then is how fast does $\gamma\_t$ converge to $0$? For this, we can make the ansatz $\mu\_t = \mu\_t^\star + \nu\_t$ for some remainder $\nu\_t$ such that $\langle\nu\_t, \mu\_t^\star\rangle = 0$ ($L^2$ scalar product of densities). This last condition allows in principle to uniquely determine $\gamma\_t$, at least in the regime where $\nu\_t$ is small. We can then derive a coupled system of ODEs for $\nu\_t$ and $\gamma\_t$, but this looks very messy.


If we only keep the dominant terms however, the equation for $\gamma$ looks like
$$
\dot \gamma\_t = - C \gamma\_t^2 \langle h^\*(\nu\_t \* \mu\_t^\star), \mu\_t^\star\rangle
$$
for some constant $C$ which can in principle be computed explicitly. At this point, we can make a leap of faith and assume that the dominant effect of $\nu\_t$ is to 'chop the tails' off $\mu\_t^\star$ to force it to be a distribution on $[0,1]$. Unless I made a mistake in my back-of-the-envelope calculation, this finally leads to
$$
\dot \gamma\_t \approx - \tilde C \gamma\_t^4\;,
$$
which suggests that one should have $\gamma\_t \propto t^{-1/3}$ for large $t$.
In other words, after $N$ steps, $\mu$ should be quite well approximated by a Cauchy distribution centred at some random location (whose law is an unknown distribution plotted in Bill's fourth graph which is unlikely to have a nice closed-form expression), with scale about $(\log N)^{-1/3}$, and tails cut off to fit into $[0,1]$.


In particular, one would expect the $k$th centred moment of this distribution to be of order $(\log N)^{-1/3}$ for every $k$, so the ratio of the fourth cumulant to the square of the variance would be expected to grow like $(\log N)^{1/3}$.


