Laurent Najman

Junior Barrera

Petros Maragos

B S Daya Sagar

Dan Schonfeld
email: dans@uic.edu

à la diffusion de documents scientifiques de niveau

recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction to the Issue on Filtering and Segmentation with Mathematical Morphology

Introduction to the Issue on Filtering and Segmentation with Mathematical Morphology

This issue presents novel contributions and introduces the state-of-art in filtering and segmentation methods using mathematical morphology. Historically, mathematical morphology developed a theoretical framework for non-linear image analysis, which from its inception led to important theoretical and practical results. Due to its algebraic foundation and geometrical intuition, the theory of mathematical morphology is very flexible, capable of handling many diverse image types, including (but not limited to) classical multi-dimensional signals as well as general graphs and surfaces.

We received 24 submissions from all over the world, from which we selected 13 articles for publication in this issue spanning a large segment of the research in morphological filtering and segmentation. The first article is an invited paper authored by Jean Serra, who is one of the two founders of Mathematical Morphology and a recent special conference in his honor has been the inspiration for this issue. In a "Tutorial on Connective Morphology", Jean Serra explores morphological operators from the perspective of connections: set connections, connective segmentations, connected operators, hierarchies of partitions and optimal cuts are the five key concepts under scrutiny.

The article "Random Projection Depth for Multivariate Mathematical Morphology," by S. Velasco-Forero and J. Angulo, considers a projection depth map for vector ordering, whereby the authors study a statistical depth function that approximates the Mahalanobis distance from a vector cloud, and use the function in applying mathematical morphological tools to multivalued (vector) images. The article "Active Contours on Graphs: Multiscale Morphology and Graphcuts," by K. Drakopoulos and P. Maragos, proposes two methods for computing geodesic active contours on arbitrary graphs, based on several ideas stemming from differential morphology: (i) approximations to the calculation of the gradient and the divergence of vector functions defined on graphs and uses these approximations to apply the technique of geodesic active contours for object detection on graphs; (ii) appropriate weights are calculated for each edge for which the Riemannian length of a contour can be approximated by the weighted sum of intersections of the contour with the edges of the graph and applied in a graphcut-based solution to the geodesic active contour problem on graphs.

The article "First Departure Algorithms and Image Decompositions into Peaks and Wells", by F. Meyer, presents viewpoint-dependent algorithms to reversibly decompose an image function into peaks and wells based on a generalized shortest path algorithm on weighted graphs. These algorithms include reconstruction closings or floodings and are illustrated with applications.

The article "Salience Adaptive Structuring Elements," by V. Curic et al., introduces salienceadaptive structuring elements that adjust both their shape and their size according to the salience of the edges of an image. Consequently, morphological operators using these kinds of structuring elements degrade edges with high salience to a lesser extent than those with low salience.

J The short article "Fast Morphological Image Processing Open-Source Extensions for GPU processing with CUDA," by M. Thurley and V. Danell, provides a framework that allows 32bit processing of large-scale images to deliver open-source extensions for GPU enabled morphological image processing based on CUDA. The implementation of vHGW algorithms on a GPU to perform fast morphological operations on images by using analog processing is useful for analysis and processing of images with more than 8 bits.

The article "Classification of Remote Sensing Optical and LiDAR Data Using Extended Attribute Profiles," by M. Pedergnana et al., proposes a technique for feature extraction and supervised classification architecture used for classifying multisensor imagery composed of multi/hyperspectral channels and LIDAR data. This approach is novel to remote sensing and emphasizes the EAP-SVM-RF formulation: (i) to prevent Hughes' effects, the Extended Attribute Profile (EAP) approach was applied first to extract features from the multi/hyperspectral and LIDAR data, (ii) vector stacking is applied to the two sets of EAP features, and (iii) SVM and Random Forest (RF) are applied to classify in the stacked-vector space.

The article "Morphology-Based Crack Detection for Steel Slabs," by A. Landström and M. Thurley, is an application of mathematical morphology to continuous casting. Continuous casting is a highly efficient process used to produce most of the world steel production tonnage, but can crack in the semi-finished steel product. The authors present an approach for an automated on-line crack detection system, based on 3D profile data of steel slab surfaces, utilizing morphological image processing and statistical classification.

The article "Multiple Luminaire Identification in Airborne Images of Airport's Approach Lighting Using Mathematical Morphology With Variable Length Structuring Element," by S. Chowdhury et al., proposes analysis of the pattern of the approach-lighting system at an airport, used to provide identification and guidance to pilots from a distance. The proposed approach is based on morphological operators with a variable-length structuring element (VLSE). A novel technique for automatic determination of the VLSE is proposed and allows successful identification of the luminaires from the image data as verified through the use of simulated and real data.

As is apparent from the description above, the articles in this issue span a wide range of important topics using innovative methods based on morphological filters. We hope that the publication of this issue will promote further research and inspire new approaches to address interesting and important challenges using methods and techniques from mathematical morphology.

In closing, we would like to thank all the authors who submitted research manuscripts to this issue, as well as the reviewers for lending their time to evaluate manuscripts and providing valuable comments and suggestions. We would like to extend a special note of gratitude to Dr. Vikram Krishnamurthy, Editor-in-Chief of the Journal of Selected Topics in Signal Processing, and Ms. Rebecca Wollman, IEEE Publications Coordinator, for their constant support.

  .-C. Pinoli and J. Debayle authored the article "Spatially and Intensity Adaptive Morphology" in the context of General Adaptive Neighborhood Image Processing (GANIP). Structuring elements are spatially varying, and their activity is controlled according to the image intensities (level sets are processed at different scales). The advantages of these adaptive operators are highlighted on real-world application examples for image background removal, image restoration and image enhancement.

The article "Efficient Robust d-dimensional Path Operators," by F. Cokelaer et al., deals with morphological operators that use constrained paths as structuring elements. These paths-based operators are very efficient at filtering thin and/or oriented structures in images. Specifically, this article introduces a novel notion of path robustness, allowing the corresponding operators to better and more efficiently deal with noise, in 2D and 3D images. The authors illustrate applications of the proposed approach to blood vessel filtering and segmentation.

The article "One-Dimensional Openings, Granulometries and Component Trees in O(1) per Pixel," by V. Morard et al., provides a novel algorithm to compute one-dimensional openings, granulometries, and 1-D component trees. Its complexity is linear with respect to image dimension and constant with respect to the opening cycle.


Abdullah Almaksour
email: abdullah.almaksour@irisa.fr

Eric Anquetil
email: eric.anquetil@irisa.fr

Improving Premise Structure in Evolving Takagi-Sugeno Neuro-Fuzzy Classifiers

Keywords: Evolving Takagi-Sugeno Neuro-Fuzzy Classifiers Abdullah Almaksour, Eric Incremental learning, Takagi-Sugeno, neuro-fuzzy

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Classification techniques appear frequently in many application areas, and become the basic tool for almost any pattern recognition task. The main problem in classification is to induce a classifier from a set of data samples. A large amount of samples is needed to set up and evaluate a classification system that can achieve high accuracy, and it is practically very difficult to have such number of samples covering the whole feature space. Therefore, life-long classifier adaptation becomes more and more an essential point. Moreover, in many application contexts, the classifier needs to take into account new unseen classes and to integrate them in the classification process, which increases the need for "evolving" classifiers.

One good example of such applications is the use of online handwriting gesture classifiers that aims at facilitating interactions with computers using pen-based interfaces like whiteboards, tablet PCs, PDA...Etc. The main drawback in the current existing systems is that they are trained "offline" on a specific group of gestures and then implemented to operate without changing their structure during the use. This fixed structure does not allow the user to choose his own set of gestures or to add new ones to assign them to new interactive commands (according to his special needs), for example.

In our work, we aim at building a handwriting classifier, on-the-fly, from scratch and using only few data. Thus, the classifier will be incrementally adapted to achieve high recognition rates as soon as possible and to keep the system robust when introducing new unseen classes at any moment in the life-long learning process.

In order to clarify the meaning of the terms "incremental" and "online" learning, one can compare them to the "batch" and "offline" learning as follows [Kasabov (2007)]:

-Batch versus incremental learning: In a batch mode of learning a predefined training dataset is learned by the system through propagating this dataset several times through the system. Each time the system optimizes its structure based on the average value of the goal function over the whole dataset. The incremental mode of learning is concerned with learning each data sample separately and the data might exist only for a short time. After introducing each data sample, the system makes changes in its structure to optimize the goal function. An incremental learning algorithm must learn new data without fully destroying the knowledge learned from old data and without the need to retrain the system on the old and the new data. -Offline versus online learning: In an offline learning mode, the system is trained on a specific dataset and then implemented to operate in a real environment, without changing its structure during the use. While in an online learning mode, the model learns from new data during its use.

