# EMT – Simulation of the Two Child Problem

The two child problem is a problem in probability theory with a solution that seems paradox on first sight. I wrote about that problem before. Let me repeat the explanation and do a simulation to convince everybody that the solution is really correct. You can find information on that problem on Wikipedia.

The first simple version starts with this story: You meet a man in a bar and he mentions his daughter. You ask him about his children and he answers that he has two kids. The question is: What is the probability of the other kid (the one he was not talking about) to be female?

The answer depends on some important details. But let us assume you have no information if the other kid is younger or older, or any other special information about the daughter and the man is just a random man from the population. These details will turn out to be important! Then the intuitive answer to our question 1/2 is wrong, or rather, it does not make sense. The correct answer is 1/3.

In my opinion, each probabilistic problem makes sense only if we can devise an experiment that would in theory or in practice simulate the problem. Our computation should predict the outcome of the experiment. Problems that cannot be simulated in any way are not interesting to me. This includes most problems in probability or measure theory which need the help of the axiom of choice.

But our problem can easily be simulated. And setting up the simulation gives great insights into the meaning of our question and the terms „probability“ and „randomly selected“. Let us do it in Euler Math Toolbox (EMT). What we do is a Monte-Carlo simulation. We need to make the following assumption: The man has two kids with random gender, one of which is female, and the probability for a kid to be either gender is 1/2 (no diverse genders in this posting). So we draw 1000000 pairs of kids with random gender. Then we count the proportion of two female kids among all pairs with at least one female kid.

>n = 1000000
1000000
>K = intrandom(2,n,2)
Real 2 x 1000000 matrix

2             1             1             2     ...
1             2             2             2     ...
>i = nonzeros(K[1]==1 || K[2]==1); ni = length(i); ni/n
0.749596
>sum(K[1,i]==1 && K[2,i]==1) / ni
0.334372115113


The syntax may seem cryptic, but it is intuitive if you understand the Matrix language. K contains a pair of kids in each column (n=1000000 columns). „K[1]==1“ returns a vector of 0/1 with 1 (true) on each position where the vector K[1] (the first row of K) is 1. I.e., a vector indicating the pairs where the first kid is female. „||“ is short for „or“, and nonzeros() returns the indices of the non-zero elements of a vector. Thus „ni“ is the number of pairs such that either kid is female. As expected, „ni/n“ is approximately 3/4. There are four cases, and one case (tow boys) is wrong.

In the final line, we count the numbers of pairs in the „i“-columns, where both kids are female, using sum() which sums up the ones and compare that to the total number of right cases „ni“. The answer is approximately 1/3.

This should not surprise us since there are three cases in the „i“-columns: (1,2), (2,1), (1,1). Only one of these cases is the correct one.

Why does the problem depend on the details? For this, we assume that it is Tuesday and the man in the bar has her birthday today. Surprisingly, this changes the problem completely! Even if we only know that the daughter is born on a Tuesday the problem changes drastically.

Let us start with the Tuesday problem. We simulate the same now by randomly selecting weekdays for the birthday of both kids.

>n = 1000000;
>K = intrandom(2,n,2);
>D = intrandom(2,n,7)
Real 2 x 1000000 matrix

7             4             6             4     ...
7             2             5             5     ...
>i = nonzeros((K[1]==1 && D[1]==2) || (K[2]==1 && D[2]==2)); ni = length(i); ni/n
0.137799
>sum(K[1,i]==1 && K[2,i]==1) / ni
0.481679838025


The code did not change very much. But the „i“-columns now contain only columns with one Tuesday girl (2 means Tuesday above). The probability for this much lower, of course. We have 4*49=196 cases in total (4 gender pairs, 7 days for each). Of these, only 13 contain a Tuesday girl. There are 6 with the first kid a Tuesday girl and the other a girl born on another day, and 6 with the younger in the same way, and one with two Tuesday kids, plus 14 cases with one Tuesday kid and a boy. This is 27/496~0.13775.

Out of these 27 cases, we have 13 cases with two girls, as computed above. We have 13/27~0.48148. The simulation works as good as it can. The accuracy of a  Monte-Carlo simulation is only about 1/sqrt(n)=1/1000.

If we know that the birthday is today we can just take the probability as 1/2. In this case, the pick of the other child is almost independent of the pick of the first child in contrast to the original problem. And that is the true reason for the confusion. If we had known that the girl is the elder one the pick of the other child is independent and thus female with probability 1/2. But we only know that one of the kids is female. That changes the situation completely and both picks depend on each other.

I recently saw a Youtube video which I do not link here because it only adds to the confusion. The video could not explain why it makes a difference between knowing that the girl has a birthday, and knowing the precise birthday. This is easy to explain if you think of an experiment as above. Drawing a man with two kids, one of them being a daughter, is a different experiment than drawing a man with two kids, one of them is a daughter born on a Tuesday.

# Yahtzee Waiting Times

I recently was asked about waiting times in the game of Yahtzee. If you do not know the game it suffices to say that you throw 5 dice, and one of the goals is to get 5 dice with the same number, a Yahtzee. I wrote about waiting times a long time ago. But let me repeat the main points.

The trick is to use a number of states S0, S1, …, Sn, and probabilities P(i,j) to get from state Si to state Sj. S0 is the starting state, and Sn is the state we want to reach, in our case the Yahtzee. For a first attempt, we use the number 6s we have on out 5 dice. Then we have 6 states, S0 to S5. With a bit of combinatorics, we can compute the probabilities P(i,j) as

$$p_{i,j} = p^{j-i} (1-p)^{n-i-(j-i)} \binom{n-i}{j-i}$$

If we compute that with Euler Math Toolbox we get the following matrix P.

>p=1/6;
>i=(0:5)'; j=i';
% This is the matrix of probabilities to get from i sixes to j sixes.
0.4018776 0.4018776 0.1607510 0.0321502 0.0032150 0.0001286
0.0000000 0.4822531 0.3858025 0.1157407 0.0154321 0.0007716
0.0000000 0.0000000 0.5787037 0.3472222 0.0694444 0.0046296
0.0000000 0.0000000 0.0000000 0.6944444 0.2777778 0.0277778
0.0000000 0.0000000 0.0000000 0.0000000 0.8333333 0.1666667
0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 1.0000000

This matrix allows us to simulate or compute results about probabilities for 6-Yahtzees.  E.g., we can start with no 6 in S0. After one dice throw, the first row of P yields the distribution of the expected states. We can do dice throws by applying P to the right of our distribution vector.

>v=[1,0,0,0,0,0].P
0.4018776 0.4018776 0.1607510 0.0321502 0.0032150 0.0001286
% After two more throws, we get the following.
>v.P.P
0.0649055 0.2362559 0.3439886 0.2504237 0.0911542 0.0132721


This tells us that the chance to be in state S5 after three throws is just 1.3%.

How about the average time that we have to wait if we keep throwing the dice until we get a 6-Yahtzee? This can be solved by denoting the waiting time from state Si to S5 by w(i) and observe

$$w_i = p_{i,n} + \sum_{j \ne n} p_{i,j} (1+w_j)$$

where n=5 is the final state. Obviously, w(5)=0. Moreover, the sum of all probabilities in one row is 1. Thus we get

$$w_i – \sum_{j=1}^{n-1} w_j = 1$$

Let us solve this system in Euler Math Toolbox.

>B=P[1:5,1:5];
>w=(id(5)-B)\ones(5,1)
13.0236615
11.9266963
10.5554446
8.7272727
6.0000000


The average waiting time for a 6-Yahtzee is approximately 13 throws. If you already got 4 sixes, the average waiting time for the fifth is indeed 6.

We can interpret this in another way. Observe

$$w = (I-B)^{-1} \cdot 1 = (I+B+B^2+\ldots) \cdot 1$$

The sum converges because the norm of B is clearly less than 1. So we interpret this as a process

$$w_{n+1} = B \cdot (w_n + 1)$$

Suppose a waiting room. We have w(n,i) persons in state Si at time n. Now, one more person arrives at each state and we process the states according to our probabilities. On the long run, the number of persons in state S0 will become the waiting time to get out of the system into the final state.

Unfortunately, the probabilities are not that easy to compute when we ask for a Yahtzee of any kind. E.g., if we have 66432 we will keep the two sixes as 66. But we may throw 66555 and switch to keep 555. There is a recursion for the probability to get from i same dice to j same dice. I computed the following matrix which replaces the matrix P above.

 >P
0.0925926 0.6944444 0.1929012 0.0192901 0.0007716
0.0000000 0.5555556 0.3703704 0.0694444 0.0046296
0.0000000 0.0000000 0.6944444 0.2777778 0.0277778
0.0000000 0.0000000 0.0000000 0.8333333 0.1666667
0.0000000 0.0000000 0.0000000 0.0000000 1.0000000


Note that I have only the states S1 to S5 here. With the same method, I get the waiting times.

