5
$\begingroup$

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

Improve this question
$\endgroup$
14
  • 2
    $\begingroup$ Quick sanity check — you say 'family of reals', but if I'm reading this right the elements of $F$ generated by this procedure would all have to be dyadic rationals; is that correct? Also, do you mean that $F$ should be a set? (i.e., if we choose $\langle x_1, y_1\rangle = \langle 0,1\rangle$ and then $\langle x_2, y_2\rangle=\langle 0,1\rangle$ as well, then after two rounds $F$ will be just $\{0, \frac12, 1\}$?) $\endgroup$
    – Steven Stadnicki
    Commented Mar 26, 2021 at 23:10
  • 1
    $\begingroup$ That makes sense. If it's a multiset, then the 'without replacement' condition doesn't ensure any more that $x_r\neq y_r$ unless I'm missing something — in the example, you could choose the two 'different' instances of $\frac12$ in round three. $\endgroup$
    – Steven Stadnicki
    Commented Mar 26, 2021 at 23:49
  • 1
    $\begingroup$ The phrase "without replacement" usually means that the items are removed from the set and not put back. But you are just choosing two distinct elements (perhaps with the same value, but distinct in the multiset) and adding a new element without removing those two; is that right? $\endgroup$
    – Brendan McKay
    Commented Mar 27, 2021 at 6:27
  • 2
    $\begingroup$ Have you run a simulation? I'm not sure that it will always cluster around $\frac12$. Actually I'm pretty sure it won't. $\endgroup$
    – Brendan McKay
    Commented Mar 27, 2021 at 12:22
  • 1
    $\begingroup$ Interesting! I tried 2-dimensional version: start with the set of three points in the plane, then at each step add the barycenter of a random triple of points from the set. Experiments show that the standard deviation of barycenters of the final sets does hardly depend on the number of iterations. $\endgroup$
    – მამუკა ჯიბლაძე
    Commented Mar 28, 2021 at 11:49

3 Answers 3

Reset to default
6
$\begingroup$

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}$.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Thank you very much MartinHairer! Do you think that the last plot of Bill Bradley would show a concave curve if one ran a significantly larger number of experiments? Another question: do you think that the process significantly changes using a set instead of a multiset (i.e., if each number of $M$ had multiplicity $1$)? $\endgroup$
    – Penelope Benenati
    Commented Mar 29, 2021 at 16:53
  • 1
    $\begingroup$ I wonder too about allowing $x_r=y_r$ (without accounting for multiplicities) One might think that after a while it would hardly matter (one might be wrong). Just that then if $M$ was described by $f(x)=\sum p_i t_i$ where $p_i$ is the probability that a randomly chosen element of $M$ had value $t_i,$ Then the probabilities for various values of the next element are given by $f(\sqrt{x})^2.$ $\endgroup$
    – Aaron Meyerowitz
    Commented Mar 30, 2021 at 1:02
  • 1
    $\begingroup$ @AaronMeyerowitz It wouldn't matter for the analysis I gave here, but it might influence the distribution of the limiting median (last plot in Bill's simulations), since that is expected to be dominated by the small N behaviour. $\endgroup$
    – Martin Hairer
    Commented Mar 30, 2021 at 9:25
  • 1
    $\begingroup$ There is then a chance of $\frac19$ that at $n=4$ $M=\{0,0,0,1\}$ and an equal chance of $M=\{0,1,1,1\}$. Either would have a significant skewing effect for a while. $\endgroup$
    – Aaron Meyerowitz
    Commented Mar 30, 2021 at 19:16
6
$\begingroup$

This isn't an answer; I'm just sharing plots from a few numerical experiments. Each time we repeat the above process for $N$ steps, we generate a (potentially) different multiset (i.e., a different sample of the process). We will refer to the $i$-th such multiset as $M_N^i$.

