Exploring the Mean-Median Map

Math 242 Modern Computational Math

A key observation is that the important properties of this map on a starting set of three real numbers can be reduced to a single parameter between 0 and 1.

That is, we we can just start with values (0, a, 1) and study the effects of the map on the parameter a.

Also, given the numerical complexity iterating the mean-median map, it is best to work with rational numbers to avoid round-off error. You can force Mathematica to convert decimal values to rational form using its built-in Rationalize function.

Define the important functions

Just as last time, we need a function to return the next value in the sequence:

nextVal[vals_]:=Module[{next,n=Length[vals]},
next = (n+1)Median[vals]-Total[vals];
next  (* return next *)
]

Also as last time, we need a function that returns the sequence, up to the nth term.

computeSeq[vals_,length_]:=Module[{next,newvals=vals},
Do[
next=nextVal[newvals];
AppendTo[newvals,next],
length   (* number of times to repeat = length *)
];
newvals  (* return the sequence *)
]

It would be useful to have a function that computes the sequence until a value repeats. Here it is:

limitSeq[vals_]:=Module[ 
(* initialize some values as local variables *)
{prevVal=vals[[-2]],
currVal=vals[[-1]],
currSeq=vals},

While[prevVal≠currVal,
(* promote prevVal and compute nextVal *)
prevVal=currVal;
currVal=nextVal[currSeq];
(* append currVal to currSeq *)
AppendTo[currSeq,currVal];
];

(* return the sequence, which ends with two repeated values *)
currSeq
]

Try out these functions. Here is an interesting set of starting values:

vals0={1,3,17}

Compute the first 80 terms:

computeSeq[vals0,80]

Now use limitSeq:

limitSeq[vals0]

Looks good. Note that the number of steps until the sequence repeats, and the repeated value itself, may be obtained as follows:

seq0=limitSeq[vals0];
Length[seq0]
Last[seq0]

Looks like it works.

Note: If we want to write efficient code, we should think about the intermediate steps of totaling and taking the median. Both of these are computationally expensive. If we kept a running total, we wouldn’t have to add up the entire sequence every time we want to compute the next value. Also, if we kept the sequence values sorted in order, we could find the median much more quickly. We won’t worry about these enhancements right now, but you could think about them for improving your code later.

Reduction to {0,a,1}

We can simplify using the results from the paper by Chamberland and Martelli. They point out that it suffices to consider starting sets of the form {0,a,1} where (1/2)≤ a ≤ (2/3).

Write two small functions to make this easy:

limSeq0[a_]:=limitSeq[{0,a,1}]
computeSeq0[a_,n_]:=computeSeq[{0,a,1},n]

Quick test:

a0=51/100
seq0=limSeq0[a0]
computeSeq0[a0,100]

Looks good!

How many steps to hit the repeated value? What is the repeated value?

length0=Length[seq0]
lim0=Last[seq0]

To match things up with the paper, define the function m(a) as the limiting value from the starting sequence (0,a,1). We can find this easily with our function limSeq0:

m[a_]:=Last[limSeq0[a]]

m[a0]

I’ve seen this value before.....

From the paper:

"We point out that every numerical investigation contains the presence of two equal values (Subscript[x, j]) = (Subscript[x, i]) before the sequence becomes stationary. We have never found an example of convergence without this feature. We cannot explain the occurrence of the two equal values, which may be separated by a large number of iterations."

The statement could a bit stronger—the two equal values appear to the be limiting value itself!

Let's check that here. Search the sequence seq0 for occurrences of lim0 and also their positions in the sequence. (Note that we use pure functions in the implentation below.)

Select[seq0,(#==lim0&)]
Position[seq0,_?(#==lim0&)]//Flatten

Yes, the value occurs twice before!

Thus we see the limiting value occurs at the locations 27 and 52 before appearing at the end. Does this always happen? Or at least appear to always happen?

We would like to have a way of checking a large number of starting values (i.e., values of the parameter a). For each starting value, we determine if the limiting value occurs before the sequence becomes constant.

This is an important technique: taking a previous computation and "encapsulating" it as a function.

prevLocs[a_]:=Module[{lim,len,seq,pos},
seq=limSeq0[a];
len=Length[seq];
lim=Last[seq];
pos=Position[seq,_?(#==lim&)]//Flatten;
pos
]

prevLocs[a0]

Now change the starting value:

a0=Rationalize[.1252]
locs=prevLocs[a0]

We should check to make sure we found the correct values. To do this, compute the sequence and select the locations returned the previous command.

limSeq0[a0][[locs]]

Looks good.

Now apply this analysis to a large set of values. Let’s choose values between (1/2) and (2/3), with a step size of (5/1000).

aVals=Table[a,{a,1/2,2/3,5/1000}]

Here we go:

pLocs=prevLocs/@aVals

See how long each of these is.

lens=Length/@pLocs

Wait—is there a 2 in this list? What does that mean?

Get the position of the 2:

pos2=Position[lens,2]//Flatten
pos2=pos2[[1]]

Did we find a counterexample to the authors’ claim?

Let’s dig deeper. First get the starting value that produced this situation.

a2=aVals[[pos2]]

Now look at the resulting limit sequence:

limSeq0[a2]

It looks like the repeated value didn’t occur previously in the sequence! However, let’s compute more terms of the sequence:

computeSeq0[a2,40]

Now we see the problem! Our limitSeq function is stopping too soon! It’s possible that the sequence might repeat in consecutive steps, but the sequence is not (yet) constant. However, if the medians repeat in consecutive steps, then the sequence is constant. Why is this?

You may either try to fix limitSeq now, or work through the rest of this notebook and then come back to limitSeq.

Plotting Figure 1 from the paper

Figure 1 in the paper shows the effect of taking the kth median of sequences at various starting values. Let’s recreate this plot.

To begin, create  table of starting values. We can use the values we had earlier.

aVals=Table[a,{a,1/2,2/3,5/1000}]

For each of these values, compute a long sequence of values. Here we use the 20. Note the use of the pure function. The function seqVals0 has two arguments, but we only want to use the first argument and fix the second at 20.

lotsOfSeqs=(computeSeq0[#,20]&)/@aVals

We need a function which will compute the median of the first k terms.

medk[seq_,k_]:=Median[seq[[1;;k]]]

Try it out:

lotsOfSeqs[[1]]
medk[lotsOfSeqs[[1]],3]

Here’s the fun part. We can apply our new medk function, once again fixing the second argument.

As well, we make a plot of the resulting values.

vals=medk[#,10]&/@lotsOfSeqs
ListLinePlot[vals,DataRange->{First[aVals],Last[aVals]},AxesLabel->{"a","10th median"}]

Now do this for each of over the first k=4,...,20 terms of the sequence.

Do[meds[k]=medk[#,k]&/@lotsOfSeqs,{k,4,20}]

Package all these together in a single plot. This is similar to Figure 1 in the paper.

ListLinePlot[Table[meds[k],{k,4,20}],DataRange->{First[aVals],Last[aVals]},AxesLabel->{"a","median value"}]

What' s next?

If you haven’t done so already, rewrite the limitSeq function so that it stops when consecutive medians are equal.

Then, start work on the Mean-Median Project, which you can find on the course web site.

As you investigate the mean-median map, ask yourself: What do you observe? What do you wonder?