>B=PM[1:4,1:4];
>w=(id(4)-B)\ones(4,1)
11.0901554
10.4602273
8.7272727
6.0000000


Thus we need only about 11 throws to get any Yahtzee. The probability to get one in three throws is now about 4.6%.

>[1,0,0,0,0].PM.PM.PM
0.0007938 0.2560109 0.4524017 0.2447649 0.0460286


# The „One More Ace?“ Problem

Let me tackle another problem of statistics. I found it on a web page among other problems of the silly kind like „Is 0.999..<1?“. But actually, this one is a real problem that can provide insight into statistical thinking. It is elementary though. I am going to use Euler Math Toolbox (EMT) for some calculations and simulations.

Assume we have a set of Bridge cards (52 cards with 4 colors from 2 to 10, J, Q, K, Ace). This is dealt to four players, each player getting 13 cards. Answer the following two questions.

1. Assume one player gets an ace of any color. What is the probability to get a second ace?
2. Assume one player gets the ace of spades. What is the probability to get another ace?

My take on this is the following guideline:

Device a repeating experiment that mimics your problem!

In this case, we deal the 52 cards „randomly“ over and over again. We assume that every deal has the same probability. And I assume that you know that there are

$$\displaystyle \binom{52}{13} = \frac{52\cdot\ldots\cdot40}{2\cdot\ldots\cdot13} = 635’013’559’600$$

possible ways to select 13 cards from 52 cards. We have assumed that all of these deals come with the same probability. The frequentistic approach is to assume that we deal each distribution with the same frequency on the long (very long) run. So we can just assume we deal each one exactly once, and get our probability (the expected ratio of „good cases“) by

probability = „number of good cases“ / „number of all cases“

For problem 1, „all cases“ are the cases where there is one ace, and „good cases“ are the cases where there are two or more aces.  For problem 2, „all cases“ are the cases with the ace of spaces and „good cases“ the cases with the ace of spades and another ace.

We are then faced with counting, a problem of combinatorics. In this case, we can count the result. In other cases, there is only the option of a Monte-Carlo simulation. Let us do one in EMT.

>function simulate1 (N) ...
$c=1:52; n1=0; n2=0;$  loop 1 to N
$c=shuffle(c);$     k=sum(c[1:13]<=4);
$if k>=1 then$         n1=n1+1;
$if k>=2 then n2=n2+1; endif;$     endif;
$end;$   return n2/n1;
$endfunction >simulate1(1000000) 0.370341025387  EMT may be an easy language to do the programming, but it is not the fastest one. And it is easy only if you know the basic syntax of a matrix language. E.g., the command c[1:13]<=4 returns an array b of 13 zeros and ones (where we simulate the four aces by the numbers 1 to 4). Its element b[i] will be one if and only if c[i]<=4. Summing up we get the number of elements in the first 13 shuffled cards that are smaller or equal to 4. In the code, n1 is then the number of times where one ace showed up and n2 the number times where two showed up. It is important to note that we discarded all shuffles with no ace in the first 13 cards. It is like a condition in Bayesian reasoning. We are asking for the probability that there are two or more aces under the condition that there is one ace. The change for the second problem is only in one line. We have to check if a specific ace (say the number 1) is in the array. >function simulate2 (N) ...$  c=1:52; n1=0; n2=0;
$loop 1 to N$     c=shuffle(c);
$k=sum(c[1:13]<=4);$     if any(c[1:13]==1) then
$n1=n1+1;$         if k>=2 then n2=n2+1; endif;
$endif;$   end;
$return n2/n1;$  endfunction
>simulate2(1000000)
0.561442110359


The function any() returns 1 if any element of the argument is non-zero. Since c[1:13]==1 tests for the elements to be equal to 1 (the ace of spaces) we get the right answer.

As you see, the outcome of the experiment is quite different. This is why the problem is interesting.

Now we want to do the combinatorial job of counting and computing the probabilities. We expect our simulations to get „close“ to the „correct“ probabilities. The mathematics behind this is very basic combinatorics. My feeling is that this kind of reasoning should be in the mathematical class of any higher education. If you do not know the result you should learn it.

The number N(k) of deals for player number 1 which contain exactly k aces is as follows.

$$\displaystyle N_k = \binom{4}{k}\binom{48}{13-k}$$

The reasoning behind this is as follows: We have to select k from the 4 aces, and 13-k cards from the 48 cards that are no aces. If we do this we get all possible selections with 13 cards that contain exactly k aces.

Finally, the probability in problem 1 is („good cases / all couting cases“)