An incremental learning algorithm is defined in [START_REF] Polikar | Learn++: An incremental learning algorithm for supervised neural networks[END_REF]] as being one that meets the following criteria: it should be able to learn additional information from new data; it should not require access to the original data (i.e. data used to train the existing classifier); it should preserve previously acquired knowledge (it should not suffer from catastrophic forgetting, significant loss of original learned knowledge); and it should be able to accommodate new classes that may be introduced with new data. Many of the existing "incremental learning" algorithms are not truly incremental because at least one of the mentioned criteria is violated. These criteria can be briefly expressed by the so called "plasticity-stability dilemma" [START_REF] Zwickel | New Directions in Human Associative Learning[END_REF]], which is that a system must be able to learn in order to adapt to a changing environment but that constant change can lead to an unstable system that can learn new information only by forgetting everything it has so far learned.

We can distinguish two main types of incremental learning algorithms: algorithms for parameter incremental learning and algorithms for structure incremental learning. The incremental learning of parameters can be considered as "adaptation" algorithm. The structure in such systems is fixed and initialized at the beginning of the learning process, and the system parameters are learned incrementally according to newly available data. Some examples of these systems are presented in [Jang (1993)] and [START_REF] Mouchere | On-line writer adaptation for handwriting recognition using fuzzy inference systems[END_REF]].

Most of structure incremental learning algorithms are based on the principle of the ART clustering algorithm [Carpenter and Grossberg (1988)], such as [Carpenter et al. (1992)], [START_REF] Sadri | A new clustering method for improving plasticity and stability in handwritten character recognition systems[END_REF]], [START_REF] Gary | [END_REF]]. The main problem of these systems is that they are sensitive to the selection of the vigilance parameter, to noise level in the training data and to the order in which training data are presented.

In an online incremental learning algorithm, the training set is not available a priori, since the learning examples come over time. Although online learning systems can continuously update and improve their models, not all of them are necessarily based on a real incremental approach. For some systems the model is completely rebuilt at each step of learning using all available data, they are systems with "full instance memory" [START_REF] Reinke | Incremental learning of concept descriptions: A method and experimental results[END_REF]]. On the other hand, if the learning algorithm modifies the model using only the last available learning example, it is called a system with "no instance memory" [START_REF] Littlestone | Redundant noisy attributes, attribute errors, and linearthreshold learning using winnow[END_REF]]. A third category is that of systems with "partial instance memory", which select and maintain a reduced subset of learning examples to use them in the next learning step [START_REF] Aha | Instance-based learning algorithms[END_REF]].

Evolving neuro-fuzzy models [START_REF] Angelov | On-line design of takagi-sugeno models[END_REF]] have been successfully used in many modeling and classification problems during the last years. Several incremental learning algorithms have been proposed for the evolving structure identification and the recursive parameters learning (Denfis [START_REF] Song | Dynamic evolving neuro-fuzzy inference system (denfis): On-line learning and application for time-series prediction[END_REF]], Flexfis [Lughofer (2008)], eTS [Angelov and Filev (2004)], eTS+ [Angelov (2010)] etc.). One of the most used neuro-fuzzy models is the First-order Takagi-Sugeno fuzzy model. It usually has a linear consequent part and its premise part can be learned from data (no expert is needed). These models can address problems with either single output or multi-output. It consists of a set of fuzzy rules of the following form:

Rule i : IF x is close to P i THEN y 1 i = l 1 i (x), ..., y k i = l k i (x) (1)
where l m i (x) is the linear consequent function of the rule i for the class m. The Prototype P is defined by a center and a fuzzy zone of influence. A membership function must be defined in order to calculate the "closeness" degree between x and P (considering its center and its fuzzy zone of influence).

In the existing approaches, the zone of influence has often a hyper-spherical form [START_REF] Yager | Learning of fuzzy rules by mountain clustering[END_REF]] [ Angelov (2004)]. It has then the same radius for all the feature space dimensions, and it can be represented by a diagonal covariance matrix with identical values on the diagonal.

More sophisticated form is presented in [Angelov and Zhou (2006)] that allows the zone radius to have different values for the different dimensions (still diagonal covariance matrix but with different values on the diagonal). This ability results in hyper-elliptical zones where the ellipses are still parallel to the feature space axes.

In this paper, we go a step ahead in the structure of the fuzzy prototypes in first order TS models by allowing them to have a hyper-elliptical form non-parallel to the feature space axes. This form enables the model to take into consideration the correlations that can exist between the features, and is represented by normal (non-diagonal) covariance matrices. We calculate the prototype centers and covariance matrices in a recursive way with zero problem-dependent parameters. We use a membership function based on a well-known multivariate probability distribution function. By experimenting this structure on our handwritten gesture problem and other benchmark multi-class problems, we show its useful effect on the overall system output and the linear consequent learning. We also propose two different methods for incremen-tal adaptation of prototypes, and we compare experimentally these two incremental learning models with a well-known reference model.

The rest of the paper is organized as follows: We describe in Section 2 the architecture of our neuro-fuzzy classifier. Then, the different elements of the used online incremental learning algorithm are detailed in Section III. We present our experimental results in Section IV before concluding in Section 5.

System architecture

As aforementioned, our system is based on first-order Takagi-Sugeno (TS) fuzzy inference system. It consists of a set of fuzzy rules of the following form:

Rule i : IF x is close to P i THEN y 1 i = l 1 i (x), ..., y k i = l k i (x) (2)
where l m i (x) is the linear consequent function of the rule i for the class m:

l m i (x) = π m i x = a m i0 + a m i1 x 1 + a m i2 x 2 + ... + a m in x n ( 3
)
where n is the size of the input vector. Singleton consequences (l m i (x) = a m i ) are sometimes used instead of linear functions in order to get simpler models that are called zero-order TS models. To find the class of x, its membership degree β i (x) to each fuzzy prototype is first computed. After normalizing these membership degrees, the sum-product inference is used to compute the system output for each class:

y m (x) = r i=1 βi (x) l m i (x) (4)
where r is the number of fuzzy rules in the system. The winning class label is given by finding the maximal output and taking the corresponding class label as response:

class(x) = y = argmax y m (x) m = 1, .., k (5)
The membership degree can be calculated in many ways. For hyper-spherical or axes-parallel hyper-elliptical prototypes, the membership degree can be computed depending on the prototype center µ i and the radius value σ i (the same value in all the dimensions for the former, and different values for the later). In this case, the Gaussian membership function is generally used. The value of β i (x) can then be computed as follows:

β i (x) = n j=1 exp(- x -µ i 2 j 2σ 2 ij ) (6)
It must then be normalized as follows:

βi (x) = β i (x) r j=1 β j (x) (7)
In our model, where the fuzzy influence zones of the prototype take a rotated hyper-elliptical form, the membership degree is computed by the prototype center µ i and its variance-covariance matrix A i :

A i =     σ 2 1 c 12 ... c 1n c 21 σ 2 2 ... c 2n ... ... ... ... c n1 c n2 ... σ 2 n     i (8)
where c 12 (= c 21 ) is the covariance of x 1 and x 2 , and so on. In order to calculate the membership degree, we use the multivariate Cauchy probability distribution. The value of β i (x) is computed in our model as follows:

β i (x) = 1 2π |A i | 1 + (x -µ i ) t A -1 i (x -µ i ) -n+1 2 (9)
The difference between the two types of premises is illustrated in fig. 1 (data from "Iris" dataset [START_REF] Frank | UCI machine learning repository[END_REF]]).

On-line incremental learning algorithm

The incremental learning algorithm of our model consists of three different tasks: the creation of new rules, the adaptation of the existing rule's premises, and the tuning of the linear consequent parameters. These three tasks must be done in an online incremental mode and all the needed calculation must be completely recursive. For the incremental rule creation, we adapt a technique of incremental clustering. Coupling this technique with the incremental adaptation of the premises of the rules is part of the originality of our learning model. We propose two different methods for adapting the premises, and we use a modified version of the recursive least squares method for estimating the linear consequent parameters in incremental manner.We will describe in this section the different parts of the learning algorithm.

Incremental clustering

In an online incremental learning problem, training data become available continuously, and the system must be learned using the first-arrived data, and then continue to evolve in a transparent manner from the user viewpoint. Classic clustering algorithms need to know all the objects in order to perform the clustering and each time we modify the set of objects, it is necessary to cluster the whole collection again. As mentioned in the introduction, a very important criterion of an efficient incremental learning solution is avoiding (or minimizing) access to previous learning data. Thus, we need algorithm able to update the clusters each time a new data become available, without neither rebuilding the whole set of clusters, nor requiring access to previous data. When introducing a new training sample in an online learning mode, it will either reinforce the information contained in the previous data and represented by the current clustering, or bring enough information to form a new cluster or modify an existing one. The importance of a given sample in the clustering process can be evaluated by its potential value. The potential of a sample is defined as inverse of the sum of dis-tances between a data sample and all the other data samples [START_REF] Yager | Learning of fuzzy rules by mountain clustering[END_REF]]:

P ot(x(t)) = 1 1 + t-1 i=1 x(t) -x(i) 2 (10)
A recursive method for the calculation of the potential of a new sample was introduced in [Angelov and Filev (2004)], which made this technique a promised solution for any incremental clustering problem. The recursive formula avoids memorizing the whole previous data but keeps -using few variables -the density distribution in the feature space based on the previous data:

P ot(x(t)) = t -1 (t -1)α(t) + γ(t) -2ζ(t) + t -1 (11)
where

α(t) = n j=1 x 2 j (t) (12)
γ(t) = γ(t -1) + α(t -1), γ(1) = 0 (13) ζ(t) = n j=1 x j (t)η j (t), η j (t) = η j (t -1) + x j (t -1), η j (1) = 0 (14)
Introducing a new sample affects the potential values of the centers of the existing clusters, which can be recursively updated by the following equation:

P ot(µ i ) = (t -1)P ot(µ i ) t -2 + P ot(µ i ) + P ot(µ i ) n j=1 µ i -x(t -1) 2 j (15)
If the potential of the new sample is higher than the potential of the existing centers then this sample will be a center of a new cluster and a new fuzzy rule will be formed in the case of our neuro-fuzzy model. So, the center of the new prototype µ r+1 = x k and its covariance matrix A r+1 = ǫI, where I is the identity matrix of size n and ǫ is a problem-independent parameter and can generally be set to 10 -2 .

Premise adaptation

This adaptation process allows to incrementally update the prototype centers coordinations according to each new available learning data, and to recursively compute the prototype covariance matrices in order to give them the rotated hyper-elliptical form. For each new sample x k , the center and the covariance matrix of the prototype that has the highest activation degree are updated. We present two different approaches of adaptation; the first method is purely statistical, while the second takes into account the error or the confusion during the recognition of the new example.

Statistical recursive adaptation

In this method, the center coordination of the selected prototype is recalculated as follows:

µ i = s i -1 s i µ i + 1 s i x t ( 16
)
where s i represents the number of updates that have been applied so far on this prototype. The covariance matrix is recursively computed as follows:

A i = s i -1 s i A i + 1 s i -1 (x t -µ i )(x t -µ i ) T (17)
For practical issues, since the membership degree can be calculated using only

A -1 (|A| = 1 |A -1 | )
, and in order to avoid any matrix inversion, we use an updating rule for A -1 directly [De Backer and Scheunders ( 2001)]:

A -1 i = A -1 i 1 -α - α 1 -α • (A -1 i (x t -µ i )) • (A -1 i (x t -µ i )) T 1 + α((x t -µ i ) T A -1 i (x t -µ i )) (18)
where a = 1 si-1 .

Confusion-driven recursive adaptation

We propose a second method of adaptation in which the adjustment is made in "supervised" manner. To guide the adaptation of premises to improve system performance, we focus on examples that are either misrecognized, or "hardly" well recognized. This is translated by adding in the adaptation formulas a weight w that will be calculated for each new instance and it represents the importance to be assigned to this example in the process of adaption. The adaptation formulas are then given as follows:

µ i = s i -w s i µ i + w s i x t ( 19
)
A i = s i -w s i A i + w s i -1 (x t -µ i )(x t -µ i ) T (20) w = 1 -[y c -y nc ] w ∈ [0, 2] (21)
where

y nc = argmax y m m = 1..k & m = c (22)
In words, w is inversely proportional to the difference between the true label score of x t (y c ) and the highest value between the other (wrong) classes' scores (y nc ). Thus, the value of w tends to 0 when x t is "strongly" well recognized, and to 2 when it is misrecognized.

Linear consequent tuning

The tuning of the linear consequent parameters in a first-order TS model can be done using weighted Recursive Least Square method (wRLS) [Angelov et al. (2008)]. The wRLS method can be applied in either local or global manner. Consequent parameters of each rule are optimized separately in the former, while parameters of all the rules are optimized together in the latter. Let Π be the matrix of all the linear consequents parameters in first-order FIS :

Π =     π 1 1 π 2 1 ... π k 1 π 1 2 π 2 2 ... π k 2 ... ... ... ... π 1 r π 2 r ... π k r     ( 23
)
where m is the number of classes, and r is the number of fuzzy rules; It can be globally and incrementally estimated by:

Π t = Π t-1 + C t ψ t (Y t -ψ T t Π t-1 ) ; Π 1 = 0 (24) C t = C t-1 - C t-1 ψ t ψ T t C t-1 1 + ψ T t C t-1 ψ t ; C 1 = ΩI ( 25
)
where

ψ t = [β 1 (x t )x t , β 2 (x t )x t , ..., β r (x t )x t ]
is the input vector weighted by the activation levels of the prototypes, I is the identity matrix and Ω is a large positive number, typically between 100 and 10,000. Small values of Ω slow down the learning, while too large Ω can prevent parameters converging properly. The value must then be estimated to be a good compromise between these two points. Ω = 1000 is adequate for most cases. Y t is a 1-by-k target vector of the input vector x t . In classical classification problems, the value of the correct class in the target vector is set to one, and the rest to zero.

When a new rule is added to the system, its consequent parameters are initialized by the weighted average of the parameters of the other rules, while the parameters of the other rules do not change:

Π t =        π 1 1(t-1) π 2 1(t-1) ... π k 1(t-1) π 1 2(t-1) π 2 2(t-1) ... π k 2(t-1) ... ... ... ... π 1 r(t-1) π 2 r(t-1) ... π k r(t-1) π 1 (r+1)t π 2 (r+1)t ... π k (r+1)t        (26)
where

π c (r+1)t = r i=1 β i (x t )π c i(t-1) (27)
In addition, the covariance matrix C is extended as follows:

C t =       ρ C t-1 0 0   Ω ... 0 ... ... ... 0 ... Ω         (28)
where ρ = (r 2 + 1)/r 2 . The complete learning algorithm can be summarized by Algorithm 1.

Algorithm 1: First-order Takagi-Sugeno online incremental learning algorithm foreach new sample x do if x is the first sample of a new class then create a new fuzzy prototype based at x; initialize its potential by 1; add a new fuzzy rule to the system; extend the consequent parameters matrix as in ( 26) and ( 27); update and extend the covariance matrix as in (28); else calculate the activations of the fuzzy rules by ( 9); determine the winning class label by [5]; get the true class label; calculate the potential of x by (11); update the potentials of the existing prototypes centers using (15);

if P ot(x) > P ot(µ i ) ∀i ∈ [1, R] then if
x is close to an existing center µ i then let x be the new center of the prototype P i and keep the same consequents of the rule i ; else create a new fuzzy prototype based at x; initialize its potential by 1; add a new fuzzy rule to the system; extend the consequent parameters matrix as in ( 26) and ( 27); update and extend the covariance matrix as in (28); end end Apply premise adaptation according to x b by either ( 16) and ( 18) or ( 19) and ( 20); update the consequents parameters using ( 24) and ( 25); end end Fig. 2 Handwritten gestures in the dataset SIGN

Experimental results

In this section, we analyze the performance and behavior of our model of incremental learning on two types of experience. The first experiment is conducted on an incremental learning problem related to our specific application context, online handwritten gesture recognition. In the second part of this experimental study, we test our model on well-known classification benchmarks, allowing comparing our results with non-incremental learning methods, while offering the opportunity to other incremental approaches to compare with our method.

Online handwritten gesture recognition

We led the experiments on the "SIGN" database, which is a database of on-line handwritten gestures (figure 2). It is composed of 25 different gestures drawn by 11 different writers on Tablet PCs. Each writer has drawn 100 samples of each gesture, i.e. 2,500 gestures in each writer-specific dataset. The dataset (and additional information on the data collection protocol) can be found in [ SIGN (2010)]. Each gesture is described by a set of 10 features. The presented results are the average of results of 11 different tests for the 11 writers.

UCI repository datasets

Besides the SIGN dataset, we evaluate the presented neuro-fuzzy model on some benchmark problems form the UCI machine learning repository [START_REF] Frank | UCI machine learning repository[END_REF]]. We particularly focus in this work on multi-class problems. We have chosen the next datasets:

-CoverType: The aim of this problem is to predict forest cover type from 12 cartographical variables. Seven classes of forest cover types are considered in this dataset. We use in our experiment a subset of 2100 instances. -PenDigits: The objective is to classify the ten digits represented by their handwriting information from pressure sensitive tablet PC. Each digits is represented by 16 features. The dataset contains about 11000 instances. We use 75% of each dataset for the incremental learning process and the rest is used to estimate the performance during and at the end of the learning process. In order to get the results unbiased by the data order effect, we repeat the experiment for each dataset 40 times with different random data orders and the mean results and standard deviations are presented in the figures and the tables. Three incremental learning models are compared in this experiment:

(I) Euclidean + statistical adaptation: an evolving first-order TS classifier with parallel-to-axes hyper-elliptical prototypes (this model is equivalent to eClass model [Angelov et al. (2008)]), (II) Mahalanobis + statistical adaptation: our extended version with rotated hyperelliptical prototypes (covariances between features are considered). The adaptation of prototypes is purely statistical (section 3.2.1) (III) Mahalanobis + confusion-driven adaptation: prototypes adaptation is based on the degrees of confusion (section 3.2.2) in this model.

As for the premise structure, several variants of consequent functions and their learning algorithm can be used as aforementioned. We have measured the impact of these variants on the classification accuracy using the enhanced premise structure on the five datasets. Table 1 shows a comparison between zero-order and first-order model and between local consequent optimization and global consequent optimization, for a hyper-elliptical premise structure with statistical adaptation. We can note that first-order models are more performant than zero-order ones, but with the price of having more parameters. On the other hand, we note that both local and global consequents optimization give almost the same results, with the advantage of having less parameter using the local optimization. Thus, we choose in the next experiment a local consequent optimization algorithm, and we compare the three premise structures that we have mentioned above. Results are shown for both zero-order and first-order models. We note from Table 2 that our improved models outperform the reference model. With the improved structure of premises, the rate of misclassification generally decreases from 10% to 50%, according to the classification problem. We can also notice that the confusion-driven premise adaptation helps in correcting some classification errors and improves the system's performance compared to the statistical adaptation. Almost the same results are obtained for either zero-order or first-order consequents.

We show in figure 3, 4, 5 and 6 the evolution of classification performance during the incremental learning process using the three different premise structures. We show in figure the evolution of classification performance during the incremental learning process using the three different premise structures. We notice the stability and the superiority of the performance using the improved premise structure.

Conclusion

An improvement on the premise structure of Takagi-Sugeno neuro-fuzzy classifier had been presented in this paper. The rotated fuzzy zones of influence in the presented model enable the premise "layer" to consider the correlations that may exist between input features and boost the overall performance of the classifier. An appropriate recursive premise adaptation methods are coupled with a known incremental clustering technique and used with the recursive least squares method to learn the classifier in one-pass incremental mode. The proposed enhancements reduce the misclassification rate by more than 30% for our specific classification problem (online handwritten gesture recognition), and these good results are also obtained on several benchmark classification problems. on Fuzzy sets and systems. pp. 576-584. [Angelov and Filev (2004)] Angelov, P., Filev, D., Feb. 2004. An approach to online identification of takagi-sugeno fuzzy models. IEEE Transactions on Systems, Man, and Cybernetics 34 (1), 484-498. [Angelov et al. (2008)] Angelov, P., Lughofer, E., Zhou, X., 2008. Evolving fuzzy classifiers using different model architectures. Fuzzy Sets Syst. 159 (23), 3160-3182. [Angelov and Zhou (2006)] Angelov, P., Zhou, X., 2006. Evolving fuzzy systems from data streams in real-time. In: IEEE Symposium on Evolving Fuzzy Systems. [Carpenter et al. (1992)] Carpenter, G., Grossberg, S., Markuzon, N., Reynolds, J., Rosen, D., 1992. Fuzzy artmap: A neural network architecture for incremental supervised learning of analog multidimensional maps. IEEE Transactions on Neural Networks 3. [Carpenter and Grossberg (1988)] Carpenter, G. A., Grossberg, S., 1988. The art of adaptive pattern recognition by a self-organizing neural network. Computer 21 (3), 77-88.

Fig. 1

 1 Fig.1The zone of influence of a prototype is parallel to the axes in (a), (c) and (e), while the rotated zones in (b), (d) and (f) through the use a multidimensional probability distribution result in a more accurate data modeling.

-

  Segment: Each instance in the dataset represents a 3x3 region from 7 outdoor images. The aim is to find the image from which the region was taken. Each region is characterized by 19 numerical attributes. There are 2310 instances in the dataset.-Letters: The objective is to identify each of a large number of black-and-white rectangular pixel displays as one of the 26 capital letters in the English alphabet. Each letter is represented by 16 primitive numerical attributes. The dataset contains 20000 instances.

Fig. 3

 3 Fig. 3 Evolution of performance during the incremental learning process(SIGN)

Table 1

 1 Misclassification rates for different consequent types

Table 2

 2 Misclassification rates for different premise structures


E Andreotti

C Arnaboldi

F T Avignoneiii

M Balata

I Bandac

M Barucci

J W Beeman

F Bellini

C Brofferio

A Bryant

C Bucci

L Canonica

S Capelli

L Carbone

M Carrettoni

M Clemenza

O Cremonesi

R J Creswick

S Di Domizio

M J Dolinski

L Ejzak

R Faccini

H A Farach

E Ferri

E Fiorini

L Foggetta

A Giachero

L Gironi

A Giuliani

P Gorla

E Guardincerri

T D Gutierrez

E E Haller

K Kazkaz

S Kraft

L Kogler

C Maiano

R H Maruyama

C Martinez

M Martinez

S Newman

S Nisi

C Nones

E B Norman

A Nucciotti

F Orio

M Pallavicini

V Palmieri

L Pattavina

M Pavan

M Pedretti

G Pessina

S Pirro

E Previtali

L Risegari

C Rosenfeld

C Rusconi

C Salvioni

S Sangiorgio

D Schaeffer

N D Scielzo

M Sisti

A R Smith

C Tomei

G Ventura

M Vignati

F T Avignone Iii

Te Neutrinoless Double-Beta Decay with CUORICINO

We report the final result of the CUORICINO experiment. Operated between 2003 and 2008, with a total exposure of 19.75 kg•y of 130 Te, CUORICINO was able to set a lower bound on the 130 Te 0νββ half-life of 2.8 × 10 24 years at 90% C.L. The limit here reported includes the effects of systematic uncertainties that are examined in detail in the paper. The corresponding upper bound on the neutrino Majorana mass is in the range 300-710 meV, depending on the adopted nuclear matrix element evaluation.

Introduction

Double beta decay (ββ) is a rare transition between two isobars, involving a change of the nuclear charge Z by two units. There exist several naturally-occurring even-even nuclei for which this is the only allowed decay mode. While the transition involving 2 electrons and 2 (anti)neutrinos (2νββ) is allowed in any theoretical model (and has been observed for various isotopes), this is not the case for the neutrinoless channel (0νββ). Despite being energetically possible, 0νββ violates lepton number conservation by 2 units and is possible only if the neutrino is a massive Majorana particle [START_REF] Tretyak | Examples of reviews are[END_REF]. 0νββ searches have been pursued for more than half a century and today they experience a renewed interest, thanks to the discovery of neutrino oscillations [2,3,4]. Neutrino oscillations imply that neutrinos have a finite mass, but are not enough to determine the nature of that mass (Dirac or Majorana) or why it is so extraordinarily small. Several theoretical speculations point toward a mass generation mechanism that implies a Majorana character of neutrinos, and that indicates the 0νββ process as the unique tool with a discov-ery potential [START_REF] Tretyak | Examples of reviews are[END_REF]2].

The 0νββ transition could proceed through several different mechanisms, of which the simplest and most commonly cited is the exchange of a light Majorana neutrino. In this case the observation of 0νββ would not only provide evidence for lepton number violation and the Majorana character of the neutrino, but would also result in a measurement of the effective Majorana mass m ee = | m i U 2 ei | (where m i are the three mass eigenstates of neutrinos and U ei are the PMNS matrix elements) through the relation:

1 τ 0ν 1/2 = m 2 ee F N = m 2 ee G 0ν |M 0ν | 2 (1)
Here τ 0ν 1/2 is the 0νββ half-life, G 0ν is the two-body phase-space factor and M 0ν is the 0νββ nuclear matrix element (NME). The product F N = G 0ν |M 0ν | 2 is called the "nuclear factor of merit"; the name refers to the fact that, according to equation [START_REF] Tretyak | Examples of reviews are[END_REF], F N directly influences the experimental sensitivity to m ee . The main uncertainty in deriving m ee (or an upper limit on it) from the Figure 1: Anticoincidence total energy spectrum of all CUORICINO detectors (black). The most prominent peaks are labeled and come from known radioactive sources such as: e + e -annihilation (1), 214 Bi (2), 40 K (3), 208 Tl (4), 60 Co (5) and 228 Ac (6). The total energy spectrum of all CUORICINO detectors during calibration measurements is also shown (color). For convenience, it is normalized to have the same intensity of the 2615 keV line of 208 Tl as measured in the non-calibration spectrum.

experimental result on τ 0ν 1/2 comes from the NME, which is a theoretical calculation still affected by a large spread among the adopted nuclear models and their implementations [5,6,7,8]. To mitigate this uncertainty, several candidate 0νββ isotopes could be studied and the m ee results compared. Although it may not be feasible to study all the ββ candidates with reasonable sensitivity, there exists a "golden list" of isotopes that -as a compromise of cost, availability, technological approach and other factors -have been studied so far. Among them, those that have yielded the most stringent limits on m ee (within the NME spread) are 76 Ge [9,10], 100 Mo [11], 130 Te [START_REF] Arnaboldi | [END_REF] and 136 Xe [13]. In all but one [14] case, only upper bounds on the Majorana mass have been reported.

