text:
\section*{RINEX Version 3.00 Library Enhancements}

The substantial changes between RINEX Version 3.00 (herein note R3)
and RINEX Version 2 (R2) resulted in the consideration of several
upgrade paths for implementation. The major design considerations for
discriminating advantages and disadvantages between choices were as
follows.
%
\begin{itemize}
  \item Full and straightforward backwards compatibility.
  \item Ease of testing.
  \item Minimal modifications required across the remainder of the
toolkit to integrate the new features.
\end{itemize}
%

Ultimately the choice was made to create an independent set of R3
classes with distinct names, simultaneously branching the GPSTk within
its Subversion repository.  This choice allows backwards compatibility
with R2, and allows for separate known versions for comparison testing
of R3 versus R2.
%
\begin{figure*}
   \centering
   \includegraphics[width=5.5in,bb=79 313 1369 1107]{UML0709.eps}
   \caption{Unified Modeling Language (UML) diagram showing inheritance relationships among RINEX Version 3.00 classes.}
   \label{fig:rinex3}
\end{figure*}
%
Figure~\ref{fig:rinex3} depicts the class structure introduced to
support R3 capabilities.

%-------------------------------------------------------------------------------
\subsection*{Time System}

The GPSTk originally assumed that all times are given in the GPS time
system.  Since this is not true in a multi-GNSS world, a Time System
(TS) data structure was added to keep track of what time frame epochs
are actually reported in. The R3 standard specifies that epoch data
for a given GNSS be given in that GNSS's native time system.  As a
result, the GPSTk has been designed to automatically encode the TS
data structure along with the observations.

%-------------------------------------------------------------------------------
\subsection*{Coordinate System}

In the same vein, each GNSS uses its own realization of the
International Terrestrial Reference Frame (ITRF) as its coordinate
system.  For GPS this is WGS84(G1150), while for GLONASS it is PZ-90.
Galileo will use the Galileo Terrestrial Reference System (GTRF).
Japan's QZSS will use JGS, which differs from the GTRF by less than
2~cm.  The disagreement between different systems varies in general,
but the agreed-upon realization goal is 2-3~cm in the near future.

We have added support for the standard reference frames to the GPSTk,
as well as the ability to transform data between reference frames via
a Helmert transformation~\cite{Hofmann2008}, defined by
%
\begin{equation*}
{\bf X}_{A} = {\bf a} + \mu {\bf R} {\bf X}_{B}.
\end{equation*}
%
in which $\bf X$ denotes a position vector as a column matrix, $A$ and $B$ denote reference frames, $\mu$ is a scale factor and ${\bf a}$ is an offset or bias.

%-------------------------------------------------------------------------------

\subsection*{Data Storage and Access}

Observation, navigation, and meteorological data have to be stored for
access (there are only those three data types in R3).  This data
storage should be efficient.  The current implementation of R2
handling uses a ``map of maps'' (MoM) approach \cite{Stroustrup2000}.
In this approach, the actual data is stored in instances of a data
element class, which become the values in a map keyed by epoch (time).
Instances of that map are then stored as values in another map, this
one keyed by satellite vehicle ID.  The advantage of this approach is
that this creates a poor-man's ``database'' that is easily accessible
by the most identifiable tag, time.  Other possible approaches include
filling matrices, valarrays, linked lists or possibly other data
types.

For R3, we choose to stay with the MoM design based on an assessment that
the current implementation of data storage in R2 is
acceptable in usability and speed, as well as being quite flexible. However, there are
some revisions.  For
example, each GNSS has its own MoM, belonging to its own storage
class.  A new class specific to reading R3 files contains an instance
of each possible GNSS MoM, and inserts new data into the appropriate
MoM as it is read in.  Look-ups can then be targeted to a specific
GNSS, or scanned over all MoMs, depending on the user's requirements.

%-------------------------------------------------------------------------------

\subsection*{Navigation Data}

R3 expanded the navigation broadcast message (``Nav'') data type to
include general GNSS Nav messages.  These formats can be quite
different. For example, Galileo's broadcast ephemerides format is
similar to GPS \cite{Hofmann2008}. However, GLONASS, the Space-based augmentation
system (SBAS), and QZSS, use tabular ephemerides \cite{glonassicd50}.  A major
restructuring of the GPSTk ephemeris processing classes was required in order to
be able to handle the R3 standard, though it currently includes only
GLONASS.

The data storage class was more cleanly split into two
branches.  One is the parent class for all tabular-type GNSSs
(GLONASS, QZSS), while the other is for those which broadcast orbital
element data (GPS, Galileo).  It should be noted that precise
ephemerides in SP3 files are presented in tabular format.  Thus, the
SP3 classes remain a subclass of \gpstkclass{TabularEphemerisStore}.  We remark here only that separate
naming for the R3 classes allowed us to keep the R2 classes in place,
and with invisible backwards compatibility.  That is, R2 users will
not have to change their existing applications to use the new GPSTk
library.

%-------------------------------------------------------------------------------

\subsection*{Observation Data}

Observation (``Obs'') data refers to any measurement of (pseudo)range,
phase, Doppler or other signal quantity made with a GNSS receiver.  R3
Obs data files have changed significantly, providing support for
additional data, e.g. moving receivers or antenna details, but there
were subtle, significant structural changes.  Chief among these is Obs
type code list and format.  Obs files may now include multiple GNSS
data in a file, which requires more detailed header information on
what data will is present for each GNSS.  There are new data
structures for satellite ID.  Epochs are presented differently.

The new R3 Obs classes in the GPSTk had to be significantly modified
to handle these changes.  Some of the data is stored in maps (and even
maps of maps) to simplify access over different GNSSs.  The extent of
changes in the R3 standard will require Obs data users to
significantly change the way they write their applications.  It has
always been up to the user to store Obs data once read in; we continue
to not provide any inherent Obs data storage classes (e.g. MoMs) in
this upgrade.

max_stars_repo_name: mfkiwl/ICE


text:
%% Domenico Giusti
%% <domenico.giusti@uni-tuebingen.de>
%%
%% Paläoanthropologie
%% Senckenberg Center for Human Evolution and Paleoecology
%% Eberhard Karls Universität Tübingen
%% Rümelinstr. 23, 72070 Tübingen, Germany
%%
%% ERC PaGE Project
%% Palaeoanthropology at the Gates of Europe
%% Human Evolution in the Southern Balkans
%%
%% LaTeX manuscript -- English {elsarticle}

\documentclass[review,authoryear]{elsarticle} %%[preprint/review/5p]

\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
\usepackage[english]{babel}

\usepackage{lineno}
%\modulolinenumbers[5]

\usepackage{textcomp}

\usepackage{amsmath}

\usepackage{hyperref}

\journal{Journal of Archaeological Science: Reports}

\begin{document}

\begin{frontmatter}

  \title{The need for a taphonomic perspective in spatial analysis: Formation processes at the Early Pleistocene site of Pirro Nord (P13), Apricena, Italy}

  \author[tue]{D.~Giusti\corref{cor1}}
  \ead{domenico.giusti@uni-tuebingen.de}
  \cortext[cor1]{Corresponding author}

  \author[fe]{M.~Arzarello}

  \address[tue]{Paläoanthropologie, Senckenberg Center for Human Evolution and Paleoecology, Eberhard Karls Universität Tübingen, Rümelinstr. 23, 72070 Tübingen, Germany}
  \address[fe]{Dipartimento di Studi Umanistici, Università degli Studi di Ferrara, C.so Ercole I d'Este 32, 44100 Ferrara, Italy}

  \begin{abstract}
Ever since their percolation from neighbour disciplines, archaeology has employed spatial statistics to unravel, at different scales, past human behaviors from scatters of material culture. However, in the interpretation of the archaeological record, particular attention must be given to disturbance factors that operate in post-depositional processes. In this paper, we answer the need for a specific taphonomic perspective in spatial analysis by applying point pattern analysis of taphonomic alterations on the faunal and lithic assemblages from the Early Pleistocene site of Pirro~Nord~13, Italy. The site, biochronologically dated between 1.3 and 1.6 Ma BP, provides evidence for an early hominin presence in Europe. The archaeological and paleontological deposits occur as filling of a karst structure that is currently exposed. We investigated the distribution of the archaeological and paleontological assemblages, as well as the distribution of identified taphonomic features, in order to evaluate the degree and reliability of the spatial association of the lithic artifacts with the faunal remains. Our results contribute to the interpretation of the diagenetic history of Pirro~Nord~13 and support the stratigraphic integrity of the site.
  \end{abstract}

  \begin{keyword}
    Spatial analysis \sep Point pattern analysis \sep Site formation processes \sep Taphonomy \sep Early Pleistocene \sep Lower Palaeolithic \sep Pirro~Nord~13
  \end{keyword}

\end{frontmatter}

\linenumbers

\section{Introduction}

Studies of site formation processes and spatial analyses have long recognized the role of post-depositional factors in affecting the integrity of archaeological assemblages \citep{Hodder1976,Petraglia1987,Schick1984,Schick1986,Schiffer1972,Schiffer1983,Schiffer1987,Wood1978}. More recently, a number of scholars have stressed the importance of establishing the degree of disturbance to archaeological deposits to fully comprehend the archaeological record \citep{Dibble1997,Djindjian1999,Texier2000}.

Besides geoarchaeological techniques, several archaeological and paleontological methods are widely applied to characterize the processes involved in the formation of an archaeological site and to assess any post-depositional `background noise`. Taphonomy moves from its original definition \citep{Efremov1940} to a wider conceptual framework, targeting vertebrate assemblages, as well as taphonomic entities produced by human behaviour \citep{Dominguez-Rodrigo2011}. Moreover and often in joint effort, from different spatial perspectives, fabric analysis \citep{Benito-Calvo2011,Bernatchez2010,Bertran1997,Bertran1995,Dominguez-Rodrigo2014,Lenoble2004,McPherron2005,Torre2013a}; refitting analysis \citep{Lopez-Ortega2011,Sisk2008,Villa1982}; vertical \citep{Anderson2008} and size distribution analysis \citep{Bertran2006,Bertran2012,Petraglia1994} offer meaningful contributions in the unraveling of site formation and modification processes.

The importance of spatial statistics in the interpretation of archaeological sites has long been recognized \citep{Hodder1976,Whallon1974}. However, studies of spatial patterning mostly focus on the behaviour of past populations, assuming that scatters of material culture (if not disturbed) are reflections of prehistoric activities. Moreover, distribution maps still rely mainly on visual examinations and subjective interpretations \citep{Bevan2013a}. On the other hand, quantitative methods, adopted from neighbor disciplines since the early 1970s (see \cite{Hodder1976,Orton1982}; and references therein), continue to promote new impulses to archaeological spatial analyses and allow for the characterization of spatial patterns by adopting a more formal, inductive approach. Recent studies \citep{Bevan2006,Bevan2009,Bevan2013c,Bevan2013a,Bevan2013,Crema2015,Crema2010,Crema2013,Eve2014,Orton2004}, even acknowledging post-depositional effects or research biases, have continued to adopt at different scales (from intra-site to regional scales) improvements in spatial statistics to unravel past human behaviors from scatters of material culture. Yet, only a relatively limited number of scholars have applied spatial statistics to site formation and modification processes analysis \citep{Carrer2015,Dominguez-Rodrigo2014b,Dominguez-Rodrigo2014c}.

In this paper, we adopt a taphonomic perspective to spatial point pattern analysis of the lithic and faunal assemblages from the Early Pleistocene site of Pirro~Nord~13, Italy \citep{Arzarello2007,Arzarello2009,Arzarello2012,Arzarello2010}.

The site (P13) provides important contributions to the ongoing debate about the first hominin occurrence in Europe \citep{Carbonell2008,Crochet2009,Despriee2006,Despriee2009,Despriee2010,Lumley1988,Pares2006,Toro-Moyano2011,Toro-Moyano2009,Toro-Moyano2013}. A `Mode~1` lithic assemblage has been identified in stratigraphic association with late Villafranchian/early Biharian paleontological remains. Furthermore, the presence of the Arvicolinae species \emph{Allophaiomys~ruffoi} correlated to the \emph{Mymomis~savini-Mymomis~pusillus} biozone, allows for a  biochronologically refined age between 1.3 and 1.6 Ma, making P13 one of the most ancient localities with human evidence currently known in Western Europe \citep{Lopez-Garcia2015}.

The paleontological and archaeological remains are preserved inside a complex karst system, exposed and partially destroyed by mining activities of a Mesozoic limestone quarry. The fissure P13 is a vertical fracture located at the stratigraphic boundary between the Mesozoic limestone and the Pleistocene calcarenite formation. The deposit of the fissure is, at the time of writing, more than 4 m thick. Four Sedimentary Units (SUs) have been distinguished on lithological basis. From the top to the bottom of the section, units A to D are characterized by sediments of clayey-sand of increasing thickness (Fig.~\ref{fig:1}). Unit A includes few coarse gravels and a very low number of paleontological and archaeological remains. Unit B contains more gravels, while an abrupt increase in the number and dimension of clasts and large blocks of Pleistocene calcarenite is evident within units C and D. These last units show poor size sorting of angular and sub-rounded gravels, probably correlating to a low degree of reworking that took place during a short interval of time. We also record a significant increase in the number of fossils and artifacts.

\begin{figure}
  \centering
  \includegraphics[width=1\textwidth]{../artwork/Fig1.pdf}
  \caption{Geographical location of the Pirro Nord site. Picture of the fissure P13 inside the Cave Dell'Erba quarry and view of the excavated area (2013), with marked bases of the sedimentary units.}
  \label{fig:1}
\end{figure}

As a residual component of a wider karst system, it is worthwhile to assess the degree of any potential post-depositional reworking of the archaeological and paleontological remains and to evaluate the stratigraphic integrity of the site.

The main goal of our study is to use a taphonomic perspective in spatial data analysis, in order to evaluate degree and reliability of the spatial association of the lithic artifacts with the faunal remains that were used for the biochronological dating of the site.

By applying point pattern analysis of the spatial distribution of the lithic and faunal assemblages, we aim to
\begin{enumerate}
  \item investigate the processes involved in the formation of the Pirro~Nord (P13) deposit.
\end{enumerate}
A positive spatial association of the two types of find would support the assumption, base on field observations, that the deposition of the archaeological and paleontological materials occurred simultaneously, as a result of subsequent mass wasting events.

With the application of point pattern analysis to identified taphonomic features on the lithic and faunal assemblages, our ultimate objective is to
\begin{enumerate}
  \setcounter{enumi}{1}
  \item evaluate the degree of post-depositional disturbance of the site.
\end{enumerate}
Indeed, reworking and re-deposition processes could put in stratigraphic contact materials from diverse provenience. The identification of taphonomic spatial patterns allows us to model the spatial processes that produced them and thus propose a reconstruction of the agents involved in the formation and modification of the deposit.

\section{Background}

With the authors' permission, we integrate in our study unpublished \citep{Bagnus2011} and published \citep{Arzarello2012,Arzarello2015} data from previous taphonomic studies. A brief report is presented here.

\subsection{Taphonomy of macrovertebrate fossils}

Taphonomic analysis \citep{Bagnus2011} on macrovertebrate fossils evaluated biostratinomic and diagenic processes and grouped faunal remains into different sub-categories: three main taphorecords \citep[TRs, \emph{sensu}][]{Fernandez-Lopez1987} are defined according to different stages of bone surface modifications by physical and chemical agents (Tab.~\ref{tab:1}). Grouping was based mainly on weathering \citep{Behrensmeyer1978,Diez1999,Kos2003,Torres2003}, abrasion \citep{Behrensmeyer1991} and oxidation \citep{Hill1982,Lopez-Gonzalez2006,White1976,White2009}, because these alterations prevail and are widespread across all the sedimentary units.

Based on macroscopic observations of these main taphonomic features, fossils from TR2 and TR3 are interpreted as re-deposited fossils: displaced bones along the sedimentary surface before burial; whereas fossils from TR1 are considered re-elaborated \citep[\emph{sensu}][]{Fernandez-Lopez1991,Fernandez-Lopez2007,Fernandez-Lopez2011}. The higher degree of abrasion and the presence in the latter sub-group of multiple generations of oxides, non-uniformly distributed on the fossil, are explained with repeated exhumations and dislocations of previously buried elements \citep{Lopez-Gonzalez2006}.

Therefore, a hypothetical model of site formation processes has been proposed: animals died close to the karst sinkhole and the action of heavy rains transported sediments and partially articulated carcasses into the fissure. The rapid burial of fossils is confirmed by the general low degree of weathering. Karst erosional processes are responsible for the very large percentage of fractured fossils, as a result of the collapse of rock blocks from the vault. The TR1 group of fossils points to internal water-flows, reworking and transportation of already fossilized bones. Finally, manganese oxides that give the external widespread black color to all the fossils, stones and part of the lithic artifacts are products of the phreatic water fluctuation.

Although the taphonomic analysis definitely improved the interpretation of the P13 fossiliferous deposit, the interactions between bones and karst water flow have not been studied in relation to the spatial distribution and orientations of the skeletal elements.

Taking into account the inherent spatial properties of taphonomic processes, we assume that taphogenic products \citep[\emph{sensu}][]{Fernandez-Lopez2000} in space are not mutually independent and that entities which are close to each other, are likely to have followed the same genesis.

Thus, in order to tackle our second objective, we analyze the spatial distribution of Fe-Mn oxides on the fossils, since the cause of their formation may derive from the action of circulating waters. Three ordinal degrees of oxidation (low, medium and high) are recognized, based on its aspect, intensity and extension. We assume that spatial aggregation of heavily-coated faunal remains (and consequently segregation from non-oxidized ones) is an indication of interactions with karst water flow.

\begin{table}
  \caption{Contingency table of taphorecords \citep[reproduced from][]{Bagnus2011}.}
  \label{tab:1}
  \vspace{0.1in}
  \centering
  \begin{tabular}{l c c c r}
    \hline
    SU & TR1 & TR2 & TR3 & Total by SU \\
    \hline
    A & 10 & 30 & 45 & 85 \\
    B & 26 & 86 & 114 & 226 \\
    C & 34 & 69 & 179 & 282 \\
    Total by TR & 70 & 185 & 338 & 593 \\
    \hline
  \end{tabular}
\end{table}

\subsection{Taphonomy of lithic artifacts}

The degree of natural alterations (thermal, tribological and chemical) of the lithic artifact surface, as a result of contact with the sediments, is a valuable index of integrity of the depositional context and it can usefully support spatial analysis in reconstructing both the past environmental conditions and the site formation processes \citep{Burroni2002}.

According to a recent review of preliminary technological analyses \citep{Arzarello2012,Arzarello2015}, the lithic assemblage shows a general good state of preservation. If we consider the degree of patination as a good indicator of the intensity, and not necessary of the duration, of chemical processes to which the deposit has been subjected \citep{Burroni2002}, artifacts undergo non-homogeneous interactions with chemical agents. Besides fresh artifacts, many of the specimens ($35\%$) bear Fe-Mn coatings (Fig.~\ref{fig:2}a). Iron-manganese, as well as white superficial patina ($5\%$), seems to equally affect artifacts of different flint raw materials, more readily on those with a porous structure (Fig.~\ref{fig:2}b).

Macroscopic observations of tribological features on the assemblage reveal mint to sharp, not rounded, artifact ridges and edges. Post-depositional fractures affect $20\%$ of the lithic material (Fig.~\ref{fig:2}c).

No refittings were found, as it is reasonable to expect for materials in a secondary context.

As particle size distribution of lithic assemblages has great implications in interpreting site formation processes \citep{Bertran2012}, systematic screen-washing of sediments was carried out in order to guarantee recovery of lithic debris, even though a very low percentage of small-size specimens has been noted. This result can be initially explained either as a function of the mode of knapping, which did not produce a lot of debris, or is more likely due to natural post-depositional processes (winnowing of low energy agents), prior to final burial, possibly outside the karst fissure. Moreover, the dimensional analysis of the complete lithic assemblage (Fig.~\ref{fig:2}d) does not show sorting effects.

We analyze the spatial distribution of taphonomic features on the lithic assemblage, considering that various natural mechanisms, disturbing the spatial arrangement of artifacts and sediments, will produce distinctive combinations of wear features on the surfaces of lithic artifacts \citep{Burroni2002}.

As for the faunal assemblage, we focus the analysis on the distribution of Fe-Mn patinae. Three ordinal degrees of patination (absent, spotted and covering) are recognized, based on its presence and extension. In order to evaluate the impact of post-depositional processes at the site, we conduct independent and comparative taphonomic spatial analyses with the fossil remains.

\begin{figure}
  \centering
  \includegraphics[width=1\textwidth]{../artwork/Fig2.pdf}
  \caption{Frequency of patinae on the lithic assemblage, grouped by SUs (a) and their distribution on raw materials (b); frequency of fractures on the lithic assemblage, grouped by SUs (c); scatterplot of artifact dimensions (d).}
  \label{fig:2}
\end{figure}

\section{Spatial data collection and sampling}

Since 2007, systematic field investigations of the P13 fissure have been carried out by the University of Ferrara (in collaboration with the Universities of Torino and Roma Sapienza, until 2010).

From the first excavation season, a grid of $1$ square meter units has been set. Since 2010, the three-dimensional coordinates of the finds are recorded with a Total Station, which replaced the use of a water level. Orientation (dip and strike) of coordinated faunal remains (length $\geq2$~cm), geological clasts (length $\geq5$~cm) and all the lithic artifacts is estimated with a $45$~degree of accuracy, which is not precise enough for detailed fabric analysis.

In order to avoid possible sampling issues in spatial data analysis due to the variation in the recording methods, we select subsets of the lithic and faunal collection, excluding SUs A and B, because they have been excavated prior the use of the Total Station.

Focusing on SUs C and D, we scale the windows of analysis according to the extension of excavated areas for each SU, excluding the presence of the large blocks of rock. We reduce in this way the impact of the Modifiable Area Unit Problem (MAUP) in point pattern analysis \citep{Openshaw1996}, especially insidious in this study due to the particular geological setting of the site. The analyzed areas of SUs C and D are respectively $4.34~m^{2}$ and $5.82~m^{2}$.

During 6 years of excavations, more than $1600$ of $2152$ macrovertebrate fossils have been spatially recorded: $471$ from SU C and $916$ from SU D. However, \citet{Bagnus2011} conducted taphonomical analysis on fossils recovered during the 2007 to 2010 field seasons and only $593$ of these are classified in one of the three taphorecords (Tab.~\ref{tab:1}). Our sample includes $135$ coordinated elements of the $282$ analyzed fossils from SU C. From the total number of $366$ lithic artifacts collected until the 2014 field season, $147$ have been recorded with three-dimensional coordinates. Our sample includes $34$ lithics from SU C and $84$ from SU D. From the micro-mammal assemblage, we include in this study only the \emph{Allophaiomys ruffoi} species. Of the $53$ arvicoline teeth collected from the screen-washed sediments, $49$ have secure provenance attribution from SU B ($n=2$), C ($n=14$) and D ($n=33$) \citep{Lopez-Garcia2015}. However, the \emph{A.~ruffoi} point pattern does not represent the exact distribution of the remains. Indeed, we randomly displaced ($r=0.5$) each point indicating the provenience of the sieved sediment.

\section{Vertical distribution}

The vertical distribution of finds is a key factor in the analysis of site formation processes. Many processes can be well approximated by a ‘nearly’ normal distribution. However, testing the appropriateness of this assumption is an essential step in spatial data analysis. Strongly right skewed distribution would occur in case of a non-uniform vertical distribution of finds; thus requiring the analysis to acknowledge the covariate effect of gravity in the observed spatial pattern.

The vertical distribution of finds within SU C is globally unimodal, roughly symmetric (slightly left skewed), in spite of some outliers (Fig.~\ref{fig:3}a). It ‘nearly’ approximates the maximum-likelihood fitting of a normal curve (red line) with mean ($\mu$) = $-1.53$~m and standard deviation ($\sigma$) = $0.27$~m. However the Shapiro-Wilk normality test rejects the null hypothesis of a gaussian distribution ($p-value=0.0005$). On the other hand, the Q-Q plot (Fig.~\ref{fig:3}b) shows deviations from the theoretical normal distribution (red line) between one ($68.27\%$ of the sample) and two ($95.45\%$) standard deviations from the mean. The S-shaped empirical distribution recalls its left skew.

The global vertical distribution of finds resemble that of the faunal assemblage, due to the weight of the latter on the sample data ($n=471$). The distribution of the lithic artifacts is more left skewed, while the small sample of micromammals follows a multimodal distribution with a prominent peak at $-1.2$~m (Fig.~\ref{fig:3}a). Although the difference in size of the two samples, it is worth notice that the mean value of the vertical distribution of \emph{A. ruffoi} is very close to that of the lithic artifacts (Welch Two Sample t-test $p-value=0.5803$).

The vertical distribution of finds in SU D is globally unimodal, slightly left skewed, with one peak at $-2.10$~m and no outliers (Fig.~\ref{fig:3}c). Although the distribution is close to the best fitting normal curve (red line) with $\mu=-2.34$~m and $\sigma=0.33$~m, the Shapiro-Wilk test rejects the hypothesis of normality ($p-value=2.497e-12$). The Q-Q plot (Fig.~\ref{fig:3}d) shows a more dispersed distribution with respect to the former one. Its steeper line follows the theoretical normal distribution within one standard deviation from the mean ($68.27\%$ of the sample).

Compared with the global distribution, the vertical distribution of lithic artifacts slightly skews to the right. Nevertheless, the Shapiro-Wilk test fails to reject the normality hypothesis ($p-value=0.2742$). On the other hand, the micromammal distribution is multimodal and slightly shifted to the right (Fig.~\ref{fig:3}c). Its mean ($-2.284$~m) is quite close to the mean of the lithic sample ($-2.423$~m). However, the Welch t-test rejects the hypothesis of two equal sample means ($p-value=0.0141$). If we cannot state that the two distributions have the same mean, we remark the highest density of both the assemblages at around $-2.5$~m.

\begin{figure}
  \centering
  \includegraphics[width=1\textwidth]{../artwork/Fig3.pdf}
  \caption{Vertical distribution of finds in SUs C (a,b) and D (c,d).}
  \label{fig:3}
\end{figure}

As for the vertical distribution of the identified taphonomic features on the lithic and faunal assemblages, Fig.~\ref{fig:4}a,b illustrates the overall distribution of patinae on the lithic artifacts, across SUs C and D. The histogram shows the increasing number of finds between the two sedimentary units. This trend is reflected as well in the rise of Fe-Mn patinated artifacts ($41\%$ in SU C and $45\%$ in SU D), compared to non-patinated ones (respectively $50\%$ and $48\%$). The kernel density estimation (blue and green lines) shows a slightly higher occurrence of patinated artifacts at the lower part of the sequence (below $-2.5$~m), whereas in SU C (up to $-2$~m) there is no evident preference in the vertical distribution of patinae. The higher density of patinated artifacts, linked to the concentration of lithics observed in Fig.~\ref{fig:3}c at about $-2.5$~m, can be localized in a restricted spot at the bottom right corner of the excavated area (Fig.~\ref{fig:4}a).

Restricting the analysis to SU C, the vertical distribution of coordinated macrovertebrate fossils analyzed by \cite{Bagnus2011} spans $71\%$ of the elevation range of the complete assemblage from the same SU. However, being only the $29\%$ of the population, we acknowledge that our sample cannot be considered representative.

The densities of the low and medium rate of oxides resemble the general distribution (Fig.~\ref{fig:4}c,b). Low values follow a ‘nearly’ normal distribution (Shapiro-Wilk normality test $p-value=0.2186$). The density of high oxidized remains ($54\%$ of the sample) draws a left skewed distribution, with a peak at about $-1.3$~m; whereas fossils with a medium degree of oxides are skewed to the right. However there is no clear preference for oxides to occur deeper in the sequence. The mean values are very close to each other and lower values of oxides are more dense at the bottom of the SU.

As for the distribution of the three taphorecords, the prominent peak of TR1 at $-1.6$~m (Fig.~\ref{fig:4}f) contrasts with a more distributed and mixed distribution of the second and third groups of fossils. However, the very low frequency of TR1 ($n=9$) limits further analyses.

\begin{figure}
  \centering
  \includegraphics[width=1\textwidth]{../artwork/Fig4.pdf}
  \caption{3D and vertical distributions of Fe-Mn patinae on the lithic assemblage from SUs C and D (a,b); oxides (c,d) and taphorecords (e,f) from the faunal sample.}
  \label{fig:4}
\end{figure}

Although our study is constrained by the small sample of fossils and by its limited spatial extension to SU C, the analysis of the vertical distribution of Fe-Mn oxides in the faunal and lithic assemblages does not show any clear global pattern. Indeed, even taking into account the localized cluster of artifacts at the very bottom of SU D (Fig.~\ref{fig:4}a), the process responsible for the distribution of Fe-Mn oxides seems to operate indistinctly through the complete stratigraphic sequence, with no explicit preference for lower elevations.

With no evidence for strong right skewed distributions of finds in SUs C and D, we have reasons to exclude the covariate effect of gravity in the observed spatial pattern. The subsequent point pattern analyses are directed to the study of the 2D spatial distribution of fossils, lithics and their taphonomic status.

\section{Point pattern analysis}

The observed patterns of the archaeological and paleontological remains within SUs C and D, as well as the patterns of taphonomic features recognized on them, have been treated as realizations of spatial point processes, i.e. site formation and modification processes.

Indeed, a spatial point pattern is generally defined as the location of events generated by a point process, operating simultaneously at different scales: a first-order global scale and a second-order local scale \citep{Bailey1995}. The former results from the frequency (density) of events within a bounded region; the latter results from spatial dependency between points, e.g. from a tendency for values of the process at nearby locations to interact with each other. Three different types of interpoint interaction are possible: random (or Poisson); regular and cluster. Regular patterns are assumed to be the result of inhibition processes, while cluster patterns are the result of attraction processes. Therefore, two main issues of interest are explored by spatial point pattern analyses: the distribution (density) of entities in space and the existence of possible interactions between them \citep{Ord1972}.

First-order effect in the observed point-pattern is generally non-parametrically evaluated by means of kernel density estimation \citep{Diggle1985}. As an average density of points in the study region, intensity informs about uniform or inhomogeneous distribution of events.

Multiple scales of second-order patterning and the probability of a stochastic occurrence are explored by the Ripley's $K$ summary function \citep{Ripley1976,Ripley1977} and derivates, for both univariate and bivariate point patterns. The $K$ function is designed to identify the relative aggregation and segregation of point data at different scales. The univariate $K(r)$ function measures the expected number of events found up to a given distance $r$ around an arbitrary event. By comparing the estimated value $\hat{K}(r)$ to its theoretical Complete Spatial Randomness (CSR) value, it is possible to assess what kind of interaction exists between events. The bivariate function, or cross-type $K_{ij}(r)$ function seeks to evaluate, at each distance $r$, the spatial relation between two types $ij$ of observed events. In this case, the definition of the null hypothesis uses a randomization technique of either the location of one of the types (random shift hypothesis), or the type itself of the event at each point, preserving the original location (random labeling hypothesis) \citep{Goreaud2003}. The former aims to evaluate the spatial relationship between patterns of two independent processes, while the latter assumes the same process in determining the pattern for different types.

Especially in small dataset, the estimation of correlations between points is biased by edge effects, arising from the unobservability of points outside the window of analysis. In order to reduce that bias, we implement here Ripley's isotropic edge correction \citep{Ripley1988,Ohser1983}.

Monte Carlo simulations \citep{Robert2004} are used to generate pointwise critical envelopes of random expected values for the null hypotheses, providing an adequate level of statistical significance. We choose a small significance level ($\alpha=0.01$ obtained with $199$ simulations), due to the higher possibility of committing a Type 1 error by testing our hypotheses. Values of the empirical distribution (black solid line) are plotted against the theoretical Poisson distribution (red dotted line) and the simulated global envelope of significance (grey area). For $K(r)$, when the solid line of the observed distribution is above or below the shaded grey area, the pattern is significantly clustered (points are closer together than would be expected for a complete random pattern) or dispersed. For $K_{ij}(r)$, the benchmark value $\pi r^2$ is consistent with independence between the points of types $i$ and $j$, and does not imply a Poisson distribution.

\subsection{Formation processes}

In order to investigate the processes involved in the formation of the Pirro Nord deposit, we provisionally assume the deposition of each sedimentary unit to be the result of mass wasting events filling the fissure and resulting in the distribution of fossils and artifacts independently of each other.

To test the appropriateness of our working assumption, we first analyze the overall distributions of finds, treated as univariate point patterns. Applying a set of exploratory statistics, we aim to determine the nature of the depositional processes, e.g. if they raise in- or homogeneous distributions. Then, we analyze the relative patterns of the faunal and lithic assemblages from SUs C and D. In this case, we treat the two distributions as multitype point patterns.

The intensity of the lithic and faunal assemblages is non-parametrically estimated by first performing a Gaussian smoothing kernel of their distributions, for both SUs. Likelihood cross-validation bandwidth, which assumes an inhomogeneous process, is selected for each pattern. Edge correction is applied using the method of \cite{Diggle1985}. Then, Berman's $Z_2$ test is used to determine whether or not the intensity depends on a spatial covariate $Z$, assuming that the spatially varying (inhomogeneous) intensity is a function of $Z$. Thus, in order to measure the strength of dependence on the covariate, we use the Receiver Operating Characteristic (ROC) curve. Spatially adaptive smoothing, nearest-neighbour density and scan tests have been used in order to assess for the evidence of hot spots in the intensities of the unmarked point patterns. Estimations of the $K(r)$ and the Kaplan-Meier corrected empty-space $F(r)$ functions provide further methods for the interpretation of the distributions.

Multitype summary functions are used in the analysis of the dependence between points of the two assemblages. In this case, our main research question is whether different types of finds have the same spatial distribution. The cross-type $K_{ij}(r)$ function and the Kaplan-Meier corrected nearest-neighbour $G_{ij}(r)$ function are used to estimate the association between points of types $i$ and $j$, for any pair of types of finds. Positive spatial correlation between the two types of finds would suggest that lithic artifacts are more likely to be found close to fossils than would be expected for the hypothesis of \emph{independence}. It would confirm the field observations about their close stratigraphic association and further support our hypothesis that both patterns are the realization of one depositional process. On the other hand, segregation of the two patterns is equivalent to variation in the probability distribution of types. Segregation could be interpreted as the expression of preferential/differential depositional processes. In this case, more detailed analyses would be necessary.

\subsection{Post-depositional processes}

In order to evaluate the degree of post-depositional disturbance of the deposit, the spatial dependence of observed taphonomic features is assumed to be the expression of a related diagenetic process. Measured phenomena that are closer together in space, tend to be more related than those further apart \citep{Tobler1970}.

Like in applications of point pattern analysis in spatial epidemiology \citep{Diggle2003,Gatrell1996}, we distinguish between \emph{cases} and \emph{controls}. The distribution of cases of a certain taphonomic alteration can be regarded as the realization of a diagenetic point process, whereas control points refer to non-altered remains. In a conditional analysis of a spatial case-control study the locations are fixed covariates, and the taphonomic status is treated as a random variable. The simplest null model (\emph{random labelling}) is that the taphonomic status of each find is random, independent and with constant risk of occurrence.

Spatial correlations of diagenetic alterations on the lithic and faunal assemblage are explored by the $K_{ij}(r)$ function, random labelling the pair case/control of Fe-Mn oxidation. We assume in this case that an independent process (karst water circulation), subsequent to the initial event responsible for the accumulation of the finds in each SU, determined their preservation status. Positive deviations from the null hypothesis, suggest that cases are more likely to be found close to controls than would be expected if their status was randomly determined. On the other hand, negative deviations would indicate segregation between cases and controls. Thus, it would suggest that the action of post-depositional water-related processes could have locally reworked the original distribution, determining the altered status of the remains.

All the spatial analyses were performed using the \emph{spatstat} package \citep{Baddeley2015} in \textsf{R} statistical software \citep{RCoreTeam2015}.

A repository containing a compendium of data, source code and text \citep{Marwick2016} is archived at the DOI: \href{http://dx.doi.org/10.5281/zenodo.55736}{10.5281/zenodo.55736}.

\section{Results}

\subsection{Formation processes}

Fig.~\ref{fig:5}c,d shows the smoothing kernel estimation of the faunal assemblage intensity respectively in SUs C and D. Lithic artifacts and micromammals remains of \emph{A.~ruffoi} are superimposed on it. The visual assessment of the plot suggests positive spatial association between the three types of finds. Higher intensities in the distributions are evident at specific values of the \emph{x} coordinate ($6<x<7$ and $8<x<9$), in both the sedimentary units. Yet, a concentration of artifacts, already observed in Fig.~\ref{fig:3}c and \ref{fig:4}a,b, is evident at the lower right corner of SU D (Fig.~\ref{fig:5}e). Such higher densities of finds are clearly shown as well in the scatterplots of the projected third coordinate (Fig.~\ref{fig:5}a,b). Notably, the thickness of the sedimentary unit cannot be accounted to be responsible for those hot spots with higher density of finds. Neither the apparent inhomogeneous intensities along the $x$ axes are supported by the ROC curves (Fig.~\ref{fig:5}e,f). Even if Bermans's $Z_2$ tests suggest significant evidence of dependence on the $x$ covariate, the ROC curves show that it does not have strong disctiminatory power.

\begin{figure}
  \centering
  \includegraphics[width=1\textwidth]{../artwork/Fig5.pdf}
  \caption{Scatterplot of finds from SU's C (a) and D (b); Smooth density estimation of the faunal assemblage and distribution of lithic artifacts and \emph{A.ruffoi} remains in SU's C (c) and D (d); ROC curves for the covariate $x$ coordinate in SU's C (e) and D (f).}
  \label{fig:5}
\end{figure}

Fig.~\ref{fig:6}a,d shows the resulting $p-values$ of likelihood ratio scan test statistic. The test detects differences in the densities of the distributions, showing zones with high abundance of finds. The estimated homogeneous $\hat K(r)$ and $\hat F(r)$ functions are consistent with this result. For both SUs C (Fig.~\ref{fig:6}b,c) and D (Fig.~\ref{fig:6}e,f) they suggest strong deviation from the null hypothesis of CSR towards aggregation, at any scale.

\begin{figure}
  \centering
  \includegraphics[width=1\textwidth]{../artwork/Fig6.pdf}
  \caption{p-values of the likelihood ratio scan test, with logarithmic colour scale, for SUs C (a) and D (d); pointwise envelopes of the homogeneous $K(r)$ and $F(r)$ functions for unmarked finds from SUs C (b,c) and D (e,f).}
  \label{fig:6}
\end{figure}

In analyzing a point pattern, it is confounding and it may be impossible to distinguish between clustering and spatial inhomogeneity \citep{Baddeley2015}. Given the context of the site, and the results of our non-parametric analyses, we proceed considering the distributions of finds as the results of cluster homogeneous processes. The bivariate version of the homogeneous $K_{ij}(r)$ and $G_{ij}(r)$ functions allows us to statistically test the hypothesis of aggregation between the types of remains.

In Fig.~\ref{fig:7}, the top line of panels (a,b,c) shows the ordinary estimations of the $K$ function for the three types of finds (Fauna, Lithic and \emph{A.~ruffoi}) from SU C. Panel~\ref{fig:7}a resembles Fig.~\ref{fig:6}b and indicates statistical significant clustering of the faunal remains for any values of $r$. The lithic assemblage shows as well a significant cluster tendency, for $r>0.1$, while it fails to reject the null hypothesis of CSR for lower values. Instead, the estimated $\hat{K}(r)$ for the micromammals shows aggregation, but, for all values of $r$, we cannot state that the distribution is not random. This result might reflect the random displacement applied to the micromammal point pattern.

The middle and bottom lines of panels in Fig.~\ref{fig:7} show estimations of the homogeneous cross-type $K$ and $G$ functions for all pairs of types $i$ and $j$. Interestingly, Fig.~\ref{fig:7}d suggests positive spatial correlation between lithic and faunal remains at any values of $r>0.05$. The corresponding $G_{ij}(r)$ function measured the cumulative distance from each point of type $i$ (Lithic) to the nearest point of type $j$ (Fauna). It shows (Fig.~\ref{fig:7}g) that the nearest-neighbour distances are significantly shorter than expected, but we cannot reject the hypothesis of independence between fossils and artifacts. However, the short scale of the function suggests that any artifact is surrounded by fossils. This result statistically confirms the stratigraphic association of artifacts and fossils, previously based on field observations. On the other hand, deviations between the $\hat{K}_{ij}(r)$ function and the benchmark $\pi r^2$ suggest segregation between lithics and \emph{A. ruffoi} specimens, but the hypothesis of independence between the two types is more significant (Fig.~\ref{fig:7}e,h). Conversely, the small mammal assemblage is closer to the rest of the fossils than expected for independent distributions, for $r>0.2$. For lower values of $r$, the $K$ and $G$ functions fail to reject the hypothesis of independence.

\begin{figure}
  \centering
  \includegraphics[width=1\textwidth]{../artwork/Fig7.pdf}
  \caption{Pointwise envelopes of the homogeneous cross-type $K_{ij}(r)$ and $G_{ij}(r)$ functions for all pairs of types $i$ and $j$ in SU C.}
  \label{fig:7}
\end{figure}

The top line of panels in Fig.~\ref{fig:8} (a,b,c) shows estimations of the $K(r)$ function for the three types of finds from SU D. Panel~\ref{fig:8}a confirms the same clustering trend of the faunal assemblage. Analogous to the distribution of finds from SU C, the global pattern is mostly weighted on the faunal assemblage (Fig.~\ref{fig:6}e). Conversely, in SU D the distribution of lithics shows stronger significant clustering for $r>0.1$. Again, the resulting $\hat{K}(r)$ for the micromammal assemblage suggests a statistically insignificant aggregation tendency for all values of $r$, but $0.4<r<0.5$. In contrast to the previous result, estimations of the $K_{ij}(r)$ function support significant positive correlation between the lithic artifacts and the \emph{A.~ruffoi} remains (Fig.~\ref{fig:8}e). Thus, they occur closer than expected in the case of independent distributions. Panel~\ref{fig:8}f shows the same positive correlation also between micro- and macro-mammals for $r>0.2$. The panels~\ref{fig:8}d,g show as well a significant positive aggregation between lithics and fossils for values of $r>0.1$. In addiction, the estimated $\hat G_{ij}(r)$ function offers a closer view of the distribution. For values of $r<0.1$, it fails to reject the hypothesis of independence.

\begin{figure}
  \centering
  \includegraphics[width=1\textwidth]{../artwork/Fig8.pdf}
  \caption{Pointwise envelopes of the homogeneous cross-type $K_{ij}(r)$ and $G_{ij}(r)$ functions for all pairs of types $i$ and $j$ in SU D.}
  \label{fig:8}
\end{figure}

\subsection{Post-depositional processes}

To achieve our second objective, namely to evaluate the degree of post-depositional disturbance of the deposit, we first analyzed spatial distribution of oxides on the lithic and faunal assemblages independently, then we moved to a comparative analysis. We are particularly interested in the spatial distribution of Fe-Mn oxides (cases) compared with the absence of them (controls).

Fig.~\ref{fig:4}a does not suggest segregation of patinated and not-patinated lithics. If we perform random labeling of the presence of Fe-Mn (spotted and covering) with its absence in both the stratigraphic units, the outputs of the cross-type function (Fig.~\ref{fig:9}a,~b) show that the observed altered artifacts are, with a $0.01$ level of significance, randomly and independently located in SU C. The positive discrepancy between the estimated $\hat K_{ij}(r)$ and the benchmark $\pi r^2$ indicates aggregation of cases and controls, but it lies within the grey envelope of the \emph{random labeling} hypothesis. Conversely, patinated and non-patinated lithics in SU D appear to be closer to each other than expected for the null hypothesis. In this unit the observed $\hat K_{ij}(r)$ function over-exceeds the envelope at values of $r>0.4 m$, hence it indicates statistically significant aggregation. Such pattern statistically confirms the visual assessment of Fig.~\ref{fig:4}a. Consequently, oxidized and non-oxidized artifacts most probably occur in SU D well aggregated in space, while their aggregation is not statistically significant in the above unit.

We could not compare the oxidation patterns between lithics and fossils from SU D, because the taphonomic analysis of \citet{Bagnus2011} did not include fossils from this unit. Thus, we focused our analysis on SU C.

The distribution map (Fig.~\ref{fig:4}c) does not suggest any evident pattern. When we apply random labelling of the absence of oxidation with the medium and high degrees of its presence, the output of the bivariate $K_{ij}(r)$ function shows a segregation tendency between them, but it is not statistically significant. (Fig.~\ref{fig:9}c). A random and independent distribution of oxides is more plausible.

Finally, Fig.~\ref{fig:9}d shows the result of the $K_{ij}(r)$ function, random labeling the cases (medium and high degrees) and controls (absent or low degree) of Fe-Mn oxides on the lithic and faunal assemblages from SU C. The empirical values of the cross-type function are balanced on the theorerical expectation for complete spatial independence (red line). It clearly lies inside the grey envelope of significance. Therefore, our analysis shows an independent spatial distribution of Fe-Mn patinated and non-patinated lithic artifacts and fossils from SU C. In the lower unit (SU D), where Fig.~\ref{fig:4}b indicates higher density of oxidized artifacts, estimations of the cross-type $K$ function suggests that they occur closer than expected to fresh ones.

\begin{figure}
  \centering
  \includegraphics[width=1\textwidth]{../artwork/Fig9.pdf}
  \caption{Pointwise envelopes of the homogeneous bivariate $K_{ij}(r)$ function, random labeling \emph{cases}/\emph{controls} of patinated lithics in SUs C (a) and D (b); oxidated fossils in SU C (c); Fe-Mn oxides on the lithic and faunal assemblages from SU C (d).}
  \label{fig:9}
\end{figure}

\section{Discussion}

The Early Pleistocene site of Pirro~Nord (fissure P13) has yielded evidence for one of the earliest occurrences of hominins in Europe. The importance of the evidence calls for a multivariate taphonomic analysis in order to establish the nature of the processes involved in the formation of the deposit and the degree of its post-depositional disturbance. We address that need by investigating the spatial association of the archaeological and paleontological remains, as well as the spatial distribution of artifacts and fossils with diagenetic alterations. We focused our analysis on the lower stratigraphic units C and D, since they provide the most significant corpus of finds and they have been studied with the same research protocol.

\subsection{Formation processes}

Non-parametric analyses have been carried out in order to characterize the processes responsible for the formation of the deposit. Then, we accounted for the relative spatial pattern of the different types of finds.

The vertical distribution of the archaeological and paleontological assemblages does not appear to be affected by strong gravitational effects. On the other hand, it resemble a 'nearly' normal distribution and suggests a very close mean occurrence of lithic artifacts and \emph{A. ruffoi} remains, despite the small sample of micromammals (Fig.~\ref{fig:3}). A visual interpretation of the projected third coordinate (Fig.~\ref{fig:5}a,b) also suggests that the intensity of finds is not a function of the covariate $z$. Moreover, higher densities are not linked to the thickness of the stratigraphic units. They are clearly localized at values of $6<x<7$ and $8<x<9$ in both the SUs, as shown also by Fig.~\ref{fig:5}c,d and \ref{fig:6}a,d. Indeed, the Berman's $Z_2$ test for the dependence of the point process on the spatial covariate $x$ failed to reject the null hypothesis for SU C ($p-value=0.0044$) and D ($p-value=1.913e-14$). However, even if it suggests significant evidence that the intensity depends on some covariate, the effect of that covariate could still be weak. ROC curves (Fig.~\ref{fig:5}e,f) indicate that the $x$ coordinate does not have disciminatory power.

\cite{Bartlett1963} showed that it is possible to formulate a point pattern which can be equally interpreted as a Poisson inhomogeneous process, or a homogeneous cluster process. According to our non-parametric analyses, we proceeded under the assumption that the processes involved in the formation of the Pirro Nord (P13) deposit are homogeneous and clustered. The scan tests in Fig.~\ref{fig:6}a,d show hot spots of points, mostly localized between $6<x<7$ and $8<x<9$. The cluster correlation between all the finds is significantly confirmed by the estimations of the $K$ and $F$ functions (Fig.~\ref{fig:6}). The first lines of panels in Fig.~\ref{fig:7} and \ref{fig:8} offer a type-based view of these patterns. Indeed, the estimated $\hat K(r)$ functions of the faunal assemblage (Fig.~\ref{fig:7}a and \ref{fig:8}a), which constitute the bigger part of the analyzed sample of data, resemble the results for the complete populations (Fig.~\ref{fig:6}b,e). The lithic assemblage also shows significant aggregation; while the small sample of \emph{A. ruffoi} falls inside the envelope of CSR.

Faunal remains and lithic artifacts show some overlapping when evaluated by means of Gaussian smoothing kernel (Fig.~\ref{fig:5}c,d). Positive spatial association between fossils and lithics is statistically confirmed by the cross-type $K$  and $G$ functions for the examinated SUs (Fig.~\ref{fig:7} and \ref{fig:8}). Fossils and artifacts tend then to occur aggregated with each other (they are closer than expected for a independent process). Significant spatial proximity is also shown between artifacts and micromammal remains, especially in SU D.

According to the results of our analyses, the stratigraphic and spatial association between the types of remains should be considered as the result of the same formation process. Finds occur in the clayey-sand sediment together with a high number of angular to sub-rounded gravels and boulder-sized rock clasts. Such a stratigraphic setting suggests repeated mass-wasting processes (at least two events, represented by SUs C and D, which were included in this study) with a low degree of reworking in a relatively short span of time \citep{Arzarello2012}. Rapid-moving and chaotic water-laden masses, such as mud-flows or earth-flows, of soilwash and rock rubble with fossils and artifacts \citep[p.~46]{Butzer1982}, could have been triggered by intense rainfalls and became trapped in the karst sink-hole directly opening to the outside. The sedimentary fill would have derived from the top, by gravity, directed into the empty space between the large limestone blocks that made up the internal structure of the fissure. The thickness of the layers is likely correlated to the intensity of such events.

On the other hand, the clustered distribution of all the finds cannot be linked to the thickness of the stratigraphic units. The big blocks of calcarenite, which in some places transect the stratigraphic units (Fig.~\ref{fig:1}), created a complex internal structure and might have influenced the direction of sediment accumulation. However, sedimentation rate, driven by the rugged topography of the site, does not seem to be spatially associated with the localized hot spots (Fig.~\ref{fig:6}a,d). Thus, clustering might have been a correlated effect of the formation process.

The presence of partially articulated vertebrate skeletal elements and their general low degree of weathering indicate fast burial and transport of bones from nearby locations \citep{Bagnus2011}. A close spatial proximity between the original location of the finds and the karst fissure, as well as a relatively fast burial, is also corroborated by the unrounded ridges and edges of the lithic artifacts and by the technological consistency of the assemblage \citep{Arzarello2015}.

In conclusion, our spatial statistics analyses confirm the field observations about the spatial association of archaeological and paleontological remains.

\subsection{Post-depositional processes}

After dealing with our first research question (to examine the spatial distribution of finds in the context of the site formation processes), our analyses were particularly directed to test the hypothesis of post-depositional processes reworking the deposit. We assume that the spatial aggregation of taphonomic surface alterations, and their relative segregation compared to non-altered finds, would indicate the localized activity of diagenetic agents.

We focused more on the spatial distribution of oxides, because traditional explanations for the development of Fe-Mn patinas on the surface of flint refer to the deposition of various iron and manganese oxides and hydroxides out of soil water \citep{Stapert1976}. Similarly, the origin of manganese coatings on fossils, in karst environments, may derive from circulating water, or from the manganese present in the surrounding limestone rock dissolved by groundwater \citep{Hill1982}.

The vertical distribution of oxides on the lithic artifacts and the sample of faunal remains (Fig.~\ref{fig:4}) spans the complete stratigraphic sequence and apparently shows a gradual increase through the lower layers, especially in the lithic assemblage. Intensity of oxides is indeed more likely proportional to the density of finds and not related to the depth.

By applying a set of spatial statistics (namely cross-type $K_{ij}(r)$ function) to the archaeological and paleontological remains, we searched for evidence of localized areas, which might have been subjected to the presence of water, especially water-flows.

In SU C, the spatial distribution of Fe-Mn patinas on lithics and fossils is, with a certain degree of significance, the result of independent processes (Fig.~\ref{fig:9}a,c,d). We cannot state that there is aggregation (spatial proximity) between oxidized finds and fresh ones. Neither the results of the bivariate $K$ function, random labeling the cases/controls of Fe-Mn coating, show segregation, which is indicative of spatially defined diagenetic processes. In contrast, in SU D (Fig.~\ref{fig:9}b), patinated and fresh artifacts occur significantly spatially aggregated to each other for values of $r>0.4 m$. They occur closer than expected by an independent process at bigger scale. However, the pattern is, with a certain confidence level, independent.

\citet{Rottlander1975} identified a possible different cause of Fe-Mn coatings in the iron that is already present in the flint. In this light, the spatial association of flint artifacts with and without patination also depends on the chemical and microstructural composition of the raw material itself. On the other hand, the same oxides affecting a good percentage of finds, have been equally found broadly scattered on the numerous clasts of calcarenite that are included in the matrix, thus supporting an external origin of the Fe-Mn coating process.

The content of water and organic matter in the sedimentary body could be responsible for the randomly diffuse Fe-Mn patinations. In the presence of organic matter, indeed, it is likely that the release of organic acids will accelerate patination on chert \citep{Burroni2002}. Moisture of the sedimentary body could also be accounted for the wide random spread of Fe-Mn coatings.

We did not find statistically significant evidence of aggregation of oxidized records compared to non-oxidized ones (Fig.~\ref{fig:9}c); thus, we can exclude the assumption of localized concentration of water, which is included in the hypothesis advanced by \citet{Bagnus2011} for the presence of interstitial flows reworking the deposit.

Due to the small sample size, we did not apply spatial analysis to the distribution of the three taphorecords. However, Fig.~\ref{fig:4}e,f suggests that fossils from the TR1 group occur spatially aggregated with fossils from the TR2 and TR3 groups. Considering the distribution of Fe-Mn oxides on the lithic and faunal assemblages, the re-elaborated TR1 \citep[\emph{sensu}][]{Fernandez-Lopez1991,Fernandez-Lopez2007,Fernandez-Lopez2011} might not be associated with the reworking action of water-flows and might be more likely correlated to random and limited rearrangement of parts of the sedimentary matrix.

A possible cause of some localized movement of sediments could be the rock falls from the vault of the karst fissure, during the deposition of SUs C and D. As showed in Fig.~\ref{fig:1}, an abrupt increase in the number of boulder-sized rocks is observed within the lower layers. Moreover, rock falls caused most of the post-depositional fractures on the faunal assemblage \citep{Bagnus2011}. Such intense erosional process could most likely be correlated to the seismic activity of the region \citep{Bertok2013}.

Results of our analyses suggest that post-depositional taphonomic alterations occurred with a certain significance as a result of independent processes. However, keeping a cautious approach to spatial analysis, a documented point pattern can be most realistically thought of as the result of multiple processes heterogeneously working at different scales \citep{Bevan2013}. Multiple or repeated post-depositional processes could obliterate contemporaneous or preceding patterns, resulting in a final random distribution of the record.

Moreover, karst site formation processes are highly dependent on the structure and extension of the overall karstic system, as well as on the surrounding environment. The lack of information about the original characteristics of the system and the reduced area of excavation strongly limit the analysis.

Furthermore, although the need for considering three dimensional distributions in site formation processes study, spatial point pattern statistics are at the moment not fully equipped to analyze three-dimensional patterns, especially when the study-area corresponds to a three-dimensional volume with a complex shape such as a karstic structure.

On the other hand, "one must look to non-spatial evidence to corroborate or disprove theories about spatial processes" \citep[p.~8]{Hodder1976}. The integration with other taphonomic disciplines reinforces the results of spatial analyses and outlines new opportunities for point pattern analyses. As recently remarked \citep{Cobo-Sanchez2014}, taphonomic research should be multivariate \citep{Dominguez-Rodrigo2010} and it should include spatial analysis as a heuristic tool in the interpretation of site integrity. This is especially demanding when the research questions deal with past human behaviour and even more so when site dating is based on the stratigraphic association of artifacts and fossils.

\section{Conclusions}

The Early Pleistocene site of Pirro~Nord~13 provides evidence of the earliest human presence in Western Europe. Lithic artifacts have been found in a karst fissure filling, together with late Villafranchian/early Biharian paleontological remains.

The main goals of our study were:
\begin{enumerate}
  \item to investigate the depositional processes involved in the formation of the deposit and
  \item to assess the degree of any potential post-depositional reworking of the archaeological and paleontological remains.
\end{enumerate}

The integration of spatial point pattern analyses with previous taphonomic studies on the faunal and lithic assemblages allowed us to test different hypotheses of site formation and modification processes.

On the basis of our analyses,
\begin{enumerate}
  \item we consider the deposit to be the result of subsequent events of some type of mass-wasting process, such as a mud-flow or earth-flow, carrying rock rubble with fossils and artifacts. The applied set of spatial analyses confirms, with an adequate level of statistical significance, the assumption, based on field observations, regarding the spatial association between the finds.
  \item Based on our taphonomic point pattern analyses of several diagenetic features on the lithic and faunal assemblages, we reject the hypothesis of a substantial post-depositional reworking and mixture of the sedimentary deposit and we corroborate the stratigraphic integrity of the Pirro Nord 13 site.
\end{enumerate}

Finally, the present study answers the need for a taphonomic perspective in spatial analysis, by applying well developed quantitative methods in spatial statistics. Point pattern analysis can be very flexible and useful in the investigation of both cultural and taphonomic processes. Until now it has found limited application on taphonomic studies, but, as our study demonstrates, it offers new analytical opportunities to the multidisciplinary study of the complex processes that operate in the formation and modification of archaeological sites. It allows analysts to test multiscalar patterns and to model the taphonomic processes underlying archaeological distributions, which are otherwise difficult to identify from the simple visualization of maps, especially for those sites characterized by complex geo-stratigraphic settings.

\section*{Acknowledgements}

DG was supported by the European Research Council Starting Grant No. 283503 (PaGE) awarded to K.~Harvati.

Fieldworks were financially supported by the University of Ferrara and the Municipality of Apricena. Excavations were possible by permission of the Soprintendenza Archeologica della Puglia and Ministero dei Beni e delle Attività Culturali; thanks to the collaboration of Dott.ssa Anna Maria Tunzi; and by courtesy of Gaetano and Franco dell'Erba.

We would like to thank all those who attended the excavation campaigns for their indispensable cooperation; Julie Arnaud for her * and contribution to the 3D modeling of the site; and Nick Thompson for English proofreading.

We are grateful to Vangelis Tourloukis and two anonymous reviewers for their critical discussions and many constructive comments that helped to improve this manuscript.

\section*{References}

\bibliographystyle{elsarticle-harv}
\bibliography{pirronord_rr}

\end{document}

ext: tex


text:
\epigraph{``Yes, it is easy to see that nearly six years of magical education has not been wasted on you ... Ghosts are transparent"}{Severus Snape}

Numerical Weather Prediction focuses on taking current observations of weather and processing this data with computer models to forecast the future state of weather. Knowing the current state of weather is just as important as the numerical computer models processing the data. Current weather observations serve as input to the numerical computer models through a process known as data assimilation to produce outputs of temperature, precipitation, and hundreds of other meteorological elements from the oceans to the top of the atmosphere\cite{nwp_introduction}.

However, weather prediction is not all that it could be. For example, you should remember Hurricane Sandy, which hit New York in 2012. The American Global Forecast System predicted it would not even reach the mainland! This, of course, was wrong, and as a result, 285 people died and 68.7 billion dollars worth of damage was done. But, why is weather forecasting so hard? Forecasting today is done by physical modelling. This means we know the equations that govern weather systems, but we can’t solve them exactly. Instead, we approximate them.  Weather agencies take reams of raw data, feed it into supercomputers, and just number crunch and number crunch. As you can imagine, this requires a significant amount of computational resources. This is not only a cost issue, it also inhibits our ability to produce forecasts of a high temporal resolution; which can be significantly important in predicting extreme weather events, such as tornadoes.

\section{Development of the Idea}
I think it is extremely difficult to pinpoint exactly where the idea originated from; on reflection, however, the spark which really ignited the flame for this project came from a deep interest in the study of the atmosphere and by extension atmospheric science. I initially became interested in this topic based on two distinct factors.

I became interested in the area of atmospheric science from watching a YouTuber named, Dr. Simon Clark. Dr. Simon Clark recently completed a PhD in theoretical atmospheric physics, researching dynamical stratosphere-troposphere coupling, and made a vlog series documenting his experiences. I was extremely intrigued by this series, with keen interest in his videos explaining his research. These videos captivated me and ultimately led me to reading his thesis (titled Quasi-geostrophic influence of the polar stratosphere on the troposphere)\cite{simonclark}.

Secondly, my interest also came from my concern over human caused climate change. I was extremely fascinated by how a gas, such as carbon dioxide, could capture so much more heat than oxygen, even though there is a higher concentration of oxygen in the atmosphere. This ultimately led to the original idea for this project: creating a computational simulation of both Earth and Venus, and using it to determine if Earth could become just as hellish and uninhabitable as Venus if the rate at which we were pumping carbon dioxide into the atmosphere remained constant. There was a few problems with this idea, primarily being that it would require a ridiculous amount of computational resources, and there was no way of ensuring it was accurate. This led me to narrow the scope of the project, and ultimately, I eventually settled on the idea of improving numerical weather prediction in some fashion. From my initial research, as I previously mentioned, atmospheric simulations required a significant amount of computational power. There have been fantastic developments in machine learning over the past decade, and I wanted to apply these advancements to the problem of weather forecasting. Through various discussions with experts in a wide variety of fields, I came to the conclusion that such an endeavour may be fruitful.

\section{Aims}
The principal aim of this project is to create a neural network architecture, which will reduce the amount of computational resources required to generate weather forecasts compared to traditional physics-based models; and one which will provide a similar level of performance when compared to physical models. Access to the source code of any given physics-based forecast system is also notoriously difficult, with enormous price tags being the normal to gain access to time critical weather information. My vision is to make better weather tools available to everyone. Once the software is created, it will be published on the open source platform, known as, GitHub. This will allow programmers, and atmospheric physicists alike to inspect the source code, and to contribute to the development of the software. The software will also be published on Anaconda Cloud. I did this in order to allow, with just a simple command, the installation of the software. I will also develop a series of documentation that will accompany the software, which will accumulate in the development of a website.

\begin{figure}[H]
    \centering
    \includegraphics[width=.75\linewidth]{Images/website}
    \caption{A screenshot of the website for the software.}
    \label{website}
\end{figure}

\section{Organisation of this Work}
The work contained in this project book is broken down into a number of different chapters:

\begin{itemize}
    \item Chapter \ref{architecture_chapter} discusses the various neural network architecture under consideration and the current implementation within the software.
    \item Chapter \ref{implementation_chapter} discusses some of the decisions made on the software development side of the project, such as, choosing an appropriate programming language for the task.
    \item Chapter \ref{benchmarking_chapter} outlines the software benchmarking experiments to be carried out on the software.
    \item Chapter \ref{results_chapter} contains the results of the benchmarking described in the previous chapter.
    \item Chapter \ref{conclusion_chapter} reviews the results and outcomes of the project, discusses the limitations of the software, and lays out a road-map for the continued development and enhancement of the software.
\end{itemize}

max_stars_repo_name: amsimp/papers

max_issues_repo_name: amsimp/papers

max_forks_repo_name: amsimp/papers


text:
\section{Classic Models}

% linear regression
\subsection{Linear Regression}

The prediction is:
\begin{equation}
    \hat{y} = \theta_0 + \sum_{i=1}^n \theta_n x_n = \columnvector{\theta}^\top \columnvector{X}
\end{equation}

The cost function is:
\begin{equation}
    J(\columnvector{\theta} ) = \text{MSE} = \frac{1}{m} \sum_{i=1}^m \left(\columnvector{\theta}^\top  \subscription{X}{i} - \subscription{Y}{i} \right)^2
\end{equation}

The close solution is:
\begin{equation}
    \hat{\theta} = (\columnvector{{X}^\top \columnvector{X}} )^{-1} \columnvector{X}^\top \columnvector{Y}
\end{equation}

The $\hat{\theta}$ is usually calculates using pseudoinverse of $\columnvector{X}$ as $\columnvector{{X}}^{+} \columnvector{Y}$ using SVD which has time complexity $O(n^2)$.

The gradient vector of $J$ is calculated as:
\begin{equation}
    \begin{aligned}
    \dpd{J(\theta)}{\theta_j} &= \frac{2}{m} \sum_{i=1}^m \left( \columnvector{\theta}^\top \subscription{X}{i} - \subscription{Y}{i} \right) \subscription{X}{i}_j \\
    \nabla_\theta J(\theta) &= \frac{2}{m} \columnvector{X}^\top (\columnvector{X} \columnvector{\theta} - \columnvector{Y} )
    \end{aligned}
\end{equation}

Linear regression model is convex, so it will be guaranteed to find global minimum using gradient descent.

% logistic regression
\subsection{Logistic Regression}

\cindex{logistic regression} is a binary classifier. It uses the \cindex{sigmoid} function $\displaystyle \sigma(t) = \frac{1}{1 + \exp{(-t)}}$ and assigns probability to $\columnvector{X}$ that $\hat{p} = \sigma(\columnvector{X}^\top \columnvector{\theta} )$. The prediction is:
\begin{equation}
	\hat{y} = \begin{cases}
		0 \text{, if } \hat{p} < 0.5 \\
		1 \text{, if } \hat{p} \geq 0.5
	\end{cases}
\end{equation}

So the prediction is $1$ if $\columnvector{X}^\top \columnvector{\theta} > 0$ and $0$ if it is negative.

The \cindex{log loss} cost function is:
\begin{equation}
	J(\theta) = - \frac{1}{m} \sum_{i=1}^m \left\{ \subscription{Y}{i} \log \subscription{\hat{p}}{i} + \left(1- \subscription{Y}{i} \right) \log \left(1 - \subscription{\hat{p}}{i} \right) \right\}
\end{equation}

The derivative is:
\begin{equation}
	\dpd{J(\theta)}{\theta_j} = \frac{1}{m} \sum_{i=1}^m \left( \sigma \left(\columnvector{\theta}^\top \subscription{X}{i} \right) - \subscription{Y}{i} \right) \subscription{X}{i}_j
\end{equation}

\subsection{Softmax Regression}

\cindex{Softmax Regression} is a multiclass version of \cindex{logistic regression}. Each class has its own parameter vector $\columnvector{\theta}_i$. The softmax score for class $k$ is $s_k (\columnvector{X} ) = \columnvector{X}^\top \columnvector{\theta}_k$ and its probability is:
\begin{equation}
    \hat{p}_k = \frac{\exp{{s_k (\columnvector{X})}}}{\displaystyle \sum_{i=1}^n \exp{s_i (\columnvector{X}) } }
\end{equation}

The prediction is:
\begin{equation}
	\hat{y} = \underset{k}{\text{argmax }} \columnvector{X}^\top \columnvector{\theta}_k
\end{equation}

The cost function is measures using \cindex{cross entropy} function:
\begin{equation}
	J(\theta)= -\frac{1}{m} \sum_{i=1}^m \sum_{k=1}^K \subscription{Y}{i}_k \log \subscription{\hat{p}}{i}_k
\end{equation}

The gradient vector is:
\begin{equation}
	\nabla_{\theta_k} J(\theta) = \frac{1}{m} \sum_{i=1}^m(\hat{p}_k^{(i)} - y_k^{(i)}) \mathbf{X}^i
\end{equation}

The cross entropy function is also called \cindex{Kullback-Leibler divergence}.

\subsection{SVM}

\cindex{Support Vector Machine} tries to separate the plane using two parallel lines. It maximize the distance between two lines. The points that these two lines touch is called \cindex{support vector}. The model is sensitive to feature scaling, so normalize the data before training. If all instances are on the opposite sides, it is called \cindex{hard margin classification}.

For each label $i$, $\subscription{Y}{i} =1$ or $\subscription{Y}{i} =-1$.The hard margin classification is:
\begin{equation}
    \begin{aligned}
        \text{minimize: } & \frac{1}{2} \columnvector{W}^\top \columnvector{W} \\
        \text{subject to: } & \subscription{Y}{i} \left( \columnvector{W}^\top \subscription{X}{i} + b \right) \geq 1
    \end{aligned}
\end{equation}

Soft margin classification uses \cindex{slack variable} $\subscription{\zeta}{i} \geq 0$ for each instance which measures how much margin violation is allowed. The optimization becomes:
\begin{equation}
    \begin{aligned}
        \text{minimize: } & \frac{1}{2} \columnvector{W}^\top \columnvector{W} \\
        \text{subject to: } & \subscription{Y}{i} \left( \columnvector{W}^\top \subscription{X}{i} + b \right) \geq 1 - \subscription{\zeta}{i} \\
        & \subscription{\zeta}{i} \geq 0
    \end{aligned}
\end{equation}

\subsection{Decision Trees}

\subsection{Ensemble Learning}

A group of predictors is called \cindex{ensemble}.

\subsubsection{Voting Classifiers}

The class is chosen by majority vote, either by the majority number of votes called \cindex{hard voting}, or by the majority probability called \cindex{soft voting}.

It requires that all models are perfectly independent and have no correlated errors. It is not true because all models are trained on the same data.

\subsubsection{Bagging and Pasting}

It generate multiple classifier by using only one training algorithm but multiple training set. If each sample data is generated using replacement, it is called \cindex{bagging}. So a instance may appear multiple times in bagging. If each sample data is generated without replacement, it is called \cindex{pasting}. All predictors could be trained in parallel and good for scale up.

If it is classification, the result is chosen by majority vote. If it is regression, the result is averaged.

In bagging, it will choose $1 - e^{-1} \approx 63.2\%$ samples. The remaining unchosen one is called \cindex{out-of-bag} instances which could be used as test test.

Bagging is more diverse so it has higher bias but lower variance. Overall bagging is better than pasting.

\subsubsection{Random Patches and Random Subspaces}

\cindex{Random Patch} will sample both training instances and model features. \cindex{Random Subspace} will run all training data but sample only model feature.

\subsubsection{Random Forest and Extra-Trees}

\cindex{Random Forest} is an ensemble of Decision Trees using bagging method (sometimes pasting). \cindex{Extra Tree} use random threshold for each Random Forest feature.

\subsection{Boosting}

\cindex{Boosting} trains multiple predictors sequentially, each trying to correct its predecessor. So it cannot be parallelized.

\subsubsection{AdaBoost}

Every instance is assigned a weight. The first algorithm trains on the data. The weight of misclassified instance will be increased. The second algorithm will run on the updated data set. Once all predictors are trained, the prediction is done by bagging or pasting except that the each predictor has different weight which is the average weight of its prediction.

\subsubsection{Gradient Boosting}

In \cindex{gradient boosting} there is a series of models $\mathbf{F}_i$:

\begin{equation}
	\begin{aligned}
		\mathbf{F}_1 &\approx (\mathbf{X},\mathbf{Y}) \\
		\mathbf{F}_2 &\approx \Big(\mathbf{X},\mathbf{Y} - \mathbf{F}_1 (\mathbf{X}) \Big ) \\
		&\dots \\
		\mathbf{F}_n &\approx \Big(\mathbf{X},\mathbf{Y} - \mathbf{F}_{n-1} (\mathbf{X}) \Big )
	\end{aligned}
\end{equation}

\begin{theorem}
    Gradient boosting is a gradient descent process.
\end{theorem}
\begin{proof}
    Suppose the loss function is:
\begin{equation}
	\mathcal{L}(\mathbf{Y},\mathbf{F}_k) = \frac{1}{2} \sum_i \Big ( y_i - \mathbf{F}_k (x_i) \Big )^2
\end{equation}

Suppose the aggregate function $\mathcal{F}_i = \sum\limits_i \mathbf{F}_i $, the derivative is:
\begin{equation}
	\begin{aligned}
	   \dpd{\mathcal{L}(\mathbf{Y},\mathbf{F}_k)}{\mathbf{F}(x_i)}	  &= \mathbf{F}_k (x_i) - y_i \\
		y_i - \mathbf{F}_k (x_i) &= - \frac{\partial \mathcal{L}(\mathbf{Y},\mathbf{F}_k )}{\partial \mathbf{F}_k (x_i)}\\
		\mathcal{F}_2 &= \mathbf{F}_1 + \mathbf{F}_2 \\
		&= \mathbf{F}_1 + y_i - \mathbf{F}_1 \\
		&= \mathbf{F}_1 - 1 \times \frac{\partial \mathcal{L}(\mathbf{Y},\mathbf{F}_1 )}{\partial \mathbf{F}_1 (x_i)} \\
		&= \mathcal{F}_1 - 1 \times \frac{\partial \mathcal{L}(\mathbf{Y},\mathcal{F}_1 )}{\partial \mathcal{F}_1 (x_i)}
	\end{aligned}
\end{equation}
\end{proof}

\subsubsection{Stacking}

Instead of using soft or hard voting to aggregate the predictions, a model called \cindex{blender} is trained to perform the aggregation. The training set is split into two sets. The first set is used to train predictors. Next the predictors run on second set. The result of second set is used as the input to the blender.

max_stars_repo_name: elvisren/machine-learning-notes

max_issues_repo_name: elvisren/machine-learning-notes

max_forks_repo_name: elvisren/machine-learning-notes