First observation: $M_N^i$ appears smooth and unimodal. For $N=1,000,000$: enter image description here Second observation: $M_N^i$ does not look normally distributed. enter image description here Third observation (per Martin Hairer's estimates): $M_N^i$ is much better approximated by a Cauchy distribution, but it still doesn't look quite right. In the figure below, we plot $M_{1,000,000}^1$ and two million-point samples of Cauchy distributions. We sampled the Cauchy twice to show the variability; in particular, they're visually nearly indistinguishable, but the distribution of $M_{1,000,000}$ is noticeably different. (Note that fitting a Cauchy is non-obvious since it has no mean value; we matched the median and the interquartile range.) enter image description here In the comments, Martin Hairer points out that the interquartile range might be too wide for a good fit. Following his suggestion, if we scale so we fit the scale between the 40th and 60 percentiles, the probability density (i.e. height) matches much better, but the discrepancy still seems marked away from the mode. Since a Cauchy only has one degree of scale freedom, it doesn't seem that we can do better, unfortunately. Here's the figure: enter image description here Fourth observation: For fixed $i$, as $N\rightarrow\infty$, it looks like $M_N^i$ converges to its mean. Plotting the distributions of $M_{1000}^1,M_{10,000}^1,M_{100,000}^1,M_{1,000,000}^1$: enter image description here Fifth observation: On the other hand, different runs appear to have substantially different means. Comparing $M_{1,000,000}^1$ to $M_{1,000,000}^2$: enter image description here Final observation: We can examine how much the means vary. Let $\mu^i$ be the sample mean of $M_{1,000}^i$. Taking $i=1,...,400,000$, we observe the following (remarkably wide) distribution of $\mu^i$: enter image description here If we graph $M_{10,000}^i$, the graph is virtually identical.

This behavior has a sort of "Polya's Urn" feel to it.

Improve this answer
$\endgroup$
14
  • 1
    $\begingroup$ I had identified the histogram for one run as probably normal, but you convinced me otherwise. Rather than the fourth central moment being 3 times the square of the variance, I'm getting 20-40 times the square of the variance. $\endgroup$
    – Brendan McKay
    Commented Mar 28, 2021 at 6:19
  • 1
    $\begingroup$ @BillBradley Would you be able to compare the simulation in your second graph to a Cauchy distribution? (See my answer below for a justification.) $\endgroup$
    – Martin Hairer
    Commented Mar 29, 2021 at 7:37
  • 1
    $\begingroup$ @MartinHairer Oh, I actually tried playing with the scale and stretched it by 1.15x; this matches the probability density of the height, but the width is still off. (Again, not terribly so, but far enough that it looks worse than just statistical noise). In any event, your suggestion is a good one. My day job's going to keep me from posting anything now, but I'll try to fit a Cauchy to a chunk of the middle region and see if I can post it later. $\endgroup$
    – Bill Bradley
    Commented Mar 30, 2021 at 13:34
  • 1
    $\begingroup$ @MartinHairer OK, plot added! I wish the figure were more encouraging :( $\endgroup$
    – Bill Bradley
    Commented Mar 31, 2021 at 4:12
  • 1
    $\begingroup$ I wonder if the peaks of the bimodal plot “are” $\frac38$ and $\frac58$ . If you still have the data then a plot with bins centered at appropriate dyadic rationals might be interesting. $\endgroup$
    – Aaron Meyerowitz
    Commented Apr 1, 2021 at 5:59
3
$\begingroup$

This is just a comment, but I wanted to include a graph. There are are $2^r-1$ possibilities for $z_r$. $z_r=\frac{p}{2^r}$ with$1 \le p \le 2^r-1.$ An interesting observation is that the denominator (in lowest terms) is unlikely to be close to $2^r$. If $z_r$ happens to have denominator $2^r$ then it is the only one with that denominator in (the current version of) $M$ and , aside from $\frac01$ and $\frac11$, each of $2,2^2,\cdots, 2^r$ is used as a denominator once. Then the denominator of $z_{r+1}$ is $2^s$ with probability $\frac{2s}{r^2+r}.$ so on average $s \approx \frac{2r}3.$

Here is the probability distribution for $z_6.$ Just under half the time it is $\frac12,\frac38,\frac58,\frac14,$ or $\frac34$ with frequencies roughly $12\%,10\%,10\%,8.5\%,8.5\%. $ enter image description here

That may not have any effect on the average behavior for very large $r$.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thank you for the interesting plot Aaron! $\endgroup$
    – Penelope Benenati
    Commented Mar 29, 2021 at 16:56

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .