Here is a computation of the exact probability. This is a bit messy and takes several minutes of computation. It would be infeasible to check all 16! ~ 21 trillion permutations. But they can be grouped as fixed-size set partitions such that either all the permutations matching a particular partition are valid, or they're all invalid. The number of permutations covered by each partition is exactly (3!)^2 * (4!)^3 = 497,664. Thus the number of cases that need to be checked is reduced to about 42 million, which is feasible.

restart:
FamN:= [3,4,4,4,1]:  #family sizes
PS:= [seq[scan= `+`]](FamN);
                    PS := [3, 7, 11, 15, 16]

FamS:= (`{}`@`$`@`..`)~([1,(PS[..-2]+~1)[]], PS)[];  #the initial set partition
       FamS := {1, 2, 3}, {4, 5, 6, 7}, {8, 9, 10, 11}, {12, 13, 14, 15}, {16}

nf:= nops(FamN):
SPFS:= Iterator:-SetPartitionFixedSize(FamN):
sp:= output(SPFS):

# For efficiency, we divide the output Array from the Iterator into 5 chunks 
# matching the pattern of the partition. There's no mathematical reason to do
# this; it's only done to save about 4 minutes of "real" computation time.

FamA:= seq(ArrayTools:-Alias(sp, [0, PS[..-2][]][i], [FamN[i]]), i= 1..nf):
ValidDraw?:= ()-> local i,x;
    not orseq(orseq(x in FamS[i], x= FamA[i]), i= 1..nf)
:
V:= 0:  #count of valid partitions (no intra-family pairings)
CodeTools:-Usage(in SPFS do if ValidDraw?() then V++ fi od):
memory used=24.64GiB, alloc change=448.64MiB, 
cpu time=15.16m, real time=5.09m, gc time=12.08m

V;
                             751650
prob:= V/Number(SPFS);
                          prob := 5011 / 280280
evalf(prob);
                         0.01787855002

So, the exact value is extremely close to my simulated estimate,

@mmcdara Thank you for your input. Okay, I understand your nuance about the correct usage of p-value. Nonetheless, would you agree with this statement: "This simulation estimates that the probability that a random drawing would have no intra-family pairings is less than 2%"?

I just finished computing the exact value of the probability, and I'm about to post it. The estimate that I made above was very good, with relative error of only 0.006.

@igor.v.proskurin Sorry if I seemed harsh. I think that you've been a great contributor to MaplePrimes considering that you've only been here for two weeks.

@igor.v.proskurin You wrote:

• So when you are doing subs(x = 1, g(x)), it means diff(g(1), 1), (not diff(g(x), x) at x= 1) and that is the problem.

That is 100% incorrect. Like the vast majority of Maple commands, subs evaluates its arguments before using them. Since g(x) evaluates to something that has no diff, there's no difference between subs(x=1, g(x)) and eval(g(x), x=1) in this case. Indeed, the OP said that subs worked.

@mmcdara I have a quibble with your expected value computation. In your procedure, 3 of the throws occur exactly 1 time. It's only the 2nd throw that has a possibility of having 6s discarded. So, the expected value is

3 + 5*sum(n/6^n, n= 1..infinity) = 4.2 (exactly).

I've verified this experimentally by putting a global counter in procedure dice and running SelectCard() 100,000 times. The throw count was 419,922.

@ Please post again in this thread if you continue to have this problem. 

Here are simple arithmetic steps for your son:

To generate a random integer n  from 1 to 120 using a 6-sided die:

1. Roll the die. Let d be the number shown.
2. Let n:= 30 × (d - 1).
3. Roll the die, perhaps repeatedly, until it shows a number less than 6. Call this number d.
4. Let n:= n + 6 × (d - 1).
5. Roll the die until it shows a number less than 5. Call the number d.
6. Let n:= n + d. 
7. n is now your card number. If that card has already been chosen, restart at step 1.

@acer The OP's intent wasn't totally clear to me. I only meant to offset the possibility that the OP mistakenly thought that relational operators (used in the normal way: exactly 2 operands) were contributing to the Boolean evaluation problem.

@digerdiga You asked:

• Is it also possible to use a different symbol (=,<=,>=,...) at each step?

From your example, I think that you're more interested in changing to symbols with different mathematical meaning rather than merely changing the typography. All of the Maple operators {=, <=, >=, etc.} (called relational operators) are inherently inert; they didn't cause your original evaluation problem. It was only the and that forced a Boolean evaluation; so it's the only thing that needs an inertness or unevaluation adjustment. The and arises from using a relational operator with more than 2 operands (in 2-D input). For example, a = b = c is parsed to a = b and b = c.

@Hahn Hahn The strange derivative output that you showed is definitely caused by the absence of a multiplication sign between (c-x) and (q-e).

@C_R The only way to remove points is for the person who gave the points to reclick the same icon before the item is deleted. After the item is deleted, points can't be changed.

@Preben Alsholm I understood your point already, and, indeed, I almost posted a similiar Question, so I'm glad that you did. I think it would be good to let it stand for a few days so that others may understand also. You can delete it when you feel it's appropriate to do so.

Any Moderator can delete a thread, but only you can deselect Best Answer to remove the 7 points.

@Christian Wolinski I also disagree, strongly, with your reasoning. Some reasons:

• subsop is a builtin (i.e. kernel) command so it can postpone automatic simplification until it exits.
• subs is capable of doing simultaneous substitutions, as shown in my Answer.
• Regardless of who wrote the subsop help page, it's clear to me that the designed intention was simultaneous substitution. Even if there were an external help-page writer, throwing in the word "simultaneous" seems an unlikely mistake.
• You wrote, "a new expression with 2 elements is formed: a*(rest of expression). There is no operand 3 and 4." 
  If (rest of expression) is type `*`, as it is in this case, I don't think it'd be considered a single operand separated from a. (And certainly that's not possible if automatic simplification is happening between substitutions (which may be happening unintentionally, due to a bug).)

So it seems that your reasoning depends on both the presence and absence of automatic simplification happening between the individual substitutions. 

@srwarner I don't know what you're talking about in your latest Reply. There is no contradiction in the first two assume commands above.

@rzarouf acer used eps = 0.00001; you used eps = 0.1. That makes his pseudospectral regions much smaller than yours. Essentially, they're only visible as the eigenvalue itself.