numWorlds = 4;

worldIndices = Range[numWorlds]

ithWorld[i_] := Normal[SparseArray[{i -> 1}, numWorlds]]

Worlds = Map[ithWorld, worldIndices]

I could also have specified Worlds as IdentityMatrix[numWorlds], but chose this because it's slightly more illuminating.

cubeCheck[x_] := Fold[ And,True,Map[Function[y, 0≤y≤ 1], x]]

probCheck[x_] := Length[x] == numWorlds && cubeCheck[x] && Total[x] == 1

incoherent[x_] := cubeCheck[x] && !probCheck[x]

Although it' s easy enough to come up with your own examples of incoherent points, you can also use Mathematica to generate more, using FindInstance as below.  I’m not sure how to use FindInstance (or anything else) to generate a list of the right dimension instead of supplying it a list of variables.

FindInstance[incoherent[{x,y,z}],{x,y,z}, Reals, 4]

c = {.2, .2, .2,.2}

incoherent[c]

ithCoherent[i_,c_] := ReplacePart[c, {i -> 1 -(Total[c]-c[[i]])}]

coherentPoints = Map[Function[i,ithCoherent[i,c]], worldIndices]

Now, Joyce' s construction makes use of an arbitrary scoring rule satisfying his assumptions.  Since our goal is not to prove the theorem, but to see how it works in various examples, I will use squared Euclidean distance to start.  To try others, just replace SquaredEuclideanDistance in the definition of f here with your favorite scoring rule.

f[i_,b_] := SquaredEuclideanDistance[b,Worlds[[i]]]-SquaredEuclideanDistance[c,Worlds[[i]]]

M[p_]:=1/numWorlds*Sum[f[i,p],{i,numWorlds}]

g[i_,p_]:=Max[0,f[i,p]-M[p]]

totalG[p_]:=Sum[g[i,p],{i,numWorlds}]

mix[i_,p_]:=g[i,p]/(1+totalG[p])

mixN[p_]:=1/(1+totalG[p])

G[p_]:=Sum[mix[i,p]*coherentPoints[[i]],{i,numWorlds}]+mixN[p]*p

Now that the continuous function G on the convex hull of coherentPoints has been defined, we know it has a fixed point.  Let' s compute it!  The second argument just tells Mathematica where to start searching for a fixed point.

dominator = FixedPoint[G,coherentPoints[[1]]]

This function tests whether a given p strictly dominates the incoherent c.  Then we check that the dominator generated by Joyce does in fact dominate c.

dominates[p_]:=Fold[And,True,Map[Function[i, f[i,p]<0],worldIndices]]

dominates[dominator]