vars = {p0000, p0001, p0010, p0011, p0100, p0101, p0110, p0111, p1000, p1001, p1010, p1011, p1100, p1101, p1110, p1111};
S1 = Simplify[-p0010*p1000-p0110*p1000+p0000*p1010+p0100*p1010-p0010*p1100-p0110*p1100+p0000*p1110+p0100*p1110==0&&-p0011*p1001-p0111*p1001+p0001*p1011+p0101*p1011-p0011*p1101-p0111*p1101+p0001*p1111+p0101*p1111==0]
S2 = Simplify[-p0010*p0100+p0000*p0110+p0110*p1000-p0100*p1010-p0010*p1100-p1010*p1100+p0000*p1110+p1000*p1110==0&&-p0011*p0101+p0001*p0111+p0111*p1001-p0101*p1011-p0011*p1101-p1011*p1101+p0001*p1111+p1001*p1111==0]

RandomDistribution[] := Module[{pvars, n=Length[vars], k, l, m, Avars, Bvars, Cvars, Dvars},
pvars = Permute[vars,RandomPermutation[n]];
k = RandomInteger[{1, n/2 - 2}];
l = RandomInteger[{1, n/2 - 1 - k}];
m = RandomInteger[{1, n/2 - k - l}];
Avars = pvars[[1;; k]];
Bvars = pvars[[k+1;; k+l]];
Cvars = pvars[[k+l+1;; k+l+m]];
Dvars = pvars[[k+l+m+1;; n]];
Catenate[{
Map[# -> 1/(k+l) &, Avars],
Map[# -> 1/(k+l) - ϵ &, Bvars],
Map[# -> l/m ϵ &, Cvars],
Map[# -> 0 &,Dvars]
}]
]

RandomGoodDistribution[] := Module[{p},
While[True, p = RandomDistribution[];
If[Not[And @@ {S1, S2} //. Catenate[{p, { ϵ -> 10^-10 }}]] && And @@ {S1, S2} //. Catenate[{p, { ϵ -> 0 }}],Break[]]];
p
]

EntropyExpr[P_] := -(Plus @@ Map[# Log[#] &, Select[P, Not[PossibleZeroQ[#]]&]]);
hA[P_] := EntropyExpr[{p0000+p0001+p0010+p0011+p0100+p0101+p0110+p0111,p1000+p1001+p1010+p1011+p1100+p1101+p1110+p1111} //. P];
hB[P_] := EntropyExpr[{p0000+p0001+p0010+p0011+p1000+p1001+p1010+p1011,p0100+p0101+p0110+p0111+p1100+p1101+p1110+p1111} //. P];
hC[P_] := EntropyExpr[{p0000+p0001+p0100+p0101+p1000+p1001+p1100+p1101,p0010+p0011+p0110+p0111+p1010+p1011+p1110+p1111} //. P];
hD[P_] := EntropyExpr[{p0000+p0010+p0100+p0110+p1000+p1010+p1100+p1110,p0001+p0011+p0101+p0111+p1001+p1011+p1101+p1111} //. P];
hAB[P_] := EntropyExpr[{p0000+p0001+p0010+p0011,p0100+p0101+p0110+p0111,p1000+p1001+p1010+p1011,p1100+p1101+p1110+p1111} //. P];
hAC[P_] := EntropyExpr[{p0000+p0001+p0100+p0101,p0010+p0011+p0110+p0111,p1000+p1001+p1100+p1101,p1010+p1011+p1110+p1111} //. P];
hAD[P_] := EntropyExpr[{p0000+p0010+p0100+p0110,p0001+p0011+p0101+p0111,p1000+p1010+p1100+p1110,p1001+p1011+p1101+p1111} //. P];
hBC[P_] := EntropyExpr[{p0000+p0001+p1000+p1001,p0010+p0011+p1010+p1011,p0100+p0101+p1100+p1101,p0110+p0111+p1110+p1111} //. P];
hBD[P_] := EntropyExpr[{p0000+p0010+p1000+p1010,p0001+p0011+p1001+p1011,p0100+p0110+p1100+p1110,p0101+p0111+p1101+p1111} //. P];
hCD[P_] := EntropyExpr[{p0000+p0100+p1000+p1100,p0001+p0101+p1001+p1101,p0010+p0110+p1010+p1110,p0011+p0111+p1011+p1111} //. P];
hABC[P_] := EntropyExpr[{p0000+p0001,p0010+p0011,p0100+p0101,p0110+p0111,p1000+p1001,p1010+p1011,p1100+p1101,p1110+p1111} //. P];
hABD[P_] := EntropyExpr[{p0000+p0010,p0001+p0011,p0100+p0110,p0101+p0111,p1000+p1010,p1001+p1011,p1100+p1110,p1101+p1111} //. P];
hACD[P_] := EntropyExpr[{p0000+p0100,p0001+p0101,p0010+p0110,p0011+p0111,p1000+p1100,p1001+p1101,p1010+p1110,p1011+p1111} //. P];
hBCD[P_] := EntropyExpr[{p0000+p1000,p0001+p1001,p0010+p1010,p0011+p1011,p0100+p1100,p0101+p1101,p0110+p1110,p0111+p1111} //. P];
hABCD[P_] := EntropyExpr[{p0000,p0100,p0001,p0101,p0010,p0110,p0011,p0111,p1000,p1100,p1001,p1101,p1010,p1110,p1011,p1111} //. P];

iCDA[P_] := hAC[P] + hAD[P] - hACD[P] - hA[P];
iCDB[P_] := hBC[P] + hBD[P] - hBCD[P] - hB[P];
iAB[P_] := hA[P] + hB[P] - hAB[P];
iCD[P_] := hC[P] + hD[P] - hCD[P];
Ingleton[P_] := iCDA[P] + iCDB[P] + iAB[P] - iCD[P];

iACD[P_] := hAD[P] + hCD[P] - hACD[P] - hD[P];
iBCD[P_] := hBD[P] + hCD[P] - hBCD[P] - hD[P];

(* Return coefficient and exponent of the lowest-order term in the power series expansion of F in x around 0. *)
SeriesOrder[F_, x_] := Module[{c=Series[Simplify[F], (*x -> 0*) {x,0,4}]},
(* Series may return a SeriesData which we can index, or a constant. *)
If[SameQ[Head[c], SeriesData],
{ c[[3,1]], c[[4]]/c[[6]] },
If[c==0, {0, Infinity},{c, 0}]
]
];
Lim[F_, x_] := Limit[F, {x -> 0}, Direction -> "FromAbove"];

SelectFirst[Range[1000000], Module[{p, ing, acd, bcd},
p = RandomGoodDistribution[];
ing = SeriesOrder[Ingleton[p], ϵ];
acd = SeriesOrder[iACD[p], ϵ];
bcd = SeriesOrder[iBCD[p], ϵ];
TimeConstrained[If[ing[[2]] > 0 && Lim[ing[[1]], ϵ] < 0 &&
((ing[[2]] < Min[acd[[2]], bcd[[2]]]) ||
(ing[[2]] == Min[acd[[2]], bcd[[2]]] && Lim[ing[[1]], ϵ] == -Infinity && Abs[Lim[acd[[1]], ϵ]] < Infinity && Abs[Lim[bcd[[1]], ϵ]] < Infinity)),
Print[p, " is good"]; True],
60,Print[p, " is unknown"]; False
]]&
];

(* This is found eventually: *)
p = {p1111->1/5,p1101->1/5,p0110->1/5,p1010->1/5,p0010->1/5-ϵ,p1100->ϵ,p0100->0,p0001->0,p0111->0,p0000->0,p1011->0,p0011->0,p0101->0,p1001->0,p1000->0,p1110->0};

Series[Ingleton[p], {ϵ,0,3}]
Simplify[iACD[p] == iBCD[p]]
Series[iACD[p], {ϵ,0,3}]