In this paper, we discuss the final 0νββ result of CUORI-CINO, which yields the most stringent bound on m ee based on 130 Te studies, and one of the best in general. CUORICINO data acquisition started in April 2003 and ended in June 2008. The data are separated into two runs (RUN I and RUN II), due to a major maintenance interruption. The data collection is summarized in table 1. A partial data-set of 11.83 kg•y of 130 Te exposure was used for the analysis reported in [START_REF] Arnaboldi | [END_REF].

The paper is organized as follows: after a short description of the experimental set-up in section 2, we present details of RUN II data analysis, discussing processing in section 3, data reduction in section 4 and efficiency evaluation in section 5. In section 6, we describe two Bayesian approaches used for the 0νββ half-life limit evaluation (one of them is that adopted in [START_REF] Arnaboldi | [END_REF]), testing the two procedures on a toy Monte Carlo and discussing their compatibility on real data. In section 7 we discuss the influence of systematic errors on the final result.

We conclude the paper with the CUORICINO final result for the 0νββ half-life limit evaluated on the entire data set. This is done treating RUN I and RUN II as two independent experiments whose likelihoods are combined. This choice was motivated by the difference in detector configuration between the two runs (increased number of active detectors, improved performance) and a presumable difference in background composition (due to detector exposure to air).

CUORICINO

For many years, the most sensitive 0νββ results for 130 Te have been obtained using bolometric detectors. A bolometer is a type of calorimeter operated at ultra-low temperature, in which the energy of incident radiation is converted to heat, raising the temperature of the detector's body. The energy released in the detector is then determined by measuring the temperature increase, in our case by using a semiconductor thermistor [START_REF] Haller | Neutron transmutation doping of semiconducting materials[END_REF] whose resistance varies exponentially with temperature. In this kind of detector, the candidate isotope is contained within the active mass of the detector itself. CUORICINO bolometers contain the isotope 130 Te (isotopic abundance 33.8%) which is a ββ candidate with a rather favorable factor of merit 1 .

About 85% of the time (see section 4), the two electrons emitted by 130 Te 0νββ would be fully contained within one crystal. The signature of the decay would therefore consist of a monochromatic peak in the energy spectra of the bolometers at an energy equal to the Q-value of the decay: 2527.518±0.013 keV [START_REF] Redshaw | [END_REF]. The difficulty of the experiment lies in the control and reduction of all the background events that could mimic such a signal. These can be non-particle signals, due to electronic or thermal noise, or particle signals, due to radioactivity and cosmic rays. The former are rejected on the basis of a pulse shape discrimination technique (see section 4). The latter are controlled during the experiment's design and construction by proper material selection, handling, and shielding [START_REF] Arnaboldi | [END_REF]18,19], and -at the stage of data-analysis -by coincidence cuts (see section 4). CUORICINO was the latest step in a long series of experiments of increasing mass, performance and sensitivity [20,21,22]. The next experiment to use this technique will be CUORE [23], which is presently under construction.

CUORICINO [START_REF] Arnaboldi | [END_REF] was a tower array of 62 TeO 2 crystals used as bolometric detectors. The array was operated underground, in a dilution refrigerator located at the Laboratori Nazionali del Gran Sasso (INFN -Italy), which provides an average coverage of 1400 m of rock (3650 m.w.e.). CUORICINO crystals can be divided into four main groups according to their mass and isotopic abundance. These are:

• the "big crystals" -44 bolometers, 5 × 5 × 5 cm 3 in size and 790 g in mass;

• the "small crystals" -14 bolometers, 3 × 3 × 6 cm 3 in size and 330 g in mass;

• the " 130 Te-enriched crystals" -2 bolometers, 3×3×6 cm 3 in size and 330 g in mass, grown with 130 Te enriched material;

• the " 128 Te-enriched crystals" -2 bolometers, with the same size and mass of the 130 Te-enriched ones but grown with 128 Te enriched material.

The 130 Te-enriched crystals have a 130 Te content corresponding to an isotopic abundance of 75% [22], while the 130 Te content of 128-enriched crystals is so low that they will not be considered for the 0νββ analysis presented in this paper.

A detailed description of the array, cryogenic set-up, shields, front-end electronics, and DAQ can be found in [START_REF] Arnaboldi | [END_REF] and references therein. Here we will describe the main steps of data acquisition and data handling which are relevant for the discussion.

The output voltage of each detector was monitored by a constant fraction trigger. When the output voltage exceeded the trigger threshold, the acquisition system recorded 512 samples (a 4 s window sampled at 125 Hz), which constitute one "event". The acquired time window fully contained the pulse development, providing an accurate description of its waveform. A pre-trigger interval just prior to the onset of the pulse was used to measure the DC level of the detector, which corresponds to the instantaneous detector temperature. The pulse amplitudes were evaluated offline for each recorded event, together with a few other characterizing parameters of the pulse. Each single CUORICINO measurement lasted about 22 hours on average, with the time between measurements (about 2 hours) dedicated to cryogenic system maintenance. A routine calibration with an external 232 Th source was performed approximately once per month, lasting for about 3 days. The accumulated data between two calibrations is referred to as a "data-set". The spectrum obtained by summing all the CUORICINO collected data (i.e. summing over detectors and data-sets) is shown in figure 1. The background recorded by the detectors is clearly dominated, in this region, by gamma emissions due to radioactive contaminations of the detector and of the surrounding apparatus. The most intense gamma lines are listed in reference [START_REF] Arnaboldi | [END_REF]. In figure 1, the spectrum corresponding to the sum of all calibration data is also shown. For convenience the calibration spectrum is normalized to have the same intensity of the 2615 keV line of 208 Tl as measured in the background spectrum.

Data processing

The analysis of CUORICINO data starts with the collection of all the triggered events. For clarity, we will model the single waveform V(t) induced by a particle interaction in the crystal as:

V(t) = V 0 s(t) + n(t) (2) V 0 = G(T ) • A(E) (3)
In the first equation, V 0 is the maximum value of the raw signal acquired at time t 0 , s(t) describes the shape of the particle signal and n(t) is an additive noise source. The second equation describes the dependence of the signal amplitude on the detector working temperature. Here we assume that the dependence of the gain on temperature, G(T ), and that of the amplitude on energy, A(E), can be factorized. This is not true in general, however it is a good approximation when dealing -as in our case -with small temperature drifts. Although it describes a naive model, this formula highlights the key points of the analysis. In order to estimate E, we need:

1. a technique able to measure V 0 , reducing the effect of n(t)

as much as possible, in order to improve our resolution (amplitude evaluation); 2. an algorithm to control for the variation of G(T ) produced by detector temperature drifts (gain instability correction); 3. a technique to measure the form of A(E) (energy calibration).

Amplitude evaluation. This is done by maximizing the signal-to-noise ratio by means of optimum filtering [24]: each waveform is convolved with a transfer function h(t) whose Fourier transform is defined as:

H(ω) = e iωt max S * (ω) N(ω) ( 4
)
where S (ω) is the Fourier transform of the average detector response function s(t), N(ω) is the spectral power density of the noise characterizing the detector, and t max is the time at which the pulse reaches its maximum. The functions s(t) and N(ω) are computed by an averaging procedure of the bolometric pulses and of the Fourier transformed baselines 2 . Figure 2 shows an example of an event due to a particle interaction and an example of a non-particle event, most likely due to an abrupt temperature increase produced by an electric disturbance or by vibrations. Each waveform is superimposed with its optimum-filtered counterpart. Once the optimum filter is applied, the amplitude of the signal is inferred from the maximum of the filtered waveform in time domain.

Gain instability correction. This correction is achieved by measuring the voltage amplitude, V re f , of a monochromatic reference pulse. This pulse is produced by depositing a fixed amount of energy into the crystal by the Joule dissipation from a heavily doped silicon resistor glued to the crystal. Because the energy deposited is fixed, any variation of V re f would be due to a variation in G(T ), which can be measured and used to correct the amplitudes of all the triggered events. For a more detailed discussion of this method we refer to [25,26].

Energy calibration. The voltage-energy relationship is reconstructed by means of routine source calibrations: two wires of thoriated tungsten are periodically inserted between the cryostat and its external lead shield. The voltage amplitude of the pulses corresponding to the main gamma lines of 232 Th are used for the determination of the parameters describing the A(E) relationship. This function is characterized by different non-linearity sources [START_REF] Arnaboldi | [END_REF], the dominant one being the dependence of the thermistor resistance on the temperature [27]. In this work, A(E) is parametrized with a third-order polynomial, which can be considered as the truncated Taylor's expansion of the real unknown calibration function. In the previous CUORI-CINO analysis [START_REF] Arnaboldi | [END_REF], a different calibration function was used. This was a second order polynomial in log(V) and log(E), based on a thermal model describing our bolometric detectors. While this function performs better at extrapolation (i.e. above the highest calibration line at 2615 keV), the third order polynomial performs better in the interpolation region (i.e. between threshold and 2.6 MeV). The difference between the two parameterizations was studied using the total background spectrum recorded by CUORICINO (figure 1); this spectrum contains several gamma peaks whose origins, and therefore nominal energies, are clearly identified. The difference between the nominal energy of each peak and its measured position (the residual) is plotted against the nominal energy in figure 3, showing the slightly better performance of the 3rd order polynomial. These residuals also provide important information concerning the precision of our calibration: their spread can be used as an estimator of the uncertainty in the energy position of a peak, including that produced by the 0νββ signal.

