Mathematica codeΒΆ

As presented in our paper, a staged tree statistical model can be characterised as the \(2\times 2\)-minors of its stage matrices. Whenever these minors are binomials, the staged tree has toric structure. However, sometimes a given parametrisation does not immediately give rise to binomials and these can only be found after an adequate reparametrisation. Even for medium sized examples, it is computationally infeasible to try out all such possible reparametrisations. We here thus follow an approach which employs random transformations on the stage matrices to overcome this computational obstacle and still provide an idea of whether a given staged tree might have toric structure. The code is essentially comprised of two functions: one which randomly generates invertible matrices with pre-specified entries and dimensions, and one which for every stage uses these matrices to create linear transformations, and subsequently checks whether these have the correct number of distinct entries (providing an isomorphism between the original specification and the new variables).

First, define a function which creates a list of invertible matrices with integer entries. This takes as a first argument a range in which to pick entries, passed to the function as a list, e.g. \(\verb|{-1,0,1,2}|\). The second argument is the dimension of the matrix, e.g. \(\verb|2|\) for a \(2\times 2\) matrix. The third argument is the length of the output list, e.g. \(\verb|100000|\). In the main body of this function, we first initalise the local variables as empty lists and lists of zeros. We then run a \(\verb|While|\) loop which uses the inbuilt \(\verb|Mathematica|\) function \(\verb|RandomChoice|\) to generate matrices with the specified entries and dimension, calculates their determinants, and deletes those which are singular until we end up with a list of the desired length.

invertibleinteger[range_, dims_, iterations_] := Module[
{integermatrices, determinants, length, finallist},

integermatrices = ConstantArray[{}, iterations];
determinants = ConstantArray[0, iterations];
length = 0;
finallist = {};

While[length < iterations,
        integermatrices = Table[RandomChoice[range, {dims, dims}], iterations];
        determinants = Map[Det, integermatrices];
        finallist = Join[finallist,Pick[integermatrices, Not /@ PossibleZeroQ[determinants]]];
        length = Length[finallist]

finallist = finallist[[1 ;; iterations]]

Second, define a function which uses the first to generate random linear transformations of a given list of matrices. These matrices are passed to the function as a first argument in the form of a list. They can, as in our application, have symbolic entries, e.g. \(\verb|{{{p1,p3},{p2,p4}},{{p1+p2,p1+p2+p3+p4,p1+p2+p3+p4+p5},{p3+p4,p5,p6}}}|\) for two \(2\times 2\) matrices representing binary stages in the first example we present. The second argument of this function is an integer specifying the number of distinct entries, or the number of variables in the joint parametrisation of the model, e.g. \(\verb|7|\). The third argument is a list of list pairs, fixing the range in which to look for transformations from the left and right, respectively, for each stage. E.g. \(\verb|ConstantArray[{Range[-1, 1], Range[-1, 1]}, 2]|\) picks transformations with entries in \(\{-1,0,1\}\) uniformly at random for both the transformations from the left and right for each stage. Weights can be specified here as an optional argument to modify the probability distribution underlying \(\verb|RandomChoice|\). The final argument is an integer specifying the number of random trials, e.g. \(\verb|100000|\). In the main body of the function, local variables are initialised as empty lists and lists of zeros. Then for each specified stage, \(\verb|lineartransform|\) calls \(\verb|invertibleinteger|\) to randomly create transformation matrices from the left and right and multiplies these by the stage matrix, finally using the inbuilt command \(\verb|Expand|\) (an alternative, albeit possibly slower, would be \(\verb|Simplify|\)) to ensure the symbolic expressions coming out of matrix multiplications are of the same form. After the \(\verb|For|\) loop, the transformed stage matrices are flattened so that the inbuilt command \(\verb|CountDistinct|\) can be applied to each transformation, checking whether the number of distinct entries equals the prespecified number. The final lines of code generate readable output, returning all involved variables for small examples with three or less iterations and a summary of the results for bigger ones.

lineartransform[stages_, distinctentries_, ranges_, iterations_] :=
{transleft, transright, stagestransformed, alltransformations, founddistinctentries, validstagetransformations, i, j, k},

transleft = ConstantArray[{}, Length[stages]];
transright = ConstantArray[{}, Length[stages]];
stagestransformed = ConstantArray[{}, Length[stages]];

For[i = 1, i <= Length[stages], i++,
        transleft[[i]] = invertibleinteger[ranges[[i]][[1]], Dimensions[stages[[i]]][[1]], iterations];
        transright[[i]] = invertibleinteger[ranges[[i]][[2]], Dimensions[stages[[i]]][[2]], iterations];
        stagestransformed[[i]] = Expand[MapThread[Dot, {transleft[[i]], ConstantArray[stages[[i]], iterations], transright[[i]]}]]

alltransformations = Map[Flatten, MapThread[Join, stagestransformed]];
founddistinctentries = Map[CountDistinct, alltransformations];
validstagetransformations = Pick[alltransformations, PossibleZeroQ[founddistinctentries - distinctentries]];

If[iterations <= 3,

For[j = 1, j <= Length[stages], j++,
Print["-STAGE ", j, "-
transformation matrices from the left: ", transleft[[j]], "
transformation matrices from the right: ", transright[[j]], "
transformations of the stage as matrices: ",
all transformations of the stages flat: ", alltransformations, "
number of distinct entries: ", founddistinctentries, "
and valid stagetransformations: ", validstagetransformations],

Print[iterations, " iterations performed. A total of ", CountDistinct[alltransformations], " different transformations have been checked.
        These had (weighted) entries in ", ranges, " for the left and right transformations of the respective ", Length[stages], " stages."];

If[Length[validstagetransformations] != 0,
Print[Length[validstagetransformations], " correct transformations found! The first ", Min[10, Length[validstagetransformations]], " of these
        are: ", validstagetransformations[[1 ;; Min[10, Length[validstagetransformations]]]]],
Print["No valid transformations with ", distinctentries, " distint entries found. ",
        Table[Length[Pick[alltransformations,PossibleZeroQ[founddistinctentries - k]]], {k,Length[Flatten[stages]]}], " of these had, respectively, ",
        Range[Length[Flatten[stages]]], " distinct entries."]]