$$\displaystyle p_1 = \frac{N_2+N_3+N_4}{N_1+N_2+N_3+N_4} \approx 0.3696$$

as computed by EMT with

>function N(k) := bin(4,k)*bin(48,13-k)
>sum(N(2:4))/sum(N(1:4))
0.369637191337


This is close enough to our Monte-Carlo simulation.

Note that N(2:4) returns the array [N(2), N(3), N(4)] automatically in a matrix language. This works since the function bin() „maps“ to elements of array arguments. The sum of N(0) to N(4) is the number N of all deals. This allows a simplification with a new insight.

$$\displaystyle p_1 = 1 – \frac{N_1}{N-N_0} = 1 – \frac{a_1}{1-a_0}$$

where a(k)=N(k)/N is the probability to get exactly k aces. Of course, 1-a(0) is the probability got the one or more aces. So the right-hand side is what we expect as the probability for our Problem 1. If you are not as frequentistic as I am you may have come up with this formula immediately.

Let us tackle the Problem 2 now. This is even easier. We are now only dealing with 51 cards (all but the ace of spades), and we deal only 12 cards. It is easy to count the number of deals with none of the remaining aces. The answer is

$$\displaystyle p_2 = 1 – \frac{\binom{48}{12}}{\binom{51}{12}} \approx 0.56115$$

This should now be obvious (after some thought), and our simulation was close enough to this.

So the probabilities are indeed different. This seems counter-intuitive at first thought. But let us try to find a good argument.

Let us look at the complementary probabilities 1-p(1) and 1-p(2). In both problems, we then have to divide the number of times we have exactly one ace („bad cases“) over the number of times we have at least one ace („all counting cases“). But both numbers change between Problem 1 and 2. The bad cases in Problem 1 are simply 4 times the bad cases in Problem 2. But the counting cases in Problem 1 are less than 4 times the counting cases of Problem 2. It is not 4 times as likely to have any 3 aces than to have the ace of spades plus any two aces. Since the numerator is smaller, 1-p(1) is larger, and p(1) is consequently smaller.

# Povray and Euler

I found a nice plot by Taha on my Google+ page. It uses a software called MathMod. Euler Math Toolbox can do this too teaming with Povray.

>load povray;
>povstart(zoom=4,distance=6,center=[1,0,0.5],height=45°);
>t=linspace(0,1,100); s=t';
>b=cos(2pi*s)*(t+0.01)+t; c=sin(2pi*s)*(t-0.01); a=t;
>X=cos(2pi*a)*(b+0.1); Y=sin(2pi*a)*(b+0.1); Z=c;
>writeln(povgrid(X,Y,Z,d=0.03,dballs=0,skip=[20,10]));
>povend();


This is all that is needed to generate the following spiraling tube.

Let me explain.

First of all, the pov…() commands generate text snippets which model graphical elements in the syntax of Povray. These text elements are written to a file by writeln(). The commands povstart() and povend() simply open and close the file and start the Povray interpreter. The output of this interpreter is then inserted into the notebook window of EMT.

For an example of such a Povray snippet, let us try the following.

 >povbox([0,0,0],[1,2,3],povlook(green))
box { <0,0,0>, <1,2,3>
texture { pigment { color rgb <0.0627451,0.564706,0.0627451> }  }
finish { ambient 0.2 }
}


This is a green box with side lengths 1,2,3.

The povgrid() command finally is able to draw the surface as a grid. We use the skip= parameter to skip most grid lines from drawing, but still get a smooth appearance of the grid.

The coordinates for the mantle surface of the tube are X,Y,Z. These three variables contain matrices. Thus X[i,j],Y[i,j],Z[i,j] is one point on the surface. To compute these matrices, we use the matrix language of EMT. First, s and t are two vectors, a row and a column vector with values from 0 to 1. Then we generate a cylindrically shaped surface with coordinates a,b,c.

>povstart(zoom=4,distance=6,angle=80°,center=[0,1,0]);
>t=linspace(0,1,100); s=t';
>b=cos(2pi*s)*(t+0.01)+t; c=sin(2pi*s)*(t+0.01); a=t;
>for i=1 to 3; writeAxis(0,1.2,axis=i,c=lightgreen); end;
>writeln(povgrid(a,b,c,d=0.03,dballs=0,skip=[50,50]));
>povend();

This cylindrical surface is then wrapped around the z-axis (point upwards). The details are a bit involved. But any mathematician should be able to understand what is going on here.

# Translating a Sage Example for EMT