Source calibration measurements are repeated for each dataset and are also used to check the detector performances over time. Figure 4 shows the distribution of all the resolutions measured in calibration for the three crystal groups (big, small, and enriched crystals).

Data reduction

The final CUORICINO spectrum is composed of events which survived two different types of data selection: global and event-based cuts.

Global cuts: these are applied following quality criteria decided a priori (e.g. an excessive noise level or an incompatibility between the two calibration measurements at the beginning and the end of a data-set). They identify bad time intervals to be discarded. This kind of cut introduces a dead time that is accounted for by properly reducing the live time of the detector of interest. The cuts are generally based on off-line checks that monitor the detector performances and flag excessive deviations from global control quantities (average resolution, average rate, etc.). The total dead time introduced by these global cuts is ∼5%. A further dead time is introduced by the rejection of a short time window centered around each reference pulse (the frequency with which the reference pulses are generated is about 3 mHz). This cut ensures the rejection of possible pile-up of a particle signal with the reference pulse (the impact of this cut is reported as an efficiency in table 2).

Event-based cuts: these are the pulse-shape and the anticoincidence cuts. The former is used to reject non-physical and pile-up events (the presence of a pile-up prevents the optimum filter algorithm from providing a correct evaluation of the pulse amplitude). The latter allows for the reduction of the background counting rate in the region of interest (ROI). The 0νββ signature we look for consists of a "single-hit" event (only one detector at a time involved), while many of the background counts in the ROI are due to "multiple-hit" events. These include events due to alpha decays on the crystal surfaces that deposit energy in two neighboring crystals and events due to gammas that Compton scatter in one crystal before interacting in another one.

The pulse shape parameters used in this analysis are the rise time and decay time of the raw waveform and a parameter that measures the consistency of one of the basic statements of optimum filter theory. This "Optimum Filter Test" parameter (OFT) is the difference (expressed in percentage of the total amplitude) between the evaluation of the pulse height in the time domain (as the maximum value of the filtered pulse) and that in the frequency domain (as the integral of the filtered-pulse power spectrum). Indeed, if the shape of an event is identical to the average detector response, the two methods yield the same result. However, if the shape of the signal is different from the expected one (such as in the case of a non-physical or a pile-up event), they differ. Figure 5 shows the scatter plot of OFT as a function of energy for a CUORICINO detector (here only one data-set is reported). The main trend reflects the change of the signal shape with energy, and has a minimum in the region where the average response was measured (1-2 MeV). This variation in the pulse shape is caused by non-linearities introduced by the thermistor. Outliers on this plot correspond to misshapen events which will be discarded; the colored vertical bars identify confidence regions evaluated automatically on this distribution in order to obtain cuts which are independent of the signal amplitude.

The anticoincidence cuts require that only one detector fires within a time window of 100 ms.

Signal efficiency

The signal efficiency is the probability that a 0νββ event is detected, its energy is reconstructed accurately, and that it passes the data selection cuts. This parameter must be accurately determined, since it is used to obtain the number of 130 Te 0νββ events. The overall signal efficiency is: (82.8 ± 1.1)% for the CUORICINO big crystals and (79.7 ± 1.4)% for the small and the 130 Te-enriched ones. These efficiencies were computed as discussed below.

There are two main sources of inefficiencies, one "physical" that can be computed by simulations, and the other "instrumental" that must be measured from the data. The mechanism of "physical" efficiency loss is the escape of a fraction of the 0νββ energy from the source crystal. Mechanisms for the "instrumental" efficiency loss are: the pulse-shape cut, the anti-coincidence cut and an incorrect assignment of the The main trend identifies "good" events. Pile-up and non-physical pulses are outside the confidence regions (identified by the colored vertical bars). At low energy, there is a high density distribution of events for OFT values between 10 -2 and 10 -1 ; these events are due to electric disturbances. Physical efficiency: the 0νββ signature is a sharp peak centered at the transition energy (Q-value) of the decay. The peak is produced by 0νββ decays fully contained within the source crystal. The containment probability was evaluated using a Geant4-based Monte Carlo simulation that takes into account all the possible energy escape mechanisms (i.e. electrons, Xrays or bremsstrahlung photons escaping from the source crystals). Since the escape probability depends on the crystal geometry, the efficiency is slightly different for the big and the small crystals (see table 2).

Instrumental efficiency: this is the product of the pulseshape cut, anti-coincidence and excess noise efficiencies. To evaluate the efficiency of the pulse shape cut, the background photopeak at 2615 keV due to 208 Tl was used as a proxy for the 0νββ peak. The 2615 keV peak was chosen because of its proximity to the 0νββ energy and its relatively high intensity. In principle, an ideal pulse-shape cut should leave the main peak untouched and should only reduce the flat background. The area of the peak can then be computed in terms of the the total number of signal events (N sig ), the signal efficiency ( PS ), the total number of background events (N bkg ) and the background efficiency ( bkg ). A simultaneous fit was done on both the spectra of accepted and rejected events. The area of the peak in the accepted events spectrum is given by PS × N sig , while for the rejected events it is (1 -PS ) × N sig . Similarly, the background yield for the accepted events is bkg × N bkg , and the background yield for the rejected events is (1bkg ) × N bkg . By including PS directly in the parametrization of the fit, correlations among the fit parameters are automatically taken into account when the error on PS is calculated. The result is PS = (98.5 ± 0.3)%, and bkg = (64 ± 2)%. The pulse shape cut is clearly very powerful, rejecting approximately 36% of the events in the continuum background while retaining 98.5% of the signal events in the peak. The events discarded in this region are mainly pile-up with real or spurious signals.

To estimate the efficiency of the anti-coincidence cut, we used the same procedure but considered the only available high intensity peak that is produced by a nuclear decay with no detectable coincidence radiation. This is the 1460 keV gamma line emitted in 40 K electron capture (the only coincident radiation, a 3 keV X-ray, is far below our threshold). Since events in this photopeak are single hits, their reduction after an anticoincidence cut can be ascribed only to random coincidences.

The last source of inefficiency is the loss of 0νββ events due to excess noise which can distort the pulse shape and introduce an error in the reconstructed energy. If such an error is greater than the resolution, the event can be considered as lost in the continuum background. In order to estimate this efficiency we compare the number of reference pulses generated during the measurements (the signals used for the gain instability correction, see section 3) with the number actually measured in the correct energy range.

0νββ analysis

The definition of the energy window used to fit the 0νββ spectrum, the hypothesis assumed for the background shape and the number of free parameters used to describe the background itself are extremely important for the choice of the analysis procedure and for the determination of its systematics. The choice of the energy window is somewhat arbitrary, but it influences the background representation. If the energy window is too wide (compared to the signal FWHM) a very precise knowledge of the background shape is necessary. Obviously there is also a minimum width necessary to be able to evaluate the background level beyond the 0νββ peak. In our case, there is a background line near the 0νββ energy, at ∼2505 keV, due to 60 Co (sum of the two photons emitted in cascade by 60 Co decay), which should also be included in the window. Given these considerations, our final choice for the fit window is 2474-2580 Figure 6: A closer view of the CUORICINO anticoincidence spectrum (presented in figure 1) near the 0νββ ROI. Note that the efficiencies listed in table 2 are not yet included.

keV. This is the widest window centered on the ββ Q-value that allows the following two background peaks to be excluded from the fit: the 2448 keV line of 214 Bi and the 2587 keV Te X-ray escape peak of the 208 Tl line. The latter peak is clearly visible in the CUORICINO calibration spectrum shown in figure 6, although -due to the lower statistics -it is not visible in the background spectrum shown in the same figure.

Within this picture, we choose the simplest possible model for the shape f i, j (E) of the spectrum (normalized to mass, live time, efficiency and isotopic abundance) of each single detector (index i) in each data-set (index j) as:

f i, j (E) = B i + Γ 60 Co i, j g i, j E -E Co + Γ 0ν g i, j E -E 0ν (5)
Here g i, j (E) is the function describing the shape of monochromatic energy lines in the i th detector, during the j th data-set, i.e. the response function that is represented by a gaussian with a width determined from calibration data 3 . E 0ν is the 130 Te ββ Qvalue, fixed at its measured value (2527.5 keV). E Co is the sum energy of the two 60 Co gamma lines (2505.7 keV). B i is the flat background component for the i th detector (here we assume a time independent background). Finally Γ 60 Co i, j and Γ 0νββ are, respectively, the 60 Co activity for the i th detector during the j th data-set ( 60 Co has a half-life of 5.27 years), and the absolute activity for 0νββ , both expressed in counts/kg/y.

