Corbelled Domes in Two and Three Dimensions: The Treasury of Atreus

Before the development of the true dome, many ancient cultures used the technique of corbelling to roof spaces. Recently, a series of related statistical models have been proposed in the literature for explaining how corbelled domes might have been constructed. The most sophisticated of these models is based on a piecewise linear structure, with an unknown number of changepoints, to guide the building process. This model is analyzed by the reversible jump Markov Chain Monte Carlo (MCMC) technique. All models considered to date have been two‐dimensional, that is, they have taken a single cross section through the dome; even when more extensive data, in the form of measurements on multiple slices through the dome, have been available, these have been averaged together for the purposes of analysis. In this paper, we extend the two‐dimensional analysis to a three‐dimensional analysis, that takes full advantage of the data collected by the archaeologists and of the rotational symmetries inherent in the structure.We also explore ways of graphically presenting the results from a complex, reversible jump MCMC implementation, in order to check convergence, good mixing, and appropriate exploration of the (high dimensional and varying dimension) parameter space. The model and the graphical techniques are demonstrated on the Treasury of Atreus in Mycenae, Greece, one of the finest extant examples of the corbelling method.


Introduction
Corbelling is an ancient technique, used before the development of the true dome, for spanning or roofing spaces. Examples of this construction method have been found in a number of prehistoric societies, ranging from tombs in Iberia and the British Isles, to buildings in Italy, to passages and chambers of the Great Pyramid of Kheops in Giza, Egypt. Most of these utilize corbelling to span relatively small distances, 2 to 3 meters at most, with a conservative slope.
One of the most impressive uses of the technique, which involves laying stones on top of each other, with successive layers being slightly offset (think of using rectangular Lego pieces to build an arch; see Figure 1), comes from Mycenae in Greece, where archaeologists have found corbelled structures spanning 8 to 14 meters. As noted by Cavanagh & Laxton (1981), only the invention of the true dome enabled builders to bridge larger distances without internal support.
Due to the prevalence of corbelling in the ancient world, and the fact that many of the structures which were built with this spanning method are still standing today, thousands of years later, Figure 1. Shape of the inner face of an idealized corbelled structure. Reprinted from Buck, Litton & Stephens (1993) with permission from the Royal Statistical Society.
archaeologists have become interested in studying this technique more closely. In particular, interest has focused on finding models that adequately describe the shape of corbelled domes, and perhaps give insight into how these structures were built in practice. Much has been written about this problem in the last twenty years, by both archaeologists and statisticians. Especially interesting from the statistical point of view is that the data lend themselves to a variety of surprisingly advanced analysis techniques, such as multiple changepoint detection and reversible jump Markov Chain Monte Carlo. In this paper, we discuss in detail the archaeological, mathematical and statistical developments that surround the analysis of corbelled domes. We also extend the current existing models to take advantage of the three dimensional nature of the structures, which to date has been ignored. The outline of the rest of the paper is as follows. Section 2 traces the development of the model, from a deterministic least squares line fit, to a multiple changepoint model cast in the Bayesian setting and analyzed using the most advanced MCMC techniques. These models are based on measurements from multiple cross-sections of the dome, that are averaged together to get a single, two-dimensional "slice". Section 3 introduces an extension of the existing models to include the third dimension of the structure. In other words, the measurements along the various cross-sections are no longer averaged together. Rotational symmetries are imposed on the individual slices, however, to model accurately the constraints that would have been placed on the builders. This three dimensional analysis is demonstrated, in some detail, on the Treasury of Atreus. Section 4 gives a discussion of the results and general conclusions.

