Category Archives: mathematica

Mental defectives, part ii

I did on a previous occasion write about the seemingly powerful relationship between the number of radios in the U.K, and the number of registered mental defectives in that country. In Google books, I recently found a 1958 book titled Principles of Statistical Techniques, which contains examples based on that remarkable dataset, the number of mental defectives in England and Wales.

The passage in question uses the “mental defectives” data from 1935 through 1946 to illustrate the ability to mislead through the use of graphs and charts. The graphs below are courtesy of Mathematica but closely replicate those in the original text.

In[]:= mentalDefectives = {{1935, 86086}, {1936, 88060}, {1938,
92299}, {1939, 99144}, {1940, 101364}, {1941, 100876}, {1942,
98125}, {1943, 98434}, {1944, 99608}, {1945, 102225}, {1946,

In[]:= ListPlot[mentalDefectives, PlotRange -> {86000, 102500},
AxesOrigin -> {1934, 86000},
AxesLabel -> {"Year", "Number of\nmental defectives"},
AspectRatio -> 2, PlotMarkers -> {Automatic, Small}]

tall skinny graph

In[]:= ListPlot[mentalDefectives, PlotRange -> {0, 120000},
AxesOrigin -> {1935, 0},
AxesLabel -> {"Year", "Number of\nmental defectives"},
PlotMarkers -> {Automatic, Small}]

wider graph

As explained in the original text, “The impression obtained from the first graph is that there was a staggering rise in the number of mental defectives in little more than a decade, whereas the second graph gives the impression of a much slower and more gradual increase. These rather different impressions are obtained, of course, by tampering with the horizontal and vertical scales used. . . . it is essential to recognize the importance of supplying graphs, as well as tables, with full and clear labelling and if possible, the source of the information.”

Frankly, their efforts seem pretty basic. We could further steepen the first graph through population adjustment, adjustment for number of doctors, and other techniques.

In googling for more on this remarkable dataset, I found letters to the editor of The British Medical Journal in which one R.A. Gibbons of London, S.W., argues that forced sterilization of mental defectives is not sufficient to protect “the race,” and that it is necessary to sterilize also those people who “from pronounced family history” are certain to “produce idiots.” Dr. Gibbons apparent delivered a paper on this topic to the Obstetrical Section of the Royal Society of Medicine and had in published in the Journal on March 18, 1922. In this letter, he refers to “over 12,000 registered mental defectives for whom [British citizens] are taxed.”

The reference to taxation feels terribly modern (“certainly their plight tugs at the heartstrings, but you just can’t care for everybody and it wouldn’t be fiscally prudent to allow them to run amok.”)

Germany is most famous for eugenics programs, but similar sterilizations of mental defectives were taking place in the United States in the early and mid 20th century. I’ll leave that for another post.

clock hands, part ii

I posted a request to mathgroup asking for guidance on why Solve[] and its brethren failed on the puzzle of the opposing clock hands. Two people have already responded with working versions from which we can draw the following conclusions,

1. We have to unwind the linedupedness[], unitTimeToHourAngle[], and unitTimeToMinuteAngle[] functions, and

2. We need to restrict the domain of the solution to the range from 0 to 1.

When we do this, we go from

In[]:= Reduce[unstraightness[unitTime] == 0 && 0 < unitTime < 1]
In[]:= Solve[Abs[Abs[4320 Mod[unitTime, 1/12] - unitTime*360] - 180] == 0 &&
0 < unitTime < 1, unitTime]

which does indeed produce the right answer. There are other ways to get there too; unwinding and restricting the solution to Reals or Rationals will work, though we get some redundant solutions that way.

Wolfram tech support had an even nicer solution, pointing out that my wrapper function (unstraightness[]) could take any kind of input, while the internal functions (unitTimeToHourAngle[] and unitTimeToMinuteAngle[]) are restricted to real numbers. This means that when these functions are introduced into the Reduce function, unitTime is not defined as a number, so the function is never evaluated. The simplest solution is just to change the definition of the internal functions to

In[]:= minuteAngle[unitTime_] := 4320 Mod[unitTime, 1/12]
In[]:= hourAngle[unitTime_] := unitTime*360 // N


In[]:=Reduce[Abs[minuteAngle[unitTime] - hourAngle[unitTime]] == 0 &&
0 < unitTime < 1]

works like a charm.

Opposing clock hands

I was recently sitting in a poorly produced play called That Hopey Changey Thing when, looking wistfully at my watch, I noticed that the hour hand and the minute hand were perfectly lined up, pointing in opposite directions. I stayed awake for the rest of the show in part thanks to the thought puzzle of how many times a day this happens, and how best to determine exactly when.

Here’s an answer using Mathematica.

First, I found it useful to think of time in terms of a unit clock cycle — the hands on our theoretical clock run from 0 to 1, rather than to 12. The following Manipulate[] command draws a clock, allowing the user to set the unit time, and translating it into both conventional time and the angle of the minute and hour hands.

In[]:= Manipulate[ Graphics[{{Opacity[0], Disk[{0, 0}, .9]}, {Opacity[.8], Thickness[.02], Line[{{0, 0}, .65 {Sin[2 Pi m], Cos[2 Pi m]}}]}, Line[{{0, 0}, .85 {Sin[12*2 Pi m], Cos[12*2 Pi m]}}]}, PlotLabel -> ToString[Floor[Mod[11 + m*12, 12] + 1]] <> " hours, " <> ToString[360.*m] <> "\[Degree]\n" <> ToString[Floor[Mod[m, 1/12]*720]] <> " minutes, " <> ToString[4320 Mod[m, 1/12]] <> "\[Degree]"], {{m, 0, "time of day"}, 0, 1}]

The code when running looks like
manipulate[] clock
This was useful when debugging the formulae for the conversions; any error was immediately apparent.

Breaking out the conversions explicitly, we have

In[]:= unitTimeToListTime[unitTime_Real] := Module[{hour, minute, seconds}, hour = Floor[Mod[11 + unitTime*12, 12] + 1]; minute = Floor[Mod[unitTime, 1/12]*720]; seconds = Mod[unitTime, 1/720]*60*720; {hour, minute, seconds} ]
In[]:= unitTimeToMinuteAngle[unitTime_Real] := 4320 Mod[unitTime, 1/12]


In[]:= unitTimeToHourAngle[unitTime_Real] := unitTime*360 // N

These in hand, I naively thought the problem was as good as solved. I defined my problem in functional form,

In[]:= unstraightness[unitTime_] := Abs[Abs[unitTimeToMinuteAngle[unitTime] - unitTimeToHourAngle[unitTime]] - 180]

checked it with

Plot[unstraightness[unitTime], {unitTime, 0, .999999}]

and we see as expected a periodic function. There are 11 solutions when “unstraightness” == 0 per day, which makes perfectly good sense when you think about it.

I thought I could determine these with NSolve[unstraightness[unitTime] == 0, unitTime]

but sadly, this yields nothing but error messages.
NSolve::nsmet: This system cannot be solved with the methods available to NSolve. >>

What to do? I decided on brute force. Slice the day into 1,000,000 pieces and identify the 11 times which produce the most opposed clock hands.

In[]:= handangle = Table[{unitTime, Abs[unitTimeToMinuteAngle[unitTime] - unitTimeToHourAngle[unitTime]]}, {unitTime, 0, .999999, .000001}];

In[]:= closeFits = Select[handangle, Abs[#[[2]] - 180] < .002 &];

In[]:= TableForm[
Map[{#, ToString[unitTimeToListTime[#]], unitTimeToHourAngle[#],
unitTimeToMinuteAngle[#]} &, Transpose[closeFits][[1]]],
TableHeadings -> {None, {"unit time", "clock time", "hour angle",
"minute angle"}}]


As looked to be the case from the graph, these solutions are clearly spaced exactly 1/11 of the clock apart from each other. This means that we can get a more perfect solution more quickly with

In[]:= Table[ToString[unitTimeToListTime[i // N][[1]]] <> ":" <> ToString[unitTimeToListTime[i // N][[2]]] <> ":" <> ToString[unitTimeToListTime[i // N][[3]]], {i, 1/22, .9999, 1/11}] // TableForm

which returns


Clearly, my bored glance at my watch had been at 8:10:54 pm.

One of the nice things about this solution method is that it easily generalizes to solve other “clock face” problems. For example, in a second we can determine that the hands of a clock are perfectly lined up, facing the same direction, at the following times:

recursion in Mathematica

I have been playing with creating pseudotexts via markov chains recently, and may post some code soon. Creating random text via markov chains has two steps — first you analyze the source material and make a tree that describes the likelihood of different words following each other, and second you create a random walk through that tree, respecting the relative probability of the different paths.

The first step is easily done via many different methods, but if you want a tree of arbitrary depth, it becomes most elegant to use recursion. This is no harder in Mathematica than any other language, and for the sake of any novices googling around for examples of using Mathematica recursively, note the following example —

If you want to compute factorials, in practice you would probably use Mathematica’s postfix operator, !

In[1]:= 4!
Out[1]= 24

You could program your own version about a thousand different ways. I find this to be natural:

In[2]:= functionalFactorial[x_] := Apply[Times, Range[x]]
In[3]:= functionalFactorial[4]
Out[3]= 24

And if you wanted a recursive version, you could do

In[4]:= recursiveFactorial[x_] := If[x == 1, 1, x*recursiveFactorial[x - 1]]
In[5]:= recursiveFactorial[4]
Out[5]= 24

As in any recursive code, you need a test in the inner loop that defines at least one stopping point (in this case when we're multiplying by 1). Without that branch, all recursions would continue indefinitely.

Either of the user-defined functions above could be compiled for additional speed.

Jewish New Year Distribution

I was talking to somebody the other day who was pleased because his birthday (September 9) landed on the Jewish New Year. He made the dubious claim that this happened something like 17% of the time.

For a quick check, I used the JewishNewYear[] function built in to Mathematica’s Calendar package. It is good from 1900 through 2100.


In[]:= Length[Select[
Map[JewishNewYear[#] &,
Range[1900, 2100]], #[[2]] == 9 && #[[3]] == 9 &]]

Out[]= 7

7/Length[Range[1900, 2100]] == .0348, so over this 201 year time span, his birthday would have landed on Rosh Hashana a bit under 3.5% of the time.

Curious about how the holiday distributed, I generated a table of frequencies (all displayed as though they were 2010 dates so we can use the built-in DateListPlot[] function).

In[]:= dateFreq =
Map[{#[[1]], Length[#]} &,
Map[{2010, JewishNewYear[#][[2]], JewishNewYear[#][[3]]} &,
Range[1900, 2100]], AbsoluteTime[#1] <= AbsoluteTime[#2] &]]]; In[]:= DateListPlot[dateFreq, Filling -> Axis,
FillingStyle -> {Automatic, {Black, AbsoluteThickness[10]}},
GridLines -> None, Frame -> False, Axes -> True, AxesStyle -> Black,
AxesOrigin -> {Automatic, 0}, PlotRange -> All,
PlotLabel -> "Rosh Hashana Distribution"]

Note the nested use of Length[], Gather[], Sort[]. That is often a powerful combination when using Mathematica to analyze the frequencies of events or patterns in data.

In this case we do not have a uniform distribution — it clearly fades at the edges, but it’s not normally distributed either. Figuring out whether or not there is a meaningful pattern here is beyond my statistical skills.