Free parameters are B i , Γ 60 Co i, j , Γ 0ν and E Co , the parameters of the response function being fixed at the values measured during calibrations. Note that the dependence of Γ 60 Co i, j on the index j is determined by the 60 Co half-life; therefore, the total number of free parameters is determined only by the number of detectors (i.e. by the index i).

Statistical approaches

The CUORICINO spectrum shows no evidence of a 0νββ signal, thus we will provide a limit for the half life of 130 Te by means of a Bayesian approach. Unlike other low statistics methods such as that of Feldman and Cousins [28], this technique does not require an exact evaluation of the expected number of background events (which is unknown). In our case, all the uncertainties are marginalized in the process of the limit computation, and our prior knowledge for the rate Γ 0ν will be represented by a flat distribution, excluding the non-physical region.

Once the statistical method is chosen, we need to decide how to model the experiment: every CUORICINO detector can in fact be imagined as an independent search for 0νββ, characterized by its own background and resolution. There are three natural approaches which can be chosen to search for a 0νββ signal: I -treat the different detectors separately; II -treat the different detectors separately, assuming an identical background for all detectors within each group (big, small, 130 Te-enriched); III -sum the spectra of all detectors belonging to the same group.

In approach I, each detector and each data-set is fit with its own function f i, j (E). In principle this is the best approach since it uses all the information available; however, the number of free parameters is huge (about a hundred).

Approach II lowers the number of free parameters by forcing the background and the 60 Co rates to be identical on detectors of the same group. In this approach each detector and data-set is still described individually by its own function f i, j (E) but the total number of free parameters is reduced since B i can assume only 3 values (for big, small and 130 Te-enriched crystals) and the same is true for Γ 60 Co i, j . This could be considered a strong assumption, but it is motivated by the fact that the low statistics prevent us from being sensitive to background variations among crystals of the same group in the 0νββ region. This method also offers the advantage of being less sensitive to fluctuations in the counting rate of a single detector over time, and takes into account the decay rate of 60 Co.

Approach III removes the background assumption of the previous model, at the price of a certain degree of information loss. The counting rate is simply averaged over all data-sets and detectors for the three mentioned groups. A variation of the background and of the 60 Co rate over time is then irrelevant, provided that the response function does not change with time. The average is done simply by summing over all the data collected with detectors belonging to the same group, thereby obtaining three spectra that can be represented by the function:

f k (E) = B k + Γ 60 Co k G k E -E Co + Γ 0ν G k E -E 0ν (6)
Here, the index k has three allowed values for big, small and 130 Te-enriched crystals while the response function G k (E) is de- fined as:

G k (E) = 1 i, j A i, j i, j A i, j √ 2πσ i, j exp        - (E -E 0 ) 2 2σ 2 i, j        (7)
where the sum over i extends on all the detectors belonging to the k th group and j runs over all the data-sets. A i, j and σ i, j are the corresponding background exposure and energy resolution measured during calibration. Note that -as is true also for the other two approaches -the response function is built using measured quantities, i.e. it does not contain any free parameters.

We discarded the first approach due to the excessive number of free parameters, and we performed two parallel (and independent) analyses on real and Monte Carlo-simulated data for the limit computation, following an unbinned likelihood technique [START_REF] Cowan | Statistical Data Analysis[END_REF] for the second approach and the standard CUORI-CINO Likelihood-Chi-Square technique [START_REF] Arnaboldi | [END_REF][START_REF] Arnaboldi | [END_REF] for the third one. The goal was to choose the most reliable procedure, checking for possible biases and comparing performances.

Comparison of methods on real data and Monte Carlo simulations

Applying the two methods on the CUORICINO RUN II data, we obtained results which were indeed very similar. Table 3 shows the best fits and limits for the two approaches: no significant difference is observed. In order evaluate the performances of the two approaches, we compared the distribution of the relevant statistical estimators for a thousand toy Monte Carlo simulated spectra for each CUORICINO detector, generated with rates corresponding to those directly measured in RUN II (see table 4) and no 0νββ signal. Several pieces of information have been obtained from this comparison:

• Both methods lead to compatible results. Figure 7 demonstrates a strong correlation between them.

• Both are unbiased. Figure 8 (left panel) shows that for both methods, the distribution of the best fits divided by their statistical errors is compatible with a gaussian centered at zero with a variance equal to one. This is an important result which is not always guaranteed by maximum likelihood methods applied to low statistics.

• Both have similar sensitivities. The distributions of the 90% confidence intervals (figure 7 top panel) show a sensitivity of nearly 2.5x10 24 y, evaluated as the median of the distribution, following the Feldmann and Cousins prescription [28].

The last point offers a nice synthesis for a better understanding of the numbers we are dealing with. A wide range of limits can be reached in experiments with the same true background level (figure 8, right panel) and therefore with the same sensitivity. In this respect, quoting only the limit reached by an experiment can be misleading if the sensitivity is not also mentioned.

Having such an impressive correspondence between the two approaches, the choice of the one or the other is somewhat arbitrary. We opted for the third method because it is consistent with the previous analysis of CUORICINO data [START_REF] Arnaboldi | [END_REF][START_REF] Arnaboldi | [END_REF] and because of its intrinsic simplicity.

Systematic uncertainties

In our analysis we have identified the following sources of systematic uncertainties:

• the calibration uncertainty;

• uncertainty in the signal efficiency;

• the background shape;

• the energy window of the fit.

In the case of the calibration uncertainty we have a direct estimate of its magnitude coming from a dedicated analysis of the residuals (see section 3): we reconstruct the position of a peak in the 0νββ region with a precision ∆E = ±0.8 keV. This systematic error has been included directly in the fit as a gaussian fluctuation in the energy position. Since this uncertainty is significantly larger than the error (0.013 keV) on the Q-value [START_REF] Redshaw | [END_REF], the systematic error discussed here also automatically includes the experimental uncertainty on the measured 130 Te transition energy.

To evaluate the uncertainty from the choice of energy window and shape of the background spectrum, we varied the model for the background (flat, linear and parabolic) at four different energy intervals centered at the Q-value of the 0νββ decay (we increased the lower and the upper bound of an energy window starting from 2527±30 keV at increasing steps of 5 keV). The results of this analysis have shown an average variation in Γ 0ν of 3 × 10 -26 y -1 .

The uncertainty on the signal efficiency is reported in table 2 to be 1.1% for the big crystals and 1.4% for the small crystals, both of which are negligible compared to the contributions from the energy scale and background parametrization uncertainties. . For each experimental result (rows 1-4) and for each of the considered NME evaluations (column 3-6), we report a m ee range. This identifies the upper bound on the neutrino Majorana mass according to the different results reported by the same author (when varying some of the parameters in the used nuclear model). In the two lower rows we compare the m ee range obtained for the 95% C.L. half-life limit on 130 Te (row 5) with the positive signal quoted by [14] (row 6). For this last case the m ee range corresponds to the 2 sigma range in the measured half-life.

The CUORICINO final result

As a next step, we added the contributions from big, small and enriched crystals from RUN I, combining their likelihoods with the RUN II data and using a similar reconstruction for the response function as described in section 6.1. The background rates are shown in table 4 while in figure 9 (left panel) we show the best fit and the corresponding 68% and 90% C.L. limits. Figure 9 (right panel) shows the logarithm of the combined likelihoods of RUN I and RUN II before and after the systematic uncertainties are included.

The resulting best fit for the 0νββ rate of 130 Te is:

Γ 0ν best = (-0.25 ± 1.44(stat) ± 0.3(syst)) × 10 -25 y -1
This result is compatible with zero, and the corresponding 90% C.L. half life lower bound is:

τ 0ν
1/2 ≥ 2.9 × 10 24 y Adding RUN I data, the limit improves more than might be naively expected from a 5% increase in exposure, since in general the limit scales with the square root of the exposure. However, as shown in figure 8, statistical fluctuations cause a spread in 90% C.L. limits that may be achieved by different experiments with the same exposure and sensitivity. In our case, in RUN I there was a downward fluctuation in the number of background counts in the region of interest, which led to a more stringent limit than could otherwise be obtained. However, this variation in the limit is well within the spread predicted by figure 8, which is why we prefer to quote the sensitivity of the experiment together with the limit. For the same reason, this final result is almost identical to the one published in [START_REF] Arnaboldi | [END_REF], although the total exposure has increased by a factor of ∼ 1.6.

It should be mentioned also that, in reference [START_REF] Arnaboldi | [END_REF], we used an older value of the 130 Te transition energy which had a much larger error and a slightly higher central value than the recently measured one [START_REF] Redshaw | [END_REF]31]. With the same data used in [START_REF] Arnaboldi | [END_REF], using the new result for the transition energy (with its smaller error) pushes the limit toward a lower half-life.

