% -*- mode: noweb; noweb-default-code-mode: R-mode; -*- \documentclass[nojss]{jss} \usepackage{dsfont} \usepackage{bbm} \usepackage{amsfonts} \usepackage{wasysym} \usepackage{xypic} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% declarations for jss.cls %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \newcommand{\cd}{\cdot} %% just as usual \author{Robin K. S. Hankin\\University of Stirling\And Luke J. West\\National Oceanography Centre, Southampton} \title{Set Partitions in \proglang{R}} %\VignetteIndexEntry{Set partitions in R} %% for pretty printing and a nice hypersummary also set: \Plainauthor{Robin K. S. Hankin, Luke G. West} \Plaintitle{Set partitions in R} \Shorttitle{Set partitions in \proglang{R}} %% an abstract and keywords \Abstract{ To cite the set partitions functionality discussed in this document, use~\cite{hankin2007a}; to cite the \pkg{partitions} package in general, use~\cite{hankin2006}. This short paper introduces a code snippet in the form of an \proglang{R} function that enumerates all possible partitions of a finite set given the sizes of the equivalence classes. Three combinatorial problems are solved using the software: one from bioinformatics, one from scheduling theory, and one from forensic science. } \Keywords{Set Partitions, Forensic Science, Enumerative Combinatorics, \proglang{R}} \Plainkeywords{Set Partitions, Forensic Science, Enumerative Combinatorics, R} %% publication information %% NOTE: This needs to filled out ONLY IF THE PAPER WAS ACCEPTED. %% If it was not (yet) accepted, leave them commented. %% \Volume{13} %% \Issue{9} %% \Month{September} %% \Year{2004} %% \Submitdate{2004-09-29} %% \Acceptdate{2004-09-29} %% The address of (at least) one author should be given %% in the following format: \Address{ Robin K. S. Hankin\\ University of Stirling\\ Scotland\\ \email{hankin.robin@gmail.com}\\ {}\\ Luke J. West\\ National Oceanography Centre, Southampton\\ European Way\\ Southampton SO14 3ZH\\ United Kingdom\\ {\rule{10mm}{0mm}}\hfill\includegraphics[width=1in]{\Sexpr{system.file("help/figures/partitions.png",package="partitions")}} } %% It is also possible to add a telephone and fax number %% before the e-mail in the following format: %% Telephone: +43/1/31336-5053 %% Fax: +43/1/31336-734 %% for those who use Sweave please include the following line (with % symbols): %% need no \usepackage{Sweave.sty} %% end of declarations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \newsymbol\leqslant 1336 \SweaveOpts{echo=FALSE} \begin{document} <>= <>= require(partitions) @ \hfill\includegraphics[width=1in]{\Sexpr{system.file("help/figures/partitions.png",package="partitions")}} \section{Introduction} A {\em partition} of a set~$S=\{1,2,\ldots,n\}$ is a family of sets~$T_1,T_2,\ldots T_k$ satisfying \begin{enumerate} \item $T_i\cup T_j=\emptyset$ if $i\neq j$ \item $\bigcup_{i=1}^k T_k=S$ \item $T_i\neq\emptyset$ for $i=1,\ldots,k$. \end{enumerate} There are exactly fifteen ways to partition a set of four elements, shown in Table~\ref{partitions.4elements} using standard notation. Here~$(13)(2)(4)$ represents the same partition as~$(2)(4)(31)$: the induced equivalence relation is identical. Either one may be represented as \code{1213}, indicating that elements one and three are in the first equivalence class (``\code{1}''), element two in the second, and element four in the third. In Table~\ref{partitions.4elements}, two partitions appear in the same column if their `shape' (that is, the size distribution of the equivalence classes, formally defined in Section~\ref{algorithms} below) is the same. The shapes correspond to integer partitions in standard form~\citep{hankin2006}, and appear in standard order~\citep{andrews1998}. This suggests a natural method for enumerating the equivalence relations on a set of cardinality~$n$: enumerate the integer partitions of~$n$, and for each of these, enumerate the distinct equivalence relations of that shape. \begin{table}[h] \centering \begin{tabular}{ccccc} (1234) & (123)(4) & (12)(34) & (12)(3)(4) & (1)(2)(3)(4)\\ & (124)(3) & (13)(24) & (13)(2)(4) & \\ & (134)(2) & (14)(23) & (14)(2)(3) & \\ & (234)(1) & & (23)(1)(4) & \\ & & & (24)(1)(3) & \\ & & & (34)(1)(2) & \\ \end{tabular} \caption{The fifteen partitions of a set of four \label{partitions.4elements} elements} \end{table} This paper introduces software that enumerates\footnote{The transitive verb ``to enumerate'' can be used in two senses: firstly, ``to calculate the number of''; and secondly, ``to list each occurrence of, as if for the purpose of counting''. This paper adopts the second usage.} all set partitions of a set of a specified (finite) size. The new functionality presented here is given by \code{setparts()}, an \proglang{R} \citep{rcore2007} wrapper for \proglang{C++} code, currently part of the \pkg{partitions} package~\citep{hankin2006,hankin2007}, version 1.7-3, available under the GPL from CRAN, \url{http://www.r-project.org/}. Function \code{setparts()} takes a single argument \code{a}, coerced to integer mode. The default case is that of \code{a} being a vector of length greater than one, in which case it returns all partitions of a set of \code{sum(a)} elements into equivalence classes of the sizes of the elements of \code{a}; in the case of \code{a} being of length one, it returns all partitions of a set of \code{a} elements. If \code{a} is a matrix, it returns all equivalence classes with sizes of the columns of \code{a}. \section{Algorithms} \label{algorithms} In this paper, we present an algorithm that lists all partitions of a given ``shape''. Formally, the {\bf shape} $\lambda=(\lambda_1,\ldots\lambda_k)$ of a partition is given by the sizes of the~$T_i$; these are conventionally given in standard (i.e. non-increasing) order. Thus~$\lambda$ is an integer partition of~$n$ \citep{hankin2006}, written~$\lambda\vdash n$. The distinct equivalence relations may be enumerated, for a given integer partition, using a recursive algorithm that operates on partially filled partitions. Consider an integer partition~$\lambda= \left(\lambda_1,\ldots\lambda_k\right)$, also written~$\left(1^{f_1},\ldots,n^{f_n}\right)$; the latter notation indicates that~$f_i$ elements of~$\lambda$ are equal to~$i$. The algorithm used here fills the equivalence classes sequentially from left to right using a technique originally due to~\citet{killworth2007}. The partitions are ordered from largest to smallest; non-empty ties are broken lexicographically. This prevents the reporting of identical equivalence classes (to see this, consider the two identical equivalence relations~$(123)(45)(67)(8)$ and~$(123)(67)(45)(8)$: the second descriptor does not occur in the algorithm because the third equivalence class precedes the second). \begin{figure}[htbp] \begin{center} \includegraphics{algorithm} \caption{The\label{algorithm} algorithm used to enumerate set partitions of a given shape. Unfilled positions are indicated by a dot. Of the partitions at a given level, the children of one are shown. The initial empty configuration is filled from left to right from a `pool' of unused elements; equivalence classes are ordered by size (largest first) and ties are broken lexicographically at runtime to suppress reporting of identical equivalence relations in which two classes are transposed. Within a class, the elements of partially- and completely- filled equivalence classes are ordered in increasing order. A cross indicates an illegal configuration in which the partially filled class cannot be completed as insufficient elements remain in the pool; the algorithm prunes such configurations as soon as they are detected. Successful termination of the algorithm, viz. a filled partition, is indicated by a tick (``check mark'') at which point the filled equivalence class is reported, and the recursion bottoms out} \end{center} \end{figure} Thus, starting with an empty partition of a given size, one may apply the algorithm detailed in Figure~\ref{algorithm} recursively until each full partition is recorded. \subsection{Performance and complexity} The total number of set partitions corresponding to a particular integer partition~$\lambda$ is given by elementary combinatorial arguments as \begin{equation}\label{combinatoric_setpartition} \frac{n!}{\prod_{i=1}^{k}\lambda_i!\cdot\prod_{j=1}^n{f_j}!}. \end{equation} However, the number of integer partitions increases rapidly with~$n$; the exact value is given by the partition function~\code{P(n)} of the package~\citep{hankin2006}, but the asymptotic form given by~\citet{hardy1918}, viz. \begin{equation}\label{asymptotic_setpartition} P(n)=\frac{e^{\pi\sqrt{2n/3}}}{4n\sqrt{3}}\cdot\left(1+O(n^{-1})\right) \end{equation} shows that the rate of growth is rapid. The number of {\em set} partitions of a set of size~$n$ is given by the Bell numbers~$B(n)$~\citep{rota1964}, which grow much more rapidly than the partition function; an asymptotic form was given by~\cite{bender1974}: \begin{equation}\label{asymptoticbell} B(n) = t^{n-t}e^{t-1}\left(\log t\right)^{-1/2}\cdot\left(1+o(n)\right) \end{equation} where~$t=e^{W_0(n-1/2)}$. Here~$W_0(\cdot)$ denotes the principal branch of the Lambert W function~\citep{corless1996}, provided with the \pkg{gsl} package~\citep{hankin2006b}. Equation~\ref{asymptoticbell} is easily evaluated and can be used to give a feel for the magnitude of the Bell numbers. Note that convergence is quite slow; for example, $B(10)=115975$ (compare about~$1.48\times 10^3$ from Equation~\ref{asymptoticbell}), and~$B(100)\simeq 4.76\times 10^{115}$ (compare~$5.43\times 10^{115}$). \section{Applications} Many problems in science and industry involve the optimization of some desideratum over a finite number of options. Such optimization is frequently soluble within the discipline of computational combinatorics: simply enumerate the possible solutions and choose the best. One advantage of this method over heuristic optimization techniques such as tabu search \citep{glover1997} is that by enumerating all possible solutions, one is \emph{certain} to find the global optimum. In this section we present three examples of \code{setparts()} in use: one from bioinformatics, one from multiprocessor scheduling, and one from forensic science. \subsection{Bioinformatics} In studies of allelic segregation under disomic inheritance, one situation encountered is that of enumerating the possible arrangements of alleles on genetic loci\footnote{A \emph{chromosome} is a double helix of DNA, typically occurring in matched pairs in cell nuclei. A \emph{locus} is an identifiable position on a chromosome pair; a \emph{gene} is a distinct sequence of DNA base pairs occurring at a locus; an \emph{allele} is one of a number of distinguishable forms for a particular gene; the \emph{genome} is the complete set of genes present in the chromosomes; a \emph{monosomic locus} contains two copies of one, and a \emph{disomic locus} two different, alleles. Two alleles at the same locus are said to be \emph{colocated}. A monosomic (disomic) allele is one that occurs at a monosomic (disomic) locus in a given genome. A \emph{gamete} is a sequence of alleles chosen from the genome: monosomic and disomic alleles are chosen with probability~1 and~0.5 respectively.}. \cite{rodzen2002} consider the white sturgeon; there is much uncertainty surrounding the organization of this species's genome. Many alleles of interest may be identified in the genome of the white sturgeon. Sometimes these are known to be either monosomic or disomic, but it is not known which alleles are colocated, nor which alleles are monosomic and which are disomic. If researchers have identified, say, five alleles (conventionally labelled 1-5), and further that two pairs of these are colocated at two disomic loci and one is monosomic, then there are a number of possible arrangements for these alleles on the chromosome. In the language of the package, we wish to enumerate the partitions of a set of five elements (the alleles) into equivalence classes of size~2,2,1 (the loci). The five alleles are arranged into two from one disomic locus, two from another disomic locus, and one from a monosomic locus. Enumerating the possible arrangements is accomplished using the \code{setparts()} function: <>= setparts(c(2,2,1)) @ Each column comprises two~\code{1}s, two~\code{2}s, and one~\code{3}. Thus the first column, \code{x=12321}, say, indicates that alleles~1 and~5 [\code{which(x==1)}] occur at the first locus (which is disomic); alleles~2 and~4 [\code{which(x==2)}] at the second (also disomic); and allele~3 [\code{which(x==3)}] at the third locus, which is monosomic. This could be written~$(15)(24)(3)$; the \proglang{R} idiom would be <>= a <- c(1,2,3,2,1) split(seq_along(a),a) @ [use \code{listParts(c(2,2,1))} to see all possible equivalence classes]. Two examples of electrophoresis gels, corresponding to the partitions of the first and last columns respectively, are shown in figure~\ref{electrophoresis_gels}, which is analogous to Figure~1 of~\citet{rodzen2002}. \begin{figure}[htbp] \begin{center} \includegraphics{strand} \caption{Two electrophoresis gels \label{electrophoresis_gels} showing the possible bands that can occur with two different arrangements of alleles. In both gels, two loci are disomic, and one is monosomic. Left, gel corresponding to $(15)(24)(3)$, or \code{12321}; this would mean that alleles~1 and~5 are at one disomic locus, alleles~2 and~4 at another, and allele~3 at a monosomic locus. Thus band combinations~1-5 and~2-4 are absent. The gel on the right corresponds to~$(1)(24)(35)$, or \code{31212}} \end{center} \end{figure} \subsection{Scheduling theory} The ``augmented multiprocessing model'' is an abstract description of a certain class of problems that occurs frequently in the field of scheduling theory~\citep{garey1975}. We consider a simple variant of the augmented multiprocessing model (specifically, the ``three partition problem'', which is known to be NP-complete) and solve it using the software presented here. Consider a set of tasks~${\mathcal T}=\{T_1,\ldots,T_m\}$; task~$i$ requires time~$\tau_i$ for its completion in the sense that if it begins executing at time~$t$ it will complete at time~$t+\tau_i$. Each of these tasks is to be executed by one of~$n$ processors. It is desired to allocate the tasks among the processors so that all tasks are completed in the minimum possible time. Thus if there are~$m=4$ tasks with~$T_i$ requiring~$\tau_i=i$ time units (``seconds'') for completion, and~$n=2$ processors, it is clear that allocating tasks~1 and~4 to one processor, and tasks~2 and~3 to the second is optimal. This would correspond to the set partition~$(14)(23)$; the time required is~$1+4=2+3=5$ seconds. Other allocations such as~$(24)(13)$ would take longer ($2+4=6$ seconds). More complex examples are not so straightforward. Consider the case where~$m=9$ and again~$\tau_i=i$; now there are~$n=3$ processors. The possible task schedules may be enumerated using the \code{restrictedparts()} function\footnote{Recall that \code{restrictedparts(a,b)} enumerates the possible partitions of \code{a} into at most \code{b} parts.} of the \pkg{partitions} package: <>= options(width=63) m <- 9 n <- 3 @ <>= jj <- setparts(restrictedparts(m,n,include.zero=FALSE)) summary(jj) @ Each column corresponds to a task allocation schedule. The first column, for example, indicates that tasks 1,2,3,4,7,8,9 are allocated to processor~1, task~5 to processor~2, and task~6 to processor~3 and the set partition might be written~$(1234789)(5)(6)$. In the matrix, the maximum entry of each column is~3, showing that all three processors are always used (setting \code{include.zero} to \code{TRUE} would relax this requirement and allow some columns to correspond to allocations in which only one or two processors are allocated a nonzero number of tasks). The objective function to be minimized is just the time for the processor which finished last: <>= tau <- 1:9 slowest <- function(x) max(tapply(tau, x, sum)) @ The minimal objective function is clearly 45/3=15. Identifying allocations that attain this bound is straightforward: <>= time.taken <- apply(jj,2,slowest) minimal.time <- sum(tau)/n jj[,time.taken == minimal.time] @ Thus~\Sexpr{sum(time.taken==minimal.time)} allocations attain the maximal possible efficiency: a conclusion difficult to achieve without enumeration. These optima may be differentiated using some secondary characteristic, such as even distribution of the {\em number} of tasks among the processors. This would prescribe either the eighth or ninth columns, which both allocate three tasks to each processor. \subsection{Forensic science} In the field of forensic science, one situation sometimes encountered is the following. Crimes~$C_j$, where~$i=1,\ldots,n$ have been committed; it is desired to infer how many perpetrators are responsible and which crimes each perpetrator is responsible for. One prominent example would be multiple unsolved homicides~\citep{roberts2003}. Evidence at each of the~$n$ crime scenes is available in the form of~$m$ Bernoulli random variables that indicate various forensic characteristics at the crime scene (examples might include presence or absence of a weapon, forced entry, etc). The evidence may be organized in a~$m$-by-$n$ matrix~$E$ of zeroes and ones indicating whether forensic characteristic~$i$ was true at crime scene~$j$, where~$1\leqslant i\leqslant m$ and~$1\leqslant j\leqslant n$. A given row of~$E$ thus corresponds to a particular type of evidence; if perpetrator~$k$ committed crime~$j$, we would have~$E_{ij}=1$ with probability~$p_{ik}$, specific to perpetrator~$k$. For each evidence~$i$, it is assumed that each perpetrator has a~$p_{ik}$ chosen at random from a beta distribution with parameters~$\alpha_i,\beta_i$ for~$1\leqslant i\leqslant m$, and that all perpetrators have the same~$\alpha_i,\beta_i$; these will be known (or at least estimated) from previous studies of similar crimes. An example is shown below. <>= E <- matrix(c(0,0,0,1,1,0,0,0,1,1,1,1,1,0,0), 5,3) dimnames(E) <- list( evidence=paste("E",1:5,sep=""), crime=paste("C",1:3,sep="") ) @ <>= E @ Thus in this example there are~$n=3$ crimes and~$m=5$ forensic characteristics. Observe that~\code{C1} and~\code{C2} have identical evidence, and that~\code{C3} has contrary evidence. This would suggest that there are exactly two perpetrators: one responsible for \code{C1} and \code{C2}, and one for \code{C3}, corresponding to partition~$(12)(3)$; the hypothesis that this is the case is written~$H_{[1,1,2]}$. Considering a single type of evidence, and a single perpetrator, the probability of observing~$a$ successes and~$b$ failures is \[ \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} \int_{p=0}^1 \left[ \frac{(a+b)!}{a!b!} p^{a}(1-p)^{b} \right] p^{\alpha-1}(1-p)^{\beta-1}\,dp= \frac{(a+b)!}{a!b!}\cdot \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\cdot \frac{\Gamma(a+\alpha)\Gamma(b+\beta)}{\Gamma(a+b+\alpha+\beta)}. \] Here the probability~$p$ is integrated over its range, weighted by its prior (the posterior would also be a beta distribution, with parameters~$\alpha+a$ and~$\beta+b$). Now consider a set partition~$\wp=T_1,\ldots,T_{r}$ of crimes; identify a perpetrator with each equivalence class of crimes and observe that a single perpetrator has~$m$ specific values of~$p_i$, $1\leqslant i\leqslant m$ and acts independently of other perpetrators. The likelihood of~$\wp$ is then \[ {\mathcal L}(\wp)= C\left(\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\right)^r \cdot\prod_{k=1}^{r}\prod_{i=1}^{m} \frac{\Gamma(a_{ik}+\alpha_i)\Gamma(b_{ik}+\beta_i)}{\Gamma(a_{ik}+b_{ik}+\alpha_i+\beta_i)} \] where~$a_{ik}=\sum_{j\in T_k}E_{ij}$ and~$b_{ik}=\sum_{j\in T_k}1-E_{ij}$ are the total successes and failures for evidence~$i$ due to perpetrator~$k$ (under~$H_\wp$). Here~$C$ denotes an arbitrary multiplicative constant. It is interesting to note that the case $\alpha=\beta=1$ is not uninformative in this context, even though this induces a uniform prior distribution on~$p$. To see this, consider the first line of the five by three matrix above [\code{0,0,1}]. The likelihood function for this row is on five partitions: \begin{eqnarray*} {\mathcal L}(H_{[1,1,1]}) &=& \frac{\Gamma(\alpha+1)\Gamma(\beta+2)}{\Gamma(\alpha+\beta+3)}\cdot \left(\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\right)^1\\ {\mathcal L}(H_{[1,2,1]}) &=& \frac{\Gamma(\alpha+1)\Gamma(\beta+1)}{\Gamma(\alpha+\beta+2)}\cdot \frac{\Gamma(\alpha)\Gamma(\beta+1)}{\Gamma(\alpha+\beta+1)}\cdot \left(\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\right)^2\\ {\mathcal L}(H_{[1,1,2]}) &=& \frac{\Gamma(\alpha)\Gamma(\beta+2)}{\Gamma(\alpha+\beta+2)}\cdot \frac{\Gamma(\alpha+1)\Gamma(\beta)}{\Gamma(\alpha+\beta+1)}\cdot \left(\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\right)^2\\ {\mathcal L}(H_{[2,1,1]}) &=& \frac{\Gamma(\alpha+1)\Gamma(\beta+1)}{\Gamma(\alpha+\beta+2)}\cdot \frac{\Gamma(\alpha)\Gamma(\beta+1)}{\Gamma(\alpha+\beta+1)}\cdot \left(\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\right)^2\\ {\mathcal L}(H_{[1,2,3]}) &=& \frac{\Gamma(\alpha)\Gamma(\beta+1)}{\Gamma(\alpha+\beta+1)}\cdot \frac{\Gamma(\alpha)\Gamma(\beta+1)}{\Gamma(\alpha+\beta+1)}\cdot \frac{\Gamma(\alpha+1)\Gamma(\beta)}{\Gamma(\alpha+\beta+1)}\cdot \left(\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\right)^3. \end{eqnarray*} Thus if~$\alpha=\beta=1$ the likelihood function is \begin{eqnarray*} {\mathcal L}(H_{[1,1,1]}) &=& C/12\\ {\mathcal L}(H_{[1,2,1]}) &=& C/12\\ {\mathcal L}(H_{[1,1,2]}) &=& C/6\\ {\mathcal L}(H_{[2,1,1]}) &=& C/12\\ {\mathcal L}(H_{[1,2,3]}) &=& C/8 \end{eqnarray*} so the highest likelihood is indeed that of~$H_{[1,1,2]}$. It is also instructive to assume~$\alpha=\beta$ and consider the limit as~$\alpha\longrightarrow 0$. This would correspond to perpetrators having either a very small~$p$, or~$p$ close to~1, intermediate values being rare. To first order in~$\alpha$: \begin{eqnarray*} {\mathcal L}(H_{[1,1,1]}) &=& C\cdot\alpha/4\\ {\mathcal L}(H_{[1,2,1]}) &=& C\cdot\alpha/4\\ {\mathcal L}(H_{[1,1,2]}) &=& C\cdot(1-\alpha)/4\\ {\mathcal L}(H_{[2,1,1]}) &=& C\cdot\alpha/4\\ {\mathcal L}(H_{[1,2,3]}) &=& C\cdot 1/8 \end{eqnarray*} so the maximum likelihood partition would again be~$H_{[1,1,2]}$; but one would not be able to reject the hypothesis that there are three distinct perpetrators (that is, $H_{[1,2,3]}$) at the two units of support level. As a final example, consider the case where~$\alpha=\beta$ and~$\alpha\longrightarrow\infty$. This corresponds to each perpetrator leaving positive evidence with probability~$\frac{1}{2}$. The five likelihoods all tend to the same limit. \subsubsection{Generalization to arbitrary factors} These ideas have a natural generalization to arbitrary factors. The beta distribution becomes a Dirichlet distribution and the posterior probability becomes \[ \frac{\left(\sum a_l\right)!}{\prod a_l!} \frac{\Gamma\left(\sum\alpha_l\right)}{\prod\Gamma\left(\alpha_l\right)} \frac{\prod\Gamma\left(a_l+\alpha_l\right)}{\Gamma\left(\sum a_l+\sum\alpha_l\right)} \] where outcome~$l$ is observed with frequency~$a_l$ and the corresponding prior distributions are Dirichlet with parameters~$\alpha_i=\left(\alpha_{i1},\ldots, \alpha_{iu}\right)$. The likelihood is then \begin{equation}\label{like.mult} {\mathcal L}(\wp)= C\prod_{i=1}^{n} \left[ \left( \frac{\Gamma\left(\sum_l\alpha_{il}\right)}{\prod_l\Gamma\left(\alpha_{il}\right)} \right)^r \prod_{k=1}^{r} %\cdot\frac{\left(\sum_{l}a_{ikl}\right)!}{\prod_la_{ikl}!} \frac{\prod_l\Gamma\left(a_{ikl} + \alpha_{il}\right)}{\Gamma\left( \sum_l\left(a_{ikl}+\alpha_{il}\right)\right)} \right] \end{equation} where~$a_{ikl}$ is the frequency of outcome~$l$ for evidence~$i$ and perpetrator~$k$, under~$H_\wp$. \paragraph{Example} <>= a <- structure( c(2, 2, 4, 4, 3, 1, 2, 1, 1, 2, 2, 1, 4, 2, 1, 2, 1, 4, 2, 4, 2, 1, 4, 4, 1, 2, 1, 4, 2, 1, 1, 2, 1, 1, 2), .Dim = c(5L, 7L), .Dimnames = list( evidence=paste("E",1:5,sep=""), crime=paste("C",1:7,sep="") ) ) @ Consider a case where seven crimes have been committed; there are five forensic characteristics. The evidence matrix is as follows: <>= a @ Each column corresponds to a crime and each row corresponds to a forensic characteristic. Rows~1 and~2 are Bernoulli: each is either~1 or~2. Rows~3-5 are Dirichlet with four possible outcomes. In this case, it is assumed that the prior is uniform, viz~$\alpha_1=\alpha_2=1$ for rows~1 and~2 and~$\alpha_1=\alpha_2=\alpha_3=\alpha_4=1$ for rows~3-5. <>= genbeta <- function(x){gamma(sum(x))/prod(gamma(x))} lp <- function(evidence, alpha, partition){ r <- length(unique(partition)) evidence <- as.factor(evidence) levels(evidence) <- seq_along(alpha) out <- 1 for(k in unique(partition)){ #k is a perp thisperp <- partition==k evidence.thisperp <- evidence[thisperp] no.of.crimes.thisperp <- sum(thisperp) out <- out/genbeta(alpha+table(evidence.thisperp)) } return(out*genbeta(alpha)^r) } sp <- setparts(7) l1 <- apply(sp,2, function(x){lp(evidence=a[1,],alpha=rep(1,2),partition=x)}) l2 <- apply(sp,2, function(x){lp(evidence=a[2,],alpha=rep(1,2),partition=x)}) l3 <- apply(sp,2, function(x){lp(evidence=a[3,],alpha=rep(1,4),partition=x)}) l4 <- apply(sp,2, function(x){lp(evidence=a[4,],alpha=rep(1,4),partition=x)}) l5 <- apply(sp,2, function(x){lp(evidence=a[5,],alpha=rep(1,4),partition=x)}) likelihood <- l1*l2*l3*l4*l5 likelihood <- likelihood/max(likelihood) support <- log(likelihood) @ \begin{figure}[htbp] \begin{center} <>= par(mfcol=c(1,2)) plot(likelihood,ylab="likelihood",xlab="hypothesis") abline(h = exp(-2)) plot(support,ylab="support",xlab="hypothesis") abline(h = -2) @ \caption{The\label{support.function} likelihood (left) and support (right) for each of the~877 partitions on a set of seven elements. The horizontal line shows two units of support less than the maximum} \end{center} \end{figure} Figure~\ref{support.function} shows the likelihood and support for each of the distinct set partitions on a set of six elements: variable \code{support} is from Equation~\ref{like.mult}. The maximum likelihood estimator is given by <>= sp <- setparts(7) <>= sp[,which.max(support)] @ which could be written~$(3456)(27)(1)$. This suggests that crimes~3-6 were committed by the same perpetrator, crimes~2 and~7 by another perpetrator, and crime~1 by a third perpetrator. Note that the maximum likelihood hypothesis is unique. The likelihood function may be used to identify those hypotheses that cannot be rejected on the grounds of insufficient support. There are~9 partitions that have fewer than two units of support less than the maximum. These are: <>= index <- order(support,decreasing=TRUE) sp <- sp[,index] support <- support[index] dimnames(sp) <- list( crime = paste("C",1:7, sep=""), partition = paste("H",1:ncol(sp),sep="") ) <>= sp[, support > -2] @ and the supports for these hypotheses are <>= support[support > -2] @ respectively. Interpretation of this result is not straightforward. One may assert that crimes~\code{C3} and~\code{C6} were committed by the same perpetrator, and that certain pairs of crimes (specifically \code{C2}-\code{C3}, \code{C2}-\code{C4}, \code{C2}-\code{C5}, \code{C2}-\code{C6}, \code{C3}-\code{C7}, \code{C4}-\code{C7}, \code{C5}-\code{C7}, and~\code{C6}-\code{C7}) were committed by different perpetrators: for any hypothesis that does not satisfy these criteria, one could gain at least two units of support by adopting instead the best supported hypothesis~$H_{[3,2,1,1,1,1,2]}$, corresponding to~$(3456)(27)(1)$. \section{Note: a gotcha} Note carefully that \code{setparts(c(2,1,1))} does \emph{not} enumerate the ways of placing four numbered balls in three labelled boxes of capacities 2,1,1. This is because there are two boxes of capacity 1, and swapping the balls between these boxes gives a different box partition but the same set partition (because sets are unordered). <>= setparts(c(2,1,1)) @ Note the absence of a column reading \code{1 1 3 2}. This would correspond to placing balls 1 and 2 in box 1 (of size 2), ball 3 in box 3 (of size 1) and ball 4 in box 2 (also of size 1). The missing column is because the two boxes of size 1 are indistinguishable, so \code{1 1 3 2} is the same \emph{set} partition as \code{1 1 2 3}, as in each case balls 3 and 4 are placed in ``a box of size 1''. If you want to enumerate the ways of choosing two workers, a secretary and a chair from four people (who are conveniently numbered 1,2,3,4), use \code{multinomial()}: <>= multinomial(c(worker=2,secretary=1,chair=1)) @ Above, see how the rows are named for their equivalence class. The first two columns differ only in the identity of the secretary and the chair (and would be identical set partitions as given by \code{setparts(c(2,1,1))}). The reason I mention this is to avoid the following gotcha (which looks plausible but is incorrect): <>= v <- c(worker=2,secretary=1,chair=1) a <- apply(setparts(v),2,order) rownames(a) <- rep(names(v),v) as.partition(a) @ The above is an incomplete enumeration. We see 1,4,2,3 [which would correspond to set partition 1,2,3,1] but not 1,4,3,2 [which would correspond to set partition 1,3,2,1]. These two differ in that the secretary and chair have swapped roles. \section{Conclusions} Enumeration of set partitions is required in several branches of computational combinatorics. The software discussed in this code snippet enumerates set partitions for finite sets, under a variety of restrictions on the sizes of the induced equivalence classes and was used to solve three problems drawn from diverse applications in computational statistics. \subsection*{Acknowledgement} We would like to acknowledge the many stimulating and insightful comments made on the \proglang{R}-help list while preparing this software. \bibliography{partitions} \end{document}