The nice Sage Math project is something I am very interested in. It is a rather complete package for mathematical computations based on Python. If you are interested in open source math this is probably one good way to start. I must add however that I prefer open source packages directly based on Python with a command line (like IPython) instead of a web based interface that covers everything with another layer. The advantage of the latter is, of course, that it may provide additional services, e.g., for cooperation.

On Windows you will have to run Sage Math in a Virtual Box. You can then open the web interface in your normal Windows browser and work there.

I will certainly come back in this blog to Sage Math. For a start, let me translate an example of Sage Math into Euler Math Toolbox (EMT). It is an example for the famous Lotka-Volterra equations, representing the relation between predator and prey. As far as I can tell, Python uses a fourth order Runge-Kutta method to solve the differential equation numerically. The equations themselves are written in symbolic form. The package contains features for symbolic computations.

### Simulates Lotka-Volterra population dynamics, producing a vector field
### Arrows indicate the direction of population evolution
### The blue contour is one potential stable loop
var('t R F R0 F0 a b c d')
## CONSTANTS
a = 0.04         # rabbit birthrate
b = 0.0005       # probability of predation per encounter
c = 0.1*0.0005   # rabbit -> fox conversion efficiency
d = 0.2          # death rate of foxes
## INITIAL CONDITIONS
R0 = 5000     #initial number of rabbits
F0 = 200      #initial number of foxes

## DIFFERENTIAL EQUATIONS
de_R = (diff(R,t) == a*R - b*R*F)
de_F = (diff(F,t) == c*R*F - d*F)

## CALCULATION PARAMETERS
end_points = 500
stepsize = 1.0
steps = end_points/stepsize
ics = [0,R0,F0]
des = [de_R.rhs(), de_F.rhs()]

## NUMERICAL SOLUTION OF DIFFERENTIAL EQUATIONS
sol = desolve_system_rk4(des,[R,F],ics,ivar=t,
end_points=end_points,step=stepsize)

## Clean up to graph
sol_t=[]; sol_R=[]; sol_F=[]
for i in range(steps):
sol_t.append(sol[i][0])
sol_R.append(sol[i][1])
sol_F.append(sol[i][2])

a = plot([],figsize=[6,4])
a += line(zip(sol_R,sol_F))
a += plot_vector_field((des[0], des[1]), (R,0,7000), (F,0,225),
xmin=1500,color='orange')
a.axes_labels(['Rabbits','Foxes']); show(a)


I entered that as is into the Sage Math browser and the following plot was produced. I had to copy the image out of the browser, because I have no idea on how to save the image with the browser interface. Let me also note that I do not understand every command above, though many of them are self explanatory.

Let us translate this too Euler Math Toolbox. We use the numerical side of EMT only, although we could have defined some functions in symbolic form. But it would not help much to do so.

>a=0.04; b=0.0005; c=0.1*b; d=0.2;
>R0=5000; F0=200;
>function dR ([R,F]) := a*R-b*R*F;
>function dF ([R,F]) := c*R*F-d*F;
>function f(x,y) := [dR(y),dF(y)];
>t=0:1:500;
>s=ode("f",t,[R0,F0]);
>{x,y}=vectorfield2("dR","dF",min(s[1]),max(s[1]),min(s[2]),max(s[2]), ...
>   scale=0.005,<plot);
>plot2d(s[1],s[2],color=blue, ...
>   xl="Rabbits",yl="Foxes");


The result is the following, very similar, plot. However, the anti-aliasing is better, but that could be due to my lack of knowledge in Sage Math.

# Python in EMT – Cross Sum of Cubes

When is the cross sum (sum of the digits) of the cube of a number equal to the number?

It is obviously true for n=0 and n=1. You might be able to find n=8 with the cube 512. Then, how do you prove that there are finitely many numbers? How long will it take to find all numbers?

For a dirty estimate we have

$$(9m)^3 \ge (a_0+\ldots+a_{m-1})^3 = a_0+a_1 10 + \ldots + a_{m-1} 10^{m-1} \ge 10^{m-1}$$

It is an easy exercise to show that m<7 is necessary. So we only have to check up to n=100. The cube of 100 has already 7 digits. This can be done by hand with some effort, but here is a Python program written in Euler Math Toolbox.

>function python qs (n) ...
$s=0$while n>0:
$s+=n%10;$   n/=10
$return s$endfunction
>function python test() ...
$v=[]$for k in range(100):
$if k == qs(k**3):$       v.append(k)
$return v$endfunction
>test()
[0,  1,  8,  17,  18,  26,  27]