The inclusion of the systematic error modifies the likelihood profile for our data as shown in the right panel of figure 9. The profile can be considered the χ 2 of our fit as a function of all the possible Γ 0ν . Thus we will refer to it as χ 2 stat . If we adopt the hypothesis that our knowledge of Γ 0ν is smeared -near the best fit values -by a gaussian systematic uncertainty of magnitude σ syst , the total χ 2 tot will be:

1 χ 2 tot = 1 χ 2 stat + 1 χ 2 syst (8)
where the simplest approximated form of χ 2 syst is:

χ 2 syst = Γ 0ν -Γ 0ν best 2 σ 2 syst ( 9
)
With this modification of χ 2 and because the systematic uncertainty is small compared to the statistical error, we obtain a slightly weaker limit on the half life:

τ 0ν
1/2 ≥ 2.8 × 10 24 y As it is a standard approach in 0νββ literature, we also present the 95% C.L. limit on τ 0ν 1/2 including systematic uncertainties: τ 0ν 1/2 ≥ 2.3 × 10 24 y

Conclusion

In this paper we have presented the CUORICINO final result on 0νββ in 130 Te , obtained with an exposure of 19.75 kg•y of 130 Te, including a detailed study of systematic errors for the first time. A half life limit of 2.8×10 24 y at 90% C.L. is obtained (2.9×10 24 y if systematic errors are not included), to be compared (as discussed in section 6.2) with an experimental sensitivity4 of ∼2.6×10 24 y. This limit can be used to extract an upper limit on m ee using the theoretical NME evaluation for 130 Te nucleus. We report here results obtained using the most recent nuclear calculations found in literature:

• 300-570 meV using the Quasiparticle Random Phase Approximation (QRPA) evaluations of reference [5] • 360-580 meV using the QRPA evaluations of reference [6] • 570-710 meV using the Shell Model (SM) evaluations of reference [7] • 350-370 meV using the Interacting Boson Model (IBM) evaluations of reference [8] Note that, for each reference, a range (and not a single value) for m ee is presented, reflecting the different results for the NME obtained by the authors when varying model parameters, such as the treatment of the short range correlations or the value of g A (the axial-vector coupling). Then, the interval 300-710 meV can be taken as the final range for the 90% C.L. upper bound on m ee (at 95% C.L. this becomes 340-780 meV).

In table 5 we compare this result with the most stringent 90% C.L. half-life lower limits present in literature. For each experimental result we report the m ee range obtained with the NME evaluations here considered. Despite the differences between the NME evaluations, it is evident that CUORICINO is one of the most sensitive experiments performed to date.

Finally the bottom two rows in table 5 compare the 95% C.L. half-life limit on 130 Te obtained in this work with the 2 sigma range corresponding to the positive signal quoted by [14] and obtained with a re-analysis of the Heidelberg-Moscow data. The two results are clearly compatible.

  /(keV*kg*y)]

Figure 2 :

 2 Figure 2: A bolometric particle event (left) and a spurious non-physical signal (right), superimposed with their optimum filter output (in color) in the time domain.

Figure 3 :

 3 Figure 3: Residuals (nominal energy -calibrated energy) vs. nominal energy evaluated on the main gamma lines identified in the CUORICINO background spectrum. Circles (in color) refer to a calibrated energy obtained with the third order polynomial, triangles (in black) with the log-polynomial.

FWHMFigure 4 :

 4 Figure 4: Distribution of the energy resolutions (FWHM) measured in calibration for the three groups of crystals during the 33 data-sets belonging to RUN II. From top to bottom: big crystals, small natural crystals and 130 Te enriched crystals.

Figure 5 :

 5 Figure5: Typical scatter plot of the OFT (deviation of filtered raw signal from the average detector response) as a function of the signal's evaluated energy. The main trend identifies "good" events. Pile-up and non-physical pulses are outside the confidence regions (identified by the colored vertical bars). At low energy, there is a high density distribution of events for OFT values between 10 -2 and 10 -1 ; these events are due to electric disturbances.

Figure 7 :

 7 Figure 7: Results of a toy Monte Carlo simulation with no 0νββ signal (1000 simulated CUORICINO-RUN II experiments). Scatter plot of the 90% C.L. limits (top panel) and of the best fits (bottom panel) obtained with the two different approaches. The colored line has slope=1 and shows the strong correlation between the two techniques.

Figure 8 :

 8 Figure8: Results of a toy Monte Carlo simulation with no 0νββ signal (approach II in black and approach III in color). Left panel: pull distributions of the obtained best fits divided by their statistical error. Right panel: distribution of the 90% confidence level limit on the decay rate.

Figure 9 :

 9 Figure 9: Left panel: best fit, 68% and 90% confidence intervals for the total statistics (RUN I+RUN II) superimposed on the CUORICINO sum spectrum of the three groups of crystals (each scaled by efficiency and exposure) in the 0νββ region. (The purpose of the plot is to give a pictorial view of the result; in fact the fit was performed separately on 6 spectra whose likelihood are combined, as described in the text). Right panel: negative profile of the combined log likelihoods of RUN I and RUN II before (black) and after (color) the systematic uncertainty is included.

Table 1 :

 1 CUORICINO crystal information and statistics. Crystal mass is the average measured mass for CUORICINO detectors.

	Crystal Type	Crystal Mass	130 Te Mass	Exposure Run II	Exposure Run I
		[g]	[g]	[kg( 130 Te)•y]	[kg( 130 Te)•y]
	big	790	217	15.80	0.94
	small	330	91	2.02	0.094
	130 Te-enriched	330	199	0.75	0.145

Table 2 :

 2 Contributions to the CUORICINO 0νββ signal efficiency.

	Source	Signal efficiency (%)
	Energy escape	87.4±1.1 (big crystals) 84.2±1.4 (small crystals)
	Pulse-shape cuts	98.5±0.3
	Anti-coincidence cut	99.3±0.1
	Noise	99.1±0.1
	Pile-up with reference pulses	97.7
	Total	82.8±1.1 (big crystals) 79.7±1.4 (small crystals)
	energy of the signal (mainly due to noise and pile-up). Their
	contributions to the total signal efficiency are summarized in
	table 2.

Table 3 :

 3 Comparison between the fit results for the two studied approaches on RUN II data.

		Method II	Method III
	Best Fit [y -1 ]	(0.2 ± 1.5)×10 -25 (0.3 ± 1.5)×10 -25
	Half-life limit [y]	2.5 × 10 24	2.4 × 10 24

Table 4 :

 4 Background and 60 Co rates expressed in counts/keV/kg/y obtained by the combined fit of RUN I and RUN II data. The absolute rates are about 20% higher because the efficiencies have not been yet included. The high background rate obtained for 130 Te-enriched detectors is mainly due to a higher intrinsic contamination of these crystals.

	Fit Parameter	big crystals	small crystals	130 Te-enriched crystals
	RUN I: flat background rate	0.20±0.02	0.20±0.02	0.8±0.4
	RUN II: flat background rate	0.153±0.006	0.17±0.02	0.35±0.05
	RUN I: 60 Co rate	4.6 ± 1.5	9 ± 6	0 ± 14
	RUN II: 60 Co rate	2.5 ± 0.3	1.7 ± 0.8	0 ± 3.5

Table 5 :

 5 We compare the most stringent 90% C.L. half-life lower limits present in literature (column 2)

A second ββ candidate with high natural abundance is the isotope 128 Te (isotopic abundance 31.7%). This isotope is not as interesting as130 Te because of its lower transition energy, which reduces the nuclear factor of merit and shifts the signal to a higher background region (therefore also lowering the achievable experimental sensitivity).

At random times during the course of data acquisition, sets of 512 samples ("baselines") were collected in anticoincidence with the trigger. These were used for the evaluation of the average noise power spectrum, as described in the text.

To evaluate the energy resolution in the 0νββ region we use the 2615 keV peak since this is the nearest peak to the Q-value clearly visible in our calibration spectra.

This is the sensitivity evaluated for the total (RUN I + RUN II) statistics.

Acknowledgements

The CUORICINO Collaboration owes many thanks to the Directors and Staff of the Laboratori Nazionali del Gran Sasso over the years of the development, construction and operation of CUORICINO, and to the technical staffs of our Laboratories. In particular we would like to thank R. Gaigher, R. Mazza, P. Nuvolone, M. Perego, B. Romualdi, L. Tatananni and A. Rotilio for continuous and constructive help in various stages of this experiment. We are grateful to our colleagues of the CUORE collaboration for valuable suggestions. Among them we would like especially to name Yury Kolomenski for help and fruitful discussions. The experiment was supported by the Instituto Nazionale di Fisica Nucleare (INFN), the Commission of the European Community under Contract No. HPRN-CT-2002-00322, by the U.S. Department of Energy under Contract No. DE-AC03-76-SF00098, and DOE W-7405-Eng-48, and by the National Science Foundation Grant Nos. PHY-0139294 and PHY-0500337.