In this vignette we show how to set up a model to perform inference for the generalised Stochastic Block Model (GSBM) using the reversible jump split-merge sampling method.
To start, load the library
In this toy example we consider a set of 50 interconnected electrical components. Each component is connected to the other and the connections are all turned on at once. We record the failure time of each connections.
This can be modelled as a network, where the edge-weights are the failure times. A standard likelihood model for failure times is the Weibull distribution (on Wikipedia). The probability density function is: \[ \frac{k}{\lambda}\left(\frac{x}{\lambda}\right)^{k-1} e^{-(x/\lambda)^{k}} \]
We posit that the components exhibit some group behaviour via the failure times, such that connections between groups tend to fail more quickly than connections within a group; however we don’t know the number of groups. The GSBM is a good model to consider for such data since it can provide a posterior distribution for the number of groups and the Weibull with both parameters unknown has no conjugate prior.
There are three parts to the GSBM: - An edge model - A block model - A parameter model (for the parameters of the edge model)
The edgemod object in SBMSplitMerge encapsulates the likelihood of the edge weights. A basic edgemod is a thin wrapper to a density function. The package insists on the density function taking the form f(x, p) where x is an edge-weight and p is a parameter 3-array. An optional random method can be added. This has the signature f(p) where p is again a vector of parameters. For the Weibull example, we can rely on ’s built in Weibull functions
In the GSBM, the prior specification for the number of blocks is explicit. To facilitate other models (such as a Chinese Restaurant Process), the  object in the package encapsulates the join distribution for the number of blocks and assignment of nodes to blocks; we refer to these as the ‘block structure’. To provide a  the following are needed: some input parameters gamma, log density for the block structure, a random method for the block structure, the conditional distribution of a single node assignment given the others and a boolean indicating if the number of blocks (kappa) is fixed. The package provides three pre-specified s: - multinom - fixed number of blocks with a multinomial assignment of nodes to blocks - crp - the Chinese restaurant process - dma - Dirichlet-Multinomial allocation
We will use dma, the first parameter (gamma) is the concentration parameter for a symmetric Dirichlet distribution and the second (delta) is the rate parameter for a Poisson. The ‘number of blocks’-1 is modelled as a Poisson(delta) and the block assignments as Dirichlet-Multinomial(gamma, K). [cf. the paper]
Finally, the prior model on the parameters of the Weibull distribution is specified. This is the most involved component since it requires some functions to map the inputs to the real line. The components of a parammod are: -  - a function(p) where p is a params object gets the probability density of a params object, -  - a function(kappa) simulates a params object for kappa blocks, -  - a function mapping the support of the parameter to the real line, -  - a function mapping the real line to the support of the parameter, -  - the jacobian of t.
For the Weibull example, \(k\) and \(\lambda\) must be positive. A reasonable prior for each is is a Gamma distribution, which can be specified like so:
param_model <- parammod(
        function(params){
                dgamma(params$theta0[1], 1, 1, log=TRUE) +
                        dgamma(params$theta0[2], 1, 1, log=TRUE) +
                        sum(dgamma(params$thetak[,1], 1, 1, log=TRUE)) +
                        sum(dgamma(params$thetak[,2], 1, 1, log=TRUE))
        },
        function(kappa){
                params(
                        c(rgamma(1, 1, 1), rgamma(1, 1, 1))
                ,
                        cbind(rgamma(kappa, 1, 1), rgamma(kappa, 1, 1))
                )
        },
        function(x){ cbind(log(x[1]), log(x[2]))},
        function(x){ cbind(exp(x[1]), exp(x[2]))},
        function(x){ -log(x[1])-log(x[2])}
        )this prior says each \(k\) and \(\lambda\) are Gamma(1,1) distributed.
We wrap all the model components up into a sbmmod object:
The data should appear in a square matrix and wrapped in an edges object. For the purpose of this example we simulate the data with some true values:
set.seed(1)
true_blocks   <- blocks(rep(c(1, 2, 3), c(10, 20, 20)))
true_params   <- params(c(1, 0.5), cbind(1, c(3,4,5)))
true_sbm      <- sbm(true_blocks, true_params)
weibull_edges <- redges(true_sbm, model$edge)The blocks object describes the block assignments of nodes to blocks:
print(true_blocks)
#> blocks object
#> kappa:  3 
#> Number of nodes:  50 
#> block sizes:  10 20 20
plot(true_blocks)We are now in a position to sample from the model:
and see some diagnostic plots
#> 
#> $post_pairs#> 
#> $post_pairs_sorted#> 
#> $blocks_trace#> 
#> $param_trace
#> $param_trace[[1]]
#> Warning: Removed 14043 row(s) containing missing values (geom_path).#> 
#> $param_trace[[2]]
#> Warning: Removed 14043 row(s) containing missing values (geom_path).#> 
#> 
#> $param_trace_sorted
#> $param_trace_sorted[[1]]
#> Warning: Removed 14043 row(s) containing missing values (geom_path).#> 
#> $param_trace_sorted[[2]]
#> Warning: Removed 14043 row(s) containing missing values (geom_path).#> 
#> 
#> $blocks_trace_sorted#> 
#> $pp
#>            [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
#>  [1,] 0.0000000 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#>  [2,] 0.9966667 0.0000000 0.9933333 0.9966667 0.9966667 0.9966667 0.9966667
#>  [3,] 0.9966667 0.9933333 0.0000000 0.9933333 0.9933333 0.9933333 0.9933333
#>  [4,] 0.9966667 0.9966667 0.9933333 0.0000000 0.9966667 0.9966667 0.9966667
#>  [5,] 0.9966667 0.9966667 0.9933333 0.9966667 0.0000000 0.9966667 0.9966667
#>  [6,] 0.9966667 0.9966667 0.9933333 0.9966667 0.9966667 0.0000000 0.9966667
#>  [7,] 0.9966667 0.9966667 0.9933333 0.9966667 0.9966667 0.9966667 0.0000000
#>  [8,] 1.0000000 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#>  [9,] 0.9966667 0.9966667 0.9933333 0.9966667 1.0000000 0.9966667 0.9966667
#> [10,] 0.9966667 0.9966667 0.9933333 0.9966667 0.9966667 0.9966667 0.9966667
#> [11,] 0.7666667 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#> [12,] 0.7633333 0.7633333 0.7600000 0.7633333 0.7633333 0.7633333 0.7633333
#> [13,] 0.7666667 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#> [14,] 0.7633333 0.7633333 0.7600000 0.7633333 0.7633333 0.7633333 0.7666667
#> [15,] 0.7633333 0.7666667 0.7600000 0.7633333 0.7633333 0.7633333 0.7633333
#> [16,] 0.7633333 0.7633333 0.7600000 0.7666667 0.7633333 0.7633333 0.7633333
#> [17,] 0.7633333 0.7633333 0.7600000 0.7666667 0.7633333 0.7633333 0.7633333
#> [18,] 0.7633333 0.7633333 0.7600000 0.7633333 0.7633333 0.7633333 0.7666667
#> [19,] 0.7633333 0.7633333 0.7600000 0.7633333 0.7633333 0.7633333 0.7666667
#> [20,] 0.7633333 0.7633333 0.7600000 0.7633333 0.7666667 0.7633333 0.7633333
#> [21,] 0.7666667 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#> [22,] 0.7666667 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#> [23,] 0.7633333 0.7633333 0.7600000 0.7633333 0.7666667 0.7633333 0.7633333
#> [24,] 0.7633333 0.7666667 0.7600000 0.7633333 0.7633333 0.7633333 0.7633333
#> [25,] 0.7666667 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#> [26,] 0.7633333 0.7633333 0.7600000 0.7633333 0.7633333 0.7633333 0.7633333
#> [27,] 0.7633333 0.7666667 0.7600000 0.7633333 0.7633333 0.7633333 0.7633333
#> [28,] 0.7633333 0.7633333 0.7600000 0.7633333 0.7633333 0.7633333 0.7633333
#> [29,] 0.7666667 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#> [30,] 0.7633333 0.7633333 0.7600000 0.7633333 0.7633333 0.7633333 0.7666667
#> [31,] 0.7633333 0.7633333 0.7600000 0.7633333 0.7666667 0.7633333 0.7633333
#> [32,] 0.7633333 0.7633333 0.7600000 0.7633333 0.7666667 0.7633333 0.7633333
#> [33,] 0.7633333 0.7633333 0.7600000 0.7633333 0.7666667 0.7633333 0.7633333
#> [34,] 0.7633333 0.7633333 0.7600000 0.7633333 0.7633333 0.7633333 0.7666667
#> [35,] 0.7666667 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#> [36,] 0.7633333 0.7666667 0.7600000 0.7633333 0.7633333 0.7633333 0.7633333
#> [37,] 0.7633333 0.7633333 0.7600000 0.7633333 0.7633333 0.7633333 0.7633333
#> [38,] 0.7633333 0.7633333 0.7600000 0.7633333 0.7633333 0.7633333 0.7633333
#> [39,] 0.7633333 0.7633333 0.7600000 0.7633333 0.7633333 0.7633333 0.7633333
#> [40,] 0.7633333 0.7633333 0.7600000 0.7666667 0.7633333 0.7633333 0.7633333
#> [41,] 0.7633333 0.7633333 0.7600000 0.7633333 0.7666667 0.7633333 0.7633333
#> [42,] 0.7633333 0.7666667 0.7600000 0.7633333 0.7633333 0.7633333 0.7633333
#> [43,] 0.7633333 0.7633333 0.7600000 0.7666667 0.7633333 0.7633333 0.7633333
#> [44,] 0.7633333 0.7633333 0.7600000 0.7633333 0.7633333 0.7633333 0.7666667
#> [45,] 0.7633333 0.7633333 0.7600000 0.7633333 0.7633333 0.7633333 0.7666667
#> [46,] 0.7633333 0.7633333 0.7600000 0.7666667 0.7633333 0.7633333 0.7633333
#> [47,] 0.7633333 0.7633333 0.7600000 0.7666667 0.7633333 0.7633333 0.7633333
#> [48,] 0.7633333 0.7633333 0.7600000 0.7633333 0.7633333 0.7633333 0.7633333
#> [49,] 0.7633333 0.7633333 0.7600000 0.7666667 0.7633333 0.7633333 0.7633333
#> [50,] 0.7633333 0.7633333 0.7600000 0.7633333 0.7633333 0.7633333 0.7633333
#>            [,8]      [,9]     [,10]     [,11]     [,12]     [,13]     [,14]
#>  [1,] 1.0000000 0.9966667 0.9966667 0.7666667 0.7633333 0.7666667 0.7633333
#>  [2,] 0.9966667 0.9966667 0.9966667 0.7633333 0.7633333 0.7633333 0.7633333
#>  [3,] 0.9966667 0.9933333 0.9933333 0.7633333 0.7600000 0.7633333 0.7600000
#>  [4,] 0.9966667 0.9966667 0.9966667 0.7633333 0.7633333 0.7633333 0.7633333
#>  [5,] 0.9966667 1.0000000 0.9966667 0.7633333 0.7633333 0.7633333 0.7633333
#>  [6,] 0.9966667 0.9966667 0.9966667 0.7633333 0.7633333 0.7633333 0.7633333
#>  [7,] 0.9966667 0.9966667 0.9966667 0.7633333 0.7633333 0.7633333 0.7666667
#>  [8,] 0.0000000 0.9966667 0.9966667 0.7666667 0.7633333 0.7666667 0.7633333
#>  [9,] 0.9966667 0.0000000 0.9966667 0.7633333 0.7633333 0.7633333 0.7633333
#> [10,] 0.9966667 0.9966667 0.0000000 0.7633333 0.7633333 0.7633333 0.7633333
#> [11,] 0.7666667 0.7633333 0.7633333 0.0000000 0.9966667 1.0000000 0.9966667
#> [12,] 0.7633333 0.7633333 0.7633333 0.9966667 0.0000000 0.9966667 0.9966667
#> [13,] 0.7666667 0.7633333 0.7633333 1.0000000 0.9966667 0.0000000 0.9966667
#> [14,] 0.7633333 0.7633333 0.7633333 0.9966667 0.9966667 0.9966667 0.0000000
#> [15,] 0.7633333 0.7633333 0.7633333 0.9966667 0.9966667 0.9966667 0.9966667
#> [16,] 0.7633333 0.7633333 0.7633333 0.9966667 0.9966667 0.9966667 0.9966667
#> [17,] 0.7633333 0.7633333 0.7633333 0.9966667 0.9966667 0.9966667 0.9966667
#> [18,] 0.7633333 0.7633333 0.7633333 0.9966667 0.9966667 0.9966667 1.0000000
#> [19,] 0.7633333 0.7633333 0.7633333 0.9966667 0.9966667 0.9966667 1.0000000
#> [20,] 0.7633333 0.7666667 0.7633333 0.9966667 0.9966667 0.9966667 0.9966667
#> [21,] 0.7666667 0.7633333 0.7633333 1.0000000 0.9966667 1.0000000 0.9966667
#> [22,] 0.7666667 0.7633333 0.7633333 1.0000000 0.9966667 1.0000000 0.9966667
#> [23,] 0.7633333 0.7666667 0.7633333 0.9966667 0.9966667 0.9966667 0.9966667
#> [24,] 0.7633333 0.7633333 0.7633333 0.9966667 0.9966667 0.9966667 0.9966667
#> [25,] 0.7666667 0.7633333 0.7633333 1.0000000 0.9966667 1.0000000 0.9966667
#> [26,] 0.7633333 0.7633333 0.7666667 0.9966667 0.9966667 0.9966667 0.9966667
#> [27,] 0.7633333 0.7633333 0.7633333 0.9966667 0.9966667 0.9966667 0.9966667
#> [28,] 0.7633333 0.7633333 0.7666667 0.9966667 0.9966667 0.9966667 0.9966667
#> [29,] 0.7666667 0.7633333 0.7633333 1.0000000 0.9966667 1.0000000 0.9966667
#> [30,] 0.7633333 0.7633333 0.7633333 0.9966667 0.9966667 0.9966667 1.0000000
#> [31,] 0.7633333 0.7666667 0.7633333 0.8366667 0.8366667 0.8366667 0.8366667
#> [32,] 0.7633333 0.7666667 0.7633333 0.8366667 0.8366667 0.8366667 0.8366667
#> [33,] 0.7633333 0.7666667 0.7633333 0.8366667 0.8366667 0.8366667 0.8366667
#> [34,] 0.7633333 0.7633333 0.7633333 0.8366667 0.8366667 0.8366667 0.8400000
#> [35,] 0.7666667 0.7633333 0.7633333 0.8400000 0.8366667 0.8400000 0.8366667
#> [36,] 0.7633333 0.7633333 0.7633333 0.8366667 0.8366667 0.8366667 0.8366667
#> [37,] 0.7633333 0.7633333 0.7633333 0.8366667 0.8366667 0.8366667 0.8366667
#> [38,] 0.7633333 0.7633333 0.7633333 0.8366667 0.8400000 0.8366667 0.8366667
#> [39,] 0.7633333 0.7633333 0.7633333 0.8366667 0.8366667 0.8366667 0.8366667
#> [40,] 0.7633333 0.7633333 0.7633333 0.8366667 0.8366667 0.8366667 0.8366667
#> [41,] 0.7633333 0.7666667 0.7633333 0.8366667 0.8366667 0.8366667 0.8366667
#> [42,] 0.7633333 0.7633333 0.7633333 0.8366667 0.8366667 0.8366667 0.8366667
#> [43,] 0.7633333 0.7633333 0.7633333 0.8366667 0.8366667 0.8366667 0.8366667
#> [44,] 0.7633333 0.7633333 0.7633333 0.8366667 0.8366667 0.8366667 0.8400000
#> [45,] 0.7633333 0.7633333 0.7633333 0.8366667 0.8366667 0.8366667 0.8400000
#> [46,] 0.7633333 0.7633333 0.7633333 0.8366667 0.8366667 0.8366667 0.8366667
#> [47,] 0.7633333 0.7633333 0.7633333 0.8366667 0.8366667 0.8366667 0.8366667
#> [48,] 0.7633333 0.7633333 0.7633333 0.8366667 0.8400000 0.8366667 0.8366667
#> [49,] 0.7633333 0.7633333 0.7633333 0.8366667 0.8366667 0.8366667 0.8366667
#> [50,] 0.7633333 0.7633333 0.7633333 0.8366667 0.8366667 0.8366667 0.8366667
#>           [,15]     [,16]     [,17]     [,18]     [,19]     [,20]     [,21]
#>  [1,] 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7666667
#>  [2,] 0.7666667 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#>  [3,] 0.7600000 0.7600000 0.7600000 0.7600000 0.7600000 0.7600000 0.7633333
#>  [4,] 0.7633333 0.7666667 0.7666667 0.7633333 0.7633333 0.7633333 0.7633333
#>  [5,] 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7666667 0.7633333
#>  [6,] 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#>  [7,] 0.7633333 0.7633333 0.7633333 0.7666667 0.7666667 0.7633333 0.7633333
#>  [8,] 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7666667
#>  [9,] 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7666667 0.7633333
#> [10,] 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#> [11,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 1.0000000
#> [12,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [13,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 1.0000000
#> [14,] 0.9966667 0.9966667 0.9966667 1.0000000 1.0000000 0.9966667 0.9966667
#> [15,] 0.0000000 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [16,] 0.9966667 0.0000000 1.0000000 0.9966667 0.9966667 0.9966667 0.9966667
#> [17,] 0.9966667 1.0000000 0.0000000 0.9966667 0.9966667 0.9966667 0.9966667
#> [18,] 0.9966667 0.9966667 0.9966667 0.0000000 1.0000000 0.9966667 0.9966667
#> [19,] 0.9966667 0.9966667 0.9966667 1.0000000 0.0000000 0.9966667 0.9966667
#> [20,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.0000000 0.9966667
#> [21,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.0000000
#> [22,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 1.0000000
#> [23,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 1.0000000 0.9966667
#> [24,] 1.0000000 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [25,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 1.0000000
#> [26,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [27,] 1.0000000 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [28,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [29,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 1.0000000
#> [30,] 0.9966667 0.9966667 0.9966667 1.0000000 1.0000000 0.9966667 0.9966667
#> [31,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8400000 0.8366667
#> [32,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8400000 0.8366667
#> [33,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8400000 0.8366667
#> [34,] 0.8366667 0.8366667 0.8366667 0.8400000 0.8400000 0.8366667 0.8366667
#> [35,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8400000
#> [36,] 0.8400000 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [37,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [38,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [39,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [40,] 0.8366667 0.8400000 0.8400000 0.8366667 0.8366667 0.8366667 0.8366667
#> [41,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8400000 0.8366667
#> [42,] 0.8400000 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [43,] 0.8366667 0.8400000 0.8400000 0.8366667 0.8366667 0.8366667 0.8366667
#> [44,] 0.8366667 0.8366667 0.8366667 0.8400000 0.8400000 0.8366667 0.8366667
#> [45,] 0.8366667 0.8366667 0.8366667 0.8400000 0.8400000 0.8366667 0.8366667
#> [46,] 0.8366667 0.8400000 0.8400000 0.8366667 0.8366667 0.8366667 0.8366667
#> [47,] 0.8366667 0.8400000 0.8400000 0.8366667 0.8366667 0.8366667 0.8366667
#> [48,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [49,] 0.8366667 0.8400000 0.8400000 0.8366667 0.8366667 0.8366667 0.8366667
#> [50,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#>           [,22]     [,23]     [,24]     [,25]     [,26]     [,27]     [,28]
#>  [1,] 0.7666667 0.7633333 0.7633333 0.7666667 0.7633333 0.7633333 0.7633333
#>  [2,] 0.7633333 0.7633333 0.7666667 0.7633333 0.7633333 0.7666667 0.7633333
#>  [3,] 0.7633333 0.7600000 0.7600000 0.7633333 0.7600000 0.7600000 0.7600000
#>  [4,] 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#>  [5,] 0.7633333 0.7666667 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#>  [6,] 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#>  [7,] 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#>  [8,] 0.7666667 0.7633333 0.7633333 0.7666667 0.7633333 0.7633333 0.7633333
#>  [9,] 0.7633333 0.7666667 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#> [10,] 0.7633333 0.7633333 0.7633333 0.7633333 0.7666667 0.7633333 0.7666667
#> [11,] 1.0000000 0.9966667 0.9966667 1.0000000 0.9966667 0.9966667 0.9966667
#> [12,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [13,] 1.0000000 0.9966667 0.9966667 1.0000000 0.9966667 0.9966667 0.9966667
#> [14,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [15,] 0.9966667 0.9966667 1.0000000 0.9966667 0.9966667 1.0000000 0.9966667
#> [16,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [17,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [18,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [19,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [20,] 0.9966667 1.0000000 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [21,] 1.0000000 0.9966667 0.9966667 1.0000000 0.9966667 0.9966667 0.9966667
#> [22,] 0.0000000 0.9966667 0.9966667 1.0000000 0.9966667 0.9966667 0.9966667
#> [23,] 0.9966667 0.0000000 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [24,] 0.9966667 0.9966667 0.0000000 0.9966667 0.9966667 1.0000000 0.9966667
#> [25,] 1.0000000 0.9966667 0.9966667 0.0000000 0.9966667 0.9966667 0.9966667
#> [26,] 0.9966667 0.9966667 0.9966667 0.9966667 0.0000000 0.9966667 1.0000000
#> [27,] 0.9966667 0.9966667 1.0000000 0.9966667 0.9966667 0.0000000 0.9966667
#> [28,] 0.9966667 0.9966667 0.9966667 0.9966667 1.0000000 0.9966667 0.0000000
#> [29,] 1.0000000 0.9966667 0.9966667 1.0000000 0.9966667 0.9966667 0.9966667
#> [30,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [31,] 0.8366667 0.8400000 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [32,] 0.8366667 0.8400000 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [33,] 0.8366667 0.8400000 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [34,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [35,] 0.8400000 0.8366667 0.8366667 0.8400000 0.8366667 0.8366667 0.8366667
#> [36,] 0.8366667 0.8366667 0.8400000 0.8366667 0.8366667 0.8400000 0.8366667
#> [37,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [38,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [39,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [40,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [41,] 0.8366667 0.8400000 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [42,] 0.8366667 0.8366667 0.8400000 0.8366667 0.8366667 0.8400000 0.8366667
#> [43,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [44,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [45,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [46,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [47,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [48,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [49,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [50,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#>           [,29]     [,30]     [,31]     [,32]     [,33]     [,34]     [,35]
#>  [1,] 0.7666667 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7666667
#>  [2,] 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#>  [3,] 0.7633333 0.7600000 0.7600000 0.7600000 0.7600000 0.7600000 0.7633333
#>  [4,] 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#>  [5,] 0.7633333 0.7633333 0.7666667 0.7666667 0.7666667 0.7633333 0.7633333
#>  [6,] 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#>  [7,] 0.7633333 0.7666667 0.7633333 0.7633333 0.7633333 0.7666667 0.7633333
#>  [8,] 0.7666667 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7666667
#>  [9,] 0.7633333 0.7633333 0.7666667 0.7666667 0.7666667 0.7633333 0.7633333
#> [10,] 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#> [11,] 1.0000000 0.9966667 0.8366667 0.8366667 0.8366667 0.8366667 0.8400000
#> [12,] 0.9966667 0.9966667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [13,] 1.0000000 0.9966667 0.8366667 0.8366667 0.8366667 0.8366667 0.8400000
#> [14,] 0.9966667 1.0000000 0.8366667 0.8366667 0.8366667 0.8400000 0.8366667
#> [15,] 0.9966667 0.9966667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [16,] 0.9966667 0.9966667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [17,] 0.9966667 0.9966667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [18,] 0.9966667 1.0000000 0.8366667 0.8366667 0.8366667 0.8400000 0.8366667
#> [19,] 0.9966667 1.0000000 0.8366667 0.8366667 0.8366667 0.8400000 0.8366667
#> [20,] 0.9966667 0.9966667 0.8400000 0.8400000 0.8400000 0.8366667 0.8366667
#> [21,] 1.0000000 0.9966667 0.8366667 0.8366667 0.8366667 0.8366667 0.8400000
#> [22,] 1.0000000 0.9966667 0.8366667 0.8366667 0.8366667 0.8366667 0.8400000
#> [23,] 0.9966667 0.9966667 0.8400000 0.8400000 0.8400000 0.8366667 0.8366667
#> [24,] 0.9966667 0.9966667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [25,] 1.0000000 0.9966667 0.8366667 0.8366667 0.8366667 0.8366667 0.8400000
#> [26,] 0.9966667 0.9966667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [27,] 0.9966667 0.9966667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [28,] 0.9966667 0.9966667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [29,] 0.0000000 0.9966667 0.8366667 0.8366667 0.8366667 0.8366667 0.8400000
#> [30,] 0.9966667 0.0000000 0.8366667 0.8366667 0.8366667 0.8400000 0.8366667
#> [31,] 0.8366667 0.8366667 0.0000000 1.0000000 1.0000000 0.9966667 0.9966667
#> [32,] 0.8366667 0.8366667 1.0000000 0.0000000 1.0000000 0.9966667 0.9966667
#> [33,] 0.8366667 0.8366667 1.0000000 1.0000000 0.0000000 0.9966667 0.9966667
#> [34,] 0.8366667 0.8400000 0.9966667 0.9966667 0.9966667 0.0000000 0.9966667
#> [35,] 0.8400000 0.8366667 0.9966667 0.9966667 0.9966667 0.9966667 0.0000000
#> [36,] 0.8366667 0.8366667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [37,] 0.8366667 0.8366667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [38,] 0.8366667 0.8366667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [39,] 0.8366667 0.8366667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [40,] 0.8366667 0.8366667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [41,] 0.8366667 0.8366667 1.0000000 1.0000000 1.0000000 0.9966667 0.9966667
#> [42,] 0.8366667 0.8366667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [43,] 0.8366667 0.8366667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [44,] 0.8366667 0.8400000 0.9966667 0.9966667 0.9966667 1.0000000 0.9966667
#> [45,] 0.8366667 0.8400000 0.9966667 0.9966667 0.9966667 1.0000000 0.9966667
#> [46,] 0.8366667 0.8366667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [47,] 0.8366667 0.8366667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [48,] 0.8366667 0.8366667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [49,] 0.8366667 0.8366667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [50,] 0.8366667 0.8366667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#>           [,36]     [,37]     [,38]     [,39]     [,40]     [,41]     [,42]
#>  [1,] 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#>  [2,] 0.7666667 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7666667
#>  [3,] 0.7600000 0.7600000 0.7600000 0.7600000 0.7600000 0.7600000 0.7600000
#>  [4,] 0.7633333 0.7633333 0.7633333 0.7633333 0.7666667 0.7633333 0.7633333
#>  [5,] 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7666667 0.7633333
#>  [6,] 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#>  [7,] 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#>  [8,] 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#>  [9,] 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7666667 0.7633333
#> [10,] 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#> [11,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [12,] 0.8366667 0.8366667 0.8400000 0.8366667 0.8366667 0.8366667 0.8366667
#> [13,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [14,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [15,] 0.8400000 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8400000
#> [16,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8400000 0.8366667 0.8366667
#> [17,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8400000 0.8366667 0.8366667
#> [18,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [19,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [20,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8400000 0.8366667
#> [21,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [22,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [23,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8400000 0.8366667
#> [24,] 0.8400000 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8400000
#> [25,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [26,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [27,] 0.8400000 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8400000
#> [28,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [29,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [30,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [31,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 1.0000000 0.9966667
#> [32,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 1.0000000 0.9966667
#> [33,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 1.0000000 0.9966667
#> [34,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [35,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [36,] 0.0000000 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 1.0000000
#> [37,] 0.9966667 0.0000000 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [38,] 0.9966667 0.9966667 0.0000000 0.9966667 0.9966667 0.9966667 0.9966667
#> [39,] 0.9966667 0.9966667 0.9966667 0.0000000 0.9966667 0.9966667 0.9966667
#> [40,] 0.9966667 0.9966667 0.9966667 0.9966667 0.0000000 0.9966667 0.9966667
#> [41,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.0000000 0.9966667
#> [42,] 1.0000000 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.0000000
#> [43,] 0.9966667 0.9966667 0.9966667 0.9966667 1.0000000 0.9966667 0.9966667
#> [44,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [45,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [46,] 0.9966667 0.9966667 0.9966667 0.9966667 1.0000000 0.9966667 0.9966667
#> [47,] 0.9966667 0.9966667 0.9966667 0.9966667 1.0000000 0.9966667 0.9966667
#> [48,] 0.9966667 0.9966667 1.0000000 0.9966667 0.9966667 0.9966667 0.9966667
#> [49,] 0.9966667 0.9966667 0.9966667 0.9966667 1.0000000 0.9966667 0.9966667
#> [50,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#>           [,43]     [,44]     [,45]     [,46]     [,47]     [,48]     [,49]
#>  [1,] 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#>  [2,] 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#>  [3,] 0.7600000 0.7600000 0.7600000 0.7600000 0.7600000 0.7600000 0.7600000
#>  [4,] 0.7666667 0.7633333 0.7633333 0.7666667 0.7666667 0.7633333 0.7666667
#>  [5,] 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#>  [6,] 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#>  [7,] 0.7633333 0.7666667 0.7666667 0.7633333 0.7633333 0.7633333 0.7633333
#>  [8,] 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#>  [9,] 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#> [10,] 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333 0.7633333
#> [11,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [12,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8400000 0.8366667
#> [13,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [14,] 0.8366667 0.8400000 0.8400000 0.8366667 0.8366667 0.8366667 0.8366667
#> [15,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [16,] 0.8400000 0.8366667 0.8366667 0.8400000 0.8400000 0.8366667 0.8400000
#> [17,] 0.8400000 0.8366667 0.8366667 0.8400000 0.8400000 0.8366667 0.8400000
#> [18,] 0.8366667 0.8400000 0.8400000 0.8366667 0.8366667 0.8366667 0.8366667
#> [19,] 0.8366667 0.8400000 0.8400000 0.8366667 0.8366667 0.8366667 0.8366667
#> [20,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [21,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [22,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [23,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [24,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [25,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [26,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [27,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [28,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [29,] 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667 0.8366667
#> [30,] 0.8366667 0.8400000 0.8400000 0.8366667 0.8366667 0.8366667 0.8366667
#> [31,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [32,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [33,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [34,] 0.9966667 1.0000000 1.0000000 0.9966667 0.9966667 0.9966667 0.9966667
#> [35,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [36,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [37,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [38,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 1.0000000 0.9966667
#> [39,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [40,] 1.0000000 0.9966667 0.9966667 1.0000000 1.0000000 0.9966667 1.0000000
#> [41,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [42,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#> [43,] 0.0000000 0.9966667 0.9966667 1.0000000 1.0000000 0.9966667 1.0000000
#> [44,] 0.9966667 0.0000000 1.0000000 0.9966667 0.9966667 0.9966667 0.9966667
#> [45,] 0.9966667 1.0000000 0.0000000 0.9966667 0.9966667 0.9966667 0.9966667
#> [46,] 1.0000000 0.9966667 0.9966667 0.0000000 1.0000000 0.9966667 1.0000000
#> [47,] 1.0000000 0.9966667 0.9966667 1.0000000 0.0000000 0.9966667 1.0000000
#> [48,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.0000000 0.9966667
#> [49,] 1.0000000 0.9966667 0.9966667 1.0000000 1.0000000 0.9966667 0.0000000
#> [50,] 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667 0.9966667
#>           [,50]
#>  [1,] 0.7633333
#>  [2,] 0.7633333
#>  [3,] 0.7600000
#>  [4,] 0.7633333
#>  [5,] 0.7633333
#>  [6,] 0.7633333
#>  [7,] 0.7633333
#>  [8,] 0.7633333
#>  [9,] 0.7633333
#> [10,] 0.7633333
#> [11,] 0.8366667
#> [12,] 0.8366667
#> [13,] 0.8366667
#> [14,] 0.8366667
#> [15,] 0.8366667
#> [16,] 0.8366667
#> [17,] 0.8366667
#> [18,] 0.8366667
#> [19,] 0.8366667
#> [20,] 0.8366667
#> [21,] 0.8366667
#> [22,] 0.8366667
#> [23,] 0.8366667
#> [24,] 0.8366667
#> [25,] 0.8366667
#> [26,] 0.8366667
#> [27,] 0.8366667
#> [28,] 0.8366667
#> [29,] 0.8366667
#> [30,] 0.8366667
#> [31,] 0.9966667
#> [32,] 0.9966667
#> [33,] 0.9966667
#> [34,] 0.9966667
#> [35,] 0.9966667
#> [36,] 0.9966667
#> [37,] 0.9966667
#> [38,] 0.9966667
#> [39,] 0.9966667
#> [40,] 0.9966667
#> [41,] 0.9966667
#> [42,] 0.9966667
#> [43,] 0.9966667
#> [44,] 0.9966667
#> [45,] 0.9966667
#> [46,] 0.9966667
#> [47,] 0.9966667
#> [48,] 0.9966667
#> [49,] 0.9966667
#> [50,] 0.0000000
#> 
#> $sortind
#>  [1]  3  6 10  2  4  5  7  9  1  8 37 39 50 12 26 28 38 48 15 24 27 36 42 14 16
#> [26] 17 18 19 20 23 30 31 32 33 34 40 41 43 44 45 46 47 49 11 13 21 22 25 29 35