Classical Analysis
In the first mathematical study, Cavanagh & Laxton (1981) considered the structural mechanics of the Mycenaean tholos tombs. Tholoi, which are also found extensively in Crete, are tombs constructed using the corbelling technique described above. By examining the structure of these tombs and the factors necessary for the vault to remain standing, they came up with the following, deterministic, model: where depth is the depth below the top of the tomb, radius is the radius of the tomb at the measured depth, and « and ¬ are constants to be determined for each tomb individually. Taking logarithms of each side of the equation, a simple linear equation results, and it is possible to find values for the two coefficients using least squares. Cavanagh and Laxton examined five tombs that varied in location, age and building style, and found that all five had a value of ¬ of around 2/3. Values of « varied, as this parameter depends on the size (specifically, the height) of each tomb. One of the tombs, the Treasury of Atreus, is almost twice the size of the other tombs in this study, and it had « ¾ , compared to « ¾ for the other four. Their model, while providing a starting point for the analysis of the Mycenaean tombs, was lacking in two important aspects, addressed in Cavanagh, Laxton & Litton (1985). First, many of the Mycenaean tholoi were built partly below ground. As described by Cavanagh, Laxton & Litton (1985), the floor of such tombs was made by cutting a cylindrical cavity out of the bedrock, or from the edge of a hillside. Thus, the Mycenaean tholoi were often either sunk into a hillside or consisted of a subterranean burial chamber with a corbelled roof; by way of contrast, the earlier Minoan circular tombs on the island of Crete were built above ground (see also http projects dartmouth edu history bronze age for historical details and comparisons). The walls in the lower portion of the tomb may have been constructed according to a different plan than the upper part. For instance, they may have been straighter, because of the constraints of the rock-cut "cylinder" or in order to accommodate the doorway. Instead of one value of « and one value of ¬ sufficing to describe the entire tomb, there might be a point at which the building regime, and the parameters guiding construction, change. The changepoint, if it exists, would correspond to where the dome and the wall meet. The first simple model considered by Cavanagh and Laxton was based only on measurements from the dome, and not the entire tomb. Thus, whereas their simple form was the basic one required for stability of the dome, it did not necessarily describe adequately the Mycenaean tombs of mainland Greece. Second, it might not be possible to measure from the apex of the dome, that is, at depth ¼. There are two reasons for this: the top of the tomb might have collapsed, or the corbelling wasn't originally continued to the top of the arch, instead a large slab was placed down on the stones, to close off the dome. Under these two modifications, and introducing a stochastic element to account for errors, the proposed model became Here, r i is the radius measured at depth d i , is the location of the changepoint, AE is an offset parameter that accounts for the fact that measurements might not be available all the way to the apex of the dome (some domes, for instance, were capped by a large stone), and the¯i are assumed to be normally distributed, with mean 0 and variance ¾ . Physical considerations require that ¬ ½ ¬ ¾ ¼. At depth , the two lines are required to intersect, that is, « ½ º AE» ¬ ½ « ¾ º AE» ¬ ¾ , so that there are only six free parameters in the problem, instead of seven. In particular, ¬ ¾ can be expressed as ¬ ½ ÐÓ º« ½ « ¾ » ÐÓ º AE» .
The authors used maximum likelihood techniques developed by Hinkley (1969Hinkley ( , 1971 for the analysis of changepoint problems, to compare the parameters of Greek tholoi and Sardinian nuraghi (Bronze Age stone towers which served as defensive structures). Based on the parameter estimates, it was possible to distinguish between the tholoi and the nuraghi. For example, for the Sardinian structures, no changepoints were apparent. It should be noted that the nuraghi, unlike the tholoi, are completely above ground, so one of the reasons for proposing a changepoint model doesn't hold, and this seemed to be confirmed by the data. For the tholoi under consideration in this study, AE was taken to be zero (that is, known and fixed), whereas for the nuraghi, it was allowed to vary, as an unknown parameter of the model. The method was also able to pick out one of the Mycenaean tombs, that at Karditsa, as different from the others in its group, since it appeared to have no changepoint. As pointed out by the authors, this tomb was unusual in two ways. One of the distinguishing features of the Karditsa tholos is that it is set entirely within a knoll. Cavanagh et al. speculated that the builders of the tomb were able to take advantage of this fact and adapted the construction regime accordingly. In addition, the tomb at Karditsa had been a source of some dispute among archaeologists, as to whether it belonged to the Bronze or Iron Age. The authors suggested that the statistical analysis could provide another, independent, criterion for dating this structure.

The Bayesian Perspective
An important step forward was taken by Buck, Litton & Stephens (1993), who used methods developed in Stephens (1994), for the analysis of changepoint problems from a Bayesian point of view. The authors applied MCMC methods to a Minoan tholos tomb at Stylos, Crete. Their base model was as in previous work. The conditions on the model were also as described above. With the assumption of normal errors, the likelihood is easily specified, although it should be noted that because of the parameter constraints and the piecewise nature of the model, sampling from the posterior is not simple, regardless of the prior specification. Buck and colleagues chose non-informative (vague) priors for ÐÓ « ½ , ÐÓ « ¾ , ¬ ½ , ¾ and ; the prior for AE was uniform between 0 and 0.34, to account for the thickness of the capping stone on this tomb, which was 0.34 meters. Under this specification, the Markov Chain Monte Carlo sampling scheme can be carried out using Gibbs sampling, which draws successive observations from the full conditional distributions, that is, the posterior distribution of each parameter given all of the others (Casella & George, 1992).
The three parameters ÐÓ « ½ , ÐÓ « ¾ and ¬ ½ have conditional posterior distributions in the normal family, while ¾ is conditionally inverse gamma. These are easily obtained. On the other hand, the full conditional densities for AE and are both of non-standard form, hence are sampled from using rejection sampling (see Buck, Litton &Stephens, 1993 andStephens, 1994 for details).
A number of advantages of the Bayesian approach are apparent for this problem. Instead of just point estimates of the parameters of interest, entire posterior distributions are obtained. This gives archaeologists more flexibility in studying and interpreting the results of the analysis. The classical analysis did not incorporate any uncertainty about AE into the inference process, whereas this was an outcome of the Bayesian analysis by its very nature. Furthermore, information about the main parameter of interest, the location of the changepoint, was refined. The 95% confidence interval reported for arising from the classical analysis (1.64-1.84), was apparently too short; the 95% highest posterior density interval from the Bayesian analysis was 1.48-2.44. As pointed out by Buck et al., this latter interval just includes the depth at which the lintel is below the capping stone of the tomb, 2.44 meters. This finding supported a supposition by archaeologists that the placing of the lintel would cause a break in the process of corbelling, with different structural forces coming in to play above and below that point. Using the classical approach, there was no support for this intuition, as the changepoint was placed too high.
This argument, however, may put too much weight on exactly where the top end of a 95% credible interval falls. The number 0.95 has nothing intrinsic to do with the problem. Its use is merely traditional; it is "treat [ed] with the same superstitious awe that is usually reserved for the number 13" (Raiffa & Schlaifer, 2000, p. vii).
Regardless of whether a single changepoint model points to a change at the level of the lintel, these considerations lead to the possibility of models with two or more changepoints, perhaps one at the lintel and one at the ground level at the time of construction. That is, are there features of a corbelled structure that might impose an adjustment in the parameters guiding construction? One could consider expanding the model described above to account for more than one changepoint. However, the appropriate number of changepoints is not known, and so it would be desirable to include this as a parameter in a new model. While seemingly a minor modification of the existing framework, in fact letting the number of building segments vary introduces nontrivial complications. The problem is that, as the number of changepoints varies across models, interpretation of the model parameters also changes. In addition, the practical issue of how to estimate the true dimension of the model, by sampling over the space of models with differing numbers of changepoints, is not an easy one to resolve.
Fortunately, at around the time we became interested in this question, a solution was proposed in the Bayesian literature, in the form of Reversible Jump MCMC (Green, 1995;Richardson & Green, 1997). The reversible jump technique designs a Markov chain that jumps between models of different dimensions, while still preserving the conditions needed to eventually reach a stationary distribution. Implementation of the method can be complicated and the details tend to be unique to each application. Denison, Mallick & Smith (1998) derived a reversible jump algorithm for dealing with multiple changepoints, however they made certain simplifying assumptions which are not appropriate for the problem we are working on.
More recently, Fan & Brooks (2000) developed methodology specifically for the analysis of corbelled structures. They extended the basic model to allow for two generalizations. First, AE was not restricted to be the same in each linear piece. Second, at lower depths, the wall of the tomb might be nearly vertical, that is, governed by a single parameter. An example of a candidate two changepoint model is The number of segments was allowed to increase to a maximum of four, by permitting as many as three changepoints. The reversible jump algorithm implemented by Fan and Brooks also explored the possibility that AE was in fact the same for all segments, and the possibility that there was no constant segment at the bottom of the tomb. In addition to reanalyzing the data from the tomb at Stylos in Crete, Fan and Brooks compared the results of their reversible jump MCMC algorithm for another Minoan tholos to those for a Mycenaean tholos and a Sardinian nuraghe. They found that for each, the method assigned slightly different orderings of the posterior probabilities, that is, indicated different underlying models as having the highest posterior probabilities. The estimated values of the ¬ parameters in particular were similar enough in all three of the buildings considered in this comparison, that the authors concluded that this information alone was not enough to distinguish among structures of differing origin. It is worth pointing out, however, that interpretation of the ¬ parameters in the comparison is complicated by the fact that different models were fit to each of the structures-a two changepoint model with a constant segment, a one changepoint model with no constant segment, and a two changepoint model with no constant segment for the Minoan tholos, the Mycenaean tholos and the nuraghe, respectively. Since the models are different, the interpretation of ¬ in each case is also different, making an overall conclusion such as that drawn by Fan & Brooks (2000) tentative at best.

Three Dimensional Data
Corbelled domes are three dimensional structures, and many, like the tholos tombs, were built rotationally symmetric, or nearly so, around a center line (though later settlement, whilst still leaving the structure stable, could distort slightly the original symmetry). For this reason, the analyses reviewed above are in fact based on averages of measurements taken on several slices starting at the apex. The model discussed below is the first, to our knowledge, to consider each measurement slice separately. Since a lintel and doorway intercept one of the slices, simple averaging, as was performed previously, is neither appropriate nor sufficient. The slice that is intersected by the lintel and the doorway has a different number of measurements than the other slices, and some of the measurements were taken at different places, as explained below. Note that under the assumption that all measurement slices contain complete data taken at the same depths, a simple sufficiency argument shows that no information is lost by averaging. However when one or more slices has missing data or has measurements taken at different locations, averaging makes little sense and results in the loss of information. From the practical point of view, furthermore, the model we develop here, which can easily handle missing data, gives the archaeologist more flexibility, since difficulties in collecting complete data along all cross-sections is not an impediment to statistical analysis.

Background and Data
The Treasury of Atreus, also known as the Tomb of Agamemnon, is a tholos tomb in Mycenae, Greece. It is considered to be the finest tomb to have survived and dates to about 1300 BCE. It is far from typical of the majority of tholos tombs because of the unusual masonry used in its construction. An interesting feature of the Treasury is that the lintel is very long, large and heavy, with a relieving triangle above. The role of this relieving triangle is to concentrate the weight of the masonry on to the jambs rather than on the lintel. This helps prevent the lintel from cracking and breaking.
A detailed architectural survey was carried out and accurate plans made. From the plans, it is possible to determine the radius of the tholos at a given depth. We have data on three separate cross-sections of the Treasury of Atreus, as shown in Figure 2. Cross-section G-H cuts through the entrance of the tomb, and thus the relieving triangle as well. The center of the tomb is denoted O. For each cross-section, the two radii were measured at regular depths from the apex of the dome. Note that for slice OH, there are no measurements for the relieving triangle, nor are there any for below the lintel. Radii OC, OD, OE, OF, OG and OH have 46, 45, 46, 45, 45 and 21 measurements respectively. The data are given in Table 1.

The Statistical Model
We used the base model proposed by Fan & Brooks (2000), with the exception that for the Treasury of Atreus, the AE parameter is identically zero, since the apex of the structure is known (very unusually, the form of the last few centimetres of the dome were carved out of the capstone of the Treasury of Atreus, giving a very precise location of where the builders placed the apex of the vault). More specifically, the data are in the form of d i y i ¡ , where d i with i ½ n j denotes the depth below the apex of the tomb, and y i ¨y i j j ½ © are corresponding radius measures at that depth, for the five slices with complete data, and one with incomplete data. The radii y i j are modelled as normal centered at r i , with unknown variance ¾ k , and k indicates the segment in the piecewise linear model to which y i j belongs. An individual observation therefore belongs to both a slice, that is, a particular cross-section of the tomb, and to a segment of the piecewise linear fit. As before, d i and r i enjoy the following basic relationship: Notice that we now need to make a distinction between the measured radii, y i j , and r i , the radius at depth d i which acts also as the theoretical mean of the observations at depth d i . Since we are in a Bayesian framework, these r i are themselves taken to be random quantities, drawn from a distribution with parameters of its own. Previous analysis did not need this hierarchy, because only a single radius measurement was used at each depth.
Given that we believe that the tomb is formed segmentwise with changepoints γ, equation (1) can likewise be extended to allow for multiple changepoints, for instance, two changepoints with a constant segment at the bottom: with continuity constraints In the presence of missing data, we denote obs mis , for "observed" and "missing" data respectively. obs , the data, are modelled via the likelihood as Normal centered at the corresponding Ö with variance ¾ . mis are treated as parameters, and samples from its predictive distribution p y i j ¡ can be easily obtained: From a sampling point of view, we first draw samples from the posterior distribution of Ö ºh» ¾ºh» ¡ , where ¨ obs ºh ½» mis © , h is the current step in the Gibbs setting. Next we use the drawn samples of Ö and ¾ to produce ºh» mis , where y i j ∼ N r i ¾ k ¡ with the corresponding i (location), j (slice), and k (segment).
The second level of the hierarchical model is defined in the following way: We chose priors to be both conjugate (hence convenient from a computational point of view, and consistent with previous analysis, such as that of Fan and Brooks) and reasonably diffuse.

Results
The data were analyzed using reversible jump MCMC (Green, 1995;Richardson & Green, 1997). The implementation of the reversible jump method followed Fan & Brooks (2000). By using the reversible jump approach, we avoid the need to explicitly calculate Bayes factors to choose among the models; the reversible jump algorithm provides information about the size of the model, since one of the parameters is the number of changepoints. We thus obtain the posterior probabilities of each model directly from the chain. It is worth emphasizing that the choice of approaching the problem by the use of reversible jump MCMC is, in fact, a choice, and not a requirement imposed by the model or the data. As pointed out by DiCiccio et al. (1997), a variety of techniques have been proposed over the years for dealing with the question of calculating Bayes factors and posterior probabilities of competing models. These methods are based on asymptotic approximations (for instance, the Laplace approximation, Tierney & Kadane, 1986), simulation (marginal likelihood, Chib (1995), Chib & Jeliazkov (2001); bridge sampling, Meng & Wong (1996), Gelman & Meng (1998)), and a combination of both (various modifications on existing methods proposed by DiCiccio et al., 1997). We explored some of these other options for our model, but found the reversible jump approach to be the most straightforward to implement and to interpret.
To analyze the data, we used a program written in C, which we modified from a program used by Fan and Brooks for the single slice case. The modified program samples all of the parameters in accordance to the formulae detailed above. It is worth noting that even with the Fan and Brooks program as a guide, writing a program to handle the three dimensional data was not a trivial undertaking. A serious amount of effort was required to institute the necessary changes and to debug the resultant collection of programs. Interested readers may acquire our program by request from the authors. We selected five very different starting values, as can be seen in Figure 3. For each model, we ran fifteen million iterations. Every five thousand iterations, we recorded which of the models the chain was visiting at that iteration and its log(posterior) value. In order to assess the stability and convergence of the chain, we looked at the traces of the average log(posterior) for each of the five runs, the cumulative probability of spending time in what turned out to be a much-visited model for each of the five runs, and the joint plot of these two measures. These are Figures 4, 5 and 6, respectively. Inspection of these Figures reveals that from all starting points the chains converged to similar places.
Since the results from the five runs, starting at five different locations, all converged to more or less the same place, we then ran a much longer chain, of sixty million iterations, discarded the first ten million as burn in, and kept the results of every five thousandth sample. We ran the chain as long as we did, and thinned it so severely, because of the high amount of correlation among the parameters, induced by the constraints in the model. See Buck et al. (1993) for discussion of this point in the context of the much simpler Bayesian model based on a single slice, with a single changepoint. The long burn in, the long chain, and the extreme thinning are all indicative of a very conservative stance, and are designed to guarantee, as far as possible, both that convergence has been attained and that posterior quantities are derived from independent samples. Table 2 shows the proportion of time spent, out of the resulting ten thousand iterations, in each of the models. As can be seen from inspection of the Table, by far the most time was spent in models with 1 or 2 changepoints. Models without a constant segment at the end were overwhelmingly preferred. Only two models were visited with any frequency-the one changepoint model with no constant segment at the bottom, and the two changepoint model with no constant segment. All other models were rarely visited. Table 3 gives the posterior summaries for the two most visited models.
The reversible jump analysis points to the model with a single changepoint and no constant segment at the bottom, as being by far the preferred model. This model matches somewhat with a structural interpretation. As described above, the OH slice has no measurements for the relieving triangle, nor below the lintel. The relieving triangle thus starts at depth 4.7 meters. The lintel starts at depth 8.1 meters. The posterior mean of ½ for the preferred model is 7.8 meters, near the start of the lintel. The other popular model places the two changepoints lower (at an average of 7 meters Table 3 Results for the two most visited models-a one changepoint model and a two changepoint model; in both cases, the bottom segment is not a constant.

changepoints 1 changepoint Mean
± HPDI Mean ± HPDI and 10.5 meters respectively) than the structural points, but the latter are included in the 95% HPD intervals.
In addition, the estimated radius measures for the missing data along the OH slice are very close to the observed values from the five complete slices. Figures 7 and 8 show the estimated variances at each observation, for the two popular models, along with 95% pointwise confidence intervals. The vertical lines give the locations of the changepoints. As can be seen from the Figures, the variance increases at greater depths, that is, near the bottom of the tomb. Our supposition is that this is due to the fact that near the top, where the angles must be steeper in order to create the characteristic shape, and to close off the space, there is less opportunity for variability in placing the stones. At the bottom, this is not the case.
One of the challenges in understanding the results from a complex model-simulation combination such as the one presented here, is to find clear ways of representing the outcome. Elsewhere (Lazar & Kadane, 2002) we have suggested the use of a statistical movie, to look at the posterior distribution output from high parameter problems, where simple trace plots might not suffice. Figure 9 shows another way of visualizing the simulation output for our particular model. In the Figure, we have plotted a typical model sampled from the posterior distribution of «s, ¬s and s, for the model with one changepoint and no constant segment at the bottom. Also on this plot is the full posterior distribution for the parameter, that is, the location of the changepoint in this model. The two together give information on where the changepoints for these models tend to be located, how variable these locations are, and what sort of fits might result from this particular placement of changepoints. While the plot doesn't have the dynamic nature of the movie, it does provide much more information than either a simple movie snapshot or a trace plot on their own could. A similar plot for the two changepoint model with no constant appears in Figure 10.
Inspection of the Figures reveals some interesting features of the posterior distributions. Much of the effort is spent in placing changepoints near the bottom, indicating that more adjustments were required during the construction of the lower layers of the tomb than of the higher ones. However, the model also considers the possibility of placing a changepoint near the top of the tomb. This can be seen by the long tail in the posterior for the location of the first changepoint in Figure 10.
Finally, it is possible to plot the output in three dimensions, as in Figure 11, which shows the data, along with the fitted model from a random iteration in the chain. A movie of this three dimensional plot shows how the fitted shape of the tomb changes as the simulation progresses. A clip from such a movie, which lasts approximately one minute, can be viewed at http www stat cmu edu nlazar movie mpeg.

Discussion
The analysis of corbelled domes, as described here, demonstrates the way in which advances in statistical thinking and methodology can interact with a scientific discipline, to further understanding of some aspect of that discipline. Indeed, in the example developed in this article, it is clear that statistical methods needed to adapt to the nature of the data and the problem-stretching the initial, deterministic models to allow for a Bayesian changepoint analysis, for example, required the development of Stephens (1994). Going in the other direction, Fan & Brooks (2000) utilized the new statistical technique of reversible jump MCMC to introduce the possibility of multiple changepoints, but of an unknown number. Here, the methodology existed, although the application was new. Finally, incorporating the true, three-dimensional nature of the structures, and the availability of such data, albeit partial in some cases by necessity, required an additional generalization of Fan & Brooks (2000). Indeed, as we have argued, the presence of missing data in the slice that passed through the doorway of the tomb meant that simply averaging the data from different cross-sections, as had been previously done, would have resulted in a loss of information. In other words, what has unfolded here is the ideal interplay between statistics and the demands of the scientific discipline, applications driving new theory, which is turned back again to the applied problem.
The process will not end at this point-nor should it. Corbelled structures are very common across the Mediterranean region, and elsewhere. Archaeologists are still uncertain as to whether the method developed independently in different eras and regions, or was spread via word of mouth and apprenticeships. While statistics alone is not likely to answer this question, the use of more complicated spatio-temporal models that incorporate the age and geography of the structures might help to narrow the possibilities down and provide further research directions for the archaeologists. This is the logical next step in the analysis of these data.