Hassane Aissaoui-Mehrez

Pascal Urien
email: pascal.urien@telecom-paristech.fr

Guy Pujolle
email: guy.pujolle@lip6.fr

Implementation Software To Secure Virtual Machines With Remote Grid of Secure Elements

Keywords: OpenID, Microcontrollers, Secure Elements, User-Centric Identity, Virtualization and Cloud Computing

Security for Future Networks (SecFuNet) is a Brazilian & European research project. The emerging Cloud of Secure Elements infrastructure is used for enforcing identity of Virtual Machines in the Cloud Computing. One of the main goals of the SecFuNet project is to develop a secure infrastructure for virtualized environments and Clouds that not only provides high availability and reliability for users, but that also provides strong isolation among virtual infrastructures. The project aims to develop a security framework for Cloud Computing and virtual environments. The goal of this paper is to describe the implementation and the experimentation of the solution for identifying users and nodes in the SecFuNet architecture. In this implementation, only authorized users are allowed to create or instantiate virtual environments. Thus, users and hypervisors are equipped with secure elements, used to open TLS secure channels with strong mutual authentication. Finally, since the physical substrates are shared by several resources (Users, VMs …), the proposed framework must ensure that one resource cannot interfere with the operations of another resource.

I. INTRODUCTION

While Cloud and virtual infrastructure services can offer great flexibility and convenience for its users, these users no longer have control over the platform on which their services are run. Not only users do not have any guarantees that their services have not leaked any sensitive information, but they may be also subject to attacks by other malicious users in the system. To address this issue, the SecFuNet project proposes to integrate the secure microcontrollers in order to introduce, among its many services, authentication and authorization functions for Cloud Computing and virtual environments. The Identity Management (IdM) system will be based on smart cards and user-centric attribute control policies. The authentication servers are composed by secure microcontrollers. The objective is to implement the framework, based on the authentication servers and EAP-TLS smart card. The proposed SecFuNet framework provides TLS secure channels for establishing trust relationships among Users, VMs, XEN and Remote Grid of Secure Elements (RGoSE). The authentication is done directly between smartcards (owned by users, XEN Hypervisor, or associated to VM) and a RGoSE arranged in a SecFuNet Identity Provider. This paper concerns a highly secure authentication server with an array of secure microcontrollers allowing VMs' strong mutual authentication with Xen hypervisor. It defines the structure and the components of the authentication server. The first part concerns the definition of the architecture of the server. The second part defines and implements the secure middleware within the Xen server to guarantee a secure communication between VMs and Xen and a response time to an EAP-TLS authentication of few seconds. This paper is organized as follows. Section II presents a concept of virtual environments and a brief state-of-art of related works on authentication with smart card and introduces the eap-tls smartcard concept. In Section III, we describe the software architecture of Xen Hypervisor and VMs layers. In section IV, we give an overview of the main features supported by our secure middleware in order to handle the remote smart cards plugged in Xen Hypervisor. Next, in section V, we present the experimental results and the performance of our architecture based on a grid of secure elements remotely accessed. Finally, we conclude this paper.

II. CONCEPTS AND PLATFORMS

The SecFuNet framework introduces state-of-the-art technologies to provide an experimental deployment and evaluation of virtual networks taking security as one of the major concerns. To address some of the security issues, the project aims to explore the application of secure elements, such as smart cards, to improve the trustworthiness of network infrastructure services for future networks. Two classes of secure microcontrollers have been studied, smart cards and TPMs (Trusted Platform Modules).

These electronics chips have different computing capabilities : a smartcard [START_REF] Jurgensen | Smartcards: The Developer's Toolkit[END_REF] is a tamper resistant device, including CPU, RAM and non volatile memory, most of these electronic chips support a Java Virtual machine (JVM) and execute software written in this programming language [START_REF] Chen | Java CardTM Technology for Smart cards: Architecture and Programmer's (The Java Series)[END_REF], and therefore are able to execute complex procedures (such as the TLS protocol), while TPMs are dedicated to the RSA algorithm. However these devices may be used in order to enforce trust for the TLS protocol or to guarantee secure storage for cryptographic keys. These security properties are directly provided by smart cards (thanks to dedicated embedded software), but require additional software components for TPMs.

In this document we call secure element a device such as smart card, which is able to totally or partially handle the TLS protocol and to realize the secure storage or computing of a cryptographic key. The main goal of this section is to briefly overview the context of the virtualization platforms and Grid of EAP-TLS smartcards protocol.

A. Architecture of a Grid of Secure Elements

In this paper we focus on arrays of EAP-TLS smartcards deployed in a Grid of Secure Elements (GoSE).

GoSE is a remote server (i.e. SecFunet Authentication Server, Xen Hypervisor…) hosting a set of secure elements (i.e. smartcards). All accesses to a GoSE require the TCP transport and are secured by the TLS protocol [START_REF] Dierks | The TLS Protocol Version 1.0[END_REF] [18] [START_REF] Dierks | The Transport Layer Security (TLS) Protocol Version 1.1[END_REF].

A GoSE offers services similar to HSM (Hardware Secure Module), but may be managed by a plurality of administrators, dealing with specific secure microcontrollers.

The goal of these platforms is to deliver trusted services over the Internet. These services are available in two functional planes:

-The user plane, which provides trusted computing and secure storage.

-The management plane, which manages the lifecycle (downloading, activation, deletion) of applications hosted by the Secure Element.

Every secure element needs a physical slot that provides power feeding and communication resources. This electrical interface is for example realized by a socket soldered on an electronic board, or a CAD (Card Acceptance Device), (i.e. a reader) supporting host buses such as USB "Fig. 1,". Within a GoSE each slot is identified by a SlotID (slot identifier) attribute, which may be a socket number or a CAD name. The SEID (Secure Element Identifier) is a unique identifier indicating that a given SE is hosted by a GoSE. It also implicitly refers the physical slot (SlotID) to which the SE is plugged.

The GoSE manages an internal table that establishes the relationship between SlotIDs and SEIDs. Therefore three parameters are needed for remote communication with secure element, the IP address of the GoSE, the associated TCP port, and the SEID.

B. EAP-TLS SmartCards Concept

EAP [START_REF]Extensible Authentication Protocol[END_REF] is a universal and flexible authentication framework. Because it can transport about any authentication protocol, it solves the interoperability concerns that their number and their disparity had risen.

Formally, EAP protocol is built on three kinds of messages: requests delivered by servers; responses returned by clients; and notifications issued by servers in order to indicate success or failure of authentication procedures.

The EAP smartcards functionality and binary encoding interface are detailed in [START_REF]EAP-Support in Smart card[END_REF]. These devices process EAP methods and act as server or client entity. They communicate via a serial link, whose throughput ranges between 9600 and 230,000 bauds.

There are two classes of operations, sending data (writing to smartcard) and receiving data (reading from smartcard); information is segmented in small blocks (up to 256 bytes) named APDUs, described by the ISO 7816 standard [20].

Security is enforced by multiple physical and logical countermeasures. The use of smartcards in TLS authentication has now a rather long history and has been largely developed according to different models. These devices run the OpenEapSmartcard JAVA open stack, introduced in [START_REF] Menon | Optimizing network virtualization in xen[END_REF] and which comprises four logical components "Fig. 1,": • The EAP Engine Object is mostly in charge of the I/O management (i.e. APDUs exchange). It is also responsible of EAP messages segmentation and reassembly. In fact, the APDU payloads maximum length is 255 bytes for input data and 256 bytes for output data, while EAP packets maximum length is about 1300 bytes of data. Consequently, EAP packets are split into several ISO 7816 units, and the Engine entity parses them in order to rebuild the proper EAP packet.

• The Credential Object holds all the credentials required by EAP-TLS method, that is to say: the Certification Authority certificate, the server Certificate and its associated private key. The EAP-TLS State Machine is reset and its according method is initialized with appropriate credentials, each time an EAP-Identity.Response message is received. This object also works as an Identity module. For now, it only holds a unique server Identity but one could possibly load different server Identities issued by several companies which, upon success, would grant different kind of access or services depending of the Client's Identity and its subscription to one or several companies' Network.

Fig. 1. EAP-TLS Smartcard components and GoSE

• The Authentication Interface object implements all services fulfilled by EAP-TLS methods, whose main procedures are initialization, packet processing or MSK key downloading.

• Lastly, the EAP-TLS object is in charge of packets processing, as specified in EAP standard [START_REF] Urien | Remote APDU Call Secure[END_REF]. Since TLS packets size may exceed the Ethernet frame capacity, EAP-TLS supports internal segmentation and reassembly mechanisms.

C. Virtualization

More recently, as technologies and applications matured the demand for novel network functionality and requirements began to rise. However, alterations to the Internet architecture have become restricted to simple incremental updates and deployment of new network technologies have become increasingly difficult [START_REF] Chowdhury | A survey of network virtualization[END_REF].

This issue became commonly known as "Internet ossification" [START_REF] Turner | Diversifying the internet[END_REF] or "Internet tussle" [START_REF] Clark | Tussle in cyberspace: defining tomorrow's internet[END_REF]. In response, network virtualization has been proposed as a means to overcome this impasse by providing diversification to the future internetworking paradigm [START_REF] Anderson | Overcoming the internet impasse through virtualization[END_REF].

Virtualization is pushing computer systems to new frontiers. From virtualized data centers to cloud providersscalability, elasticity and multi-tenancy are some of the key required or desired features. This technology has been mostly applied to manage computing and storage infrastructures.

Virtualization techniques create abstract virtual environments composed of a set of "sliced" resources, such as disk, memory, CPU, flow or network. Each virtual environment is isolated and independent from others and thus it seems to run directly over the shared underlying substrate. The virtualization abstraction is accomplished by a software layer called hypervisor, which provides virtual environment interface abstractions quite similar to underlying ones [START_REF] Fernandes | Virtual networks: Isolation, performance, and trends[END_REF].

D. Xen Hypervisor

The Xen software is a layer called hypervisor, which provides the virtual abstraction, quite similar to underlying physical node. Xen is an open-source that implements a virtualization layer over the physical hardware proposed to run on commodity equipment [START_REF] Urien | Cloud of Secure Elements, Perspectives for Mobile and Cloud Applications Security[END_REF].

The hypervisor creates a virtual environment referred to as Virtual Machine (VM). Each virtual machine accesses a resource abstraction that has its own CPU, memory, disk, and network interface, with its own operating system and applications. The hypervisor controls the access to the hardware and also manages the available resources shared by virtual machines.

A VM is isolated from others, i.e. a running virtual machine is expected not to affect the behavior of other existing virtual machines. Since Xen virtualizes hardware, a VM is an entire computer system, with its own hardware abstraction resources. Xen architecture is composed of one hypervisor located above the physical hardware and several virtual machines over the hypervisor, as shown in "Fig. 2,". Xen virtualizes the processor by assigning virtual CPUs (vCPUs) to virtual machines and implements a CPU scheduler that dynamically maps a physical CPU to each vCPU during a certain period, based on a scheduling algorithm [START_REF] Couto | Xtc: A throughput control mechanism for xen-based virtualized software routers[END_REF]. Through the paravirtualization technique, Xen enhances guest system performance by changing its behavior to call the hypervisor when necessary. Memory virtualization in Xen is static and the physical RAM is divided between virtual machines. Each virtual machine accesses a fixed amount of memory space, specified at the time of its creation.

In addition, device drivers are kept in a special virtual machine, called the driver domain. The driver domain tasks are often executed by Domain 0, a privileged virtual machine that also executes hypervisor managing tasks.

The driver domain access I/O devices directly by using its native device drivers and also performs I/O operations on behalf of virtual machines.

On its turn, a virtual machine, also called unprivileged domain or Domain U, employs virtual I/O devices controlled by virtual drivers to request Domain 0 for device access [START_REF] Menon | Optimizing network virtualization in xen[END_REF].

III. SOFTWARE ARCHITECTURE OF HYPERVISOR AND VMS

In this section, we detail components composing the infrastructure and platform layer, as illustrated by "Fig. 3,".

The SecFuNet architecture distinguishes two classes of entities: Nodes and Users. The SecFuNet Users are equipped with two classes of secure elements: EAP-TLS smartcards or PKCS#1 smartcards, used to open TLS secure channels with strong mutual authentication.

The SecFuNet nodes are divided into three layers of visibility, whose components are identified by a certificate.

The first one (the infrastructure layer) directly uses the services of a secure element, the platform layer and the application layer. Each of them is associated with an administrator in order to manage it. The infrastructure layer includes hardware resources (servers, secure element grids...) and hypervisors (HV) running in Domain0. Each hypervisor has a certificate, which instantiates its identity. It is associated to a secure microcontroller storing its private key.

TLS sessions with strong mutual authentication are mandatory for all virtual machines (VM). Thus, we implement the Remote TLS Application in VM-side and the Secure Middleware in Xen-side, in order to enforce a remote VM access to grid of secure elements (GoSE).

A. Secure Middleware API

The Secure Middleware "Fig. 4," is set of software, includes a very simple facility to access GoSE or smart card readers plugged on Xen Hypervisor, with security wrapper mechanisms. This middleware offers some features to secure connections and transactions between GoSE, Xen and VMs. However, the VM may access securely to GoSE and handles the associated remote secure element. The Secure Middleware API is split into several layers:

• PC/SC Driver: talks to GoSE hardware (i.e. readers of Smart Cards). We included the open source project (OpenSC-Tools) [14] for smart card. It implements support for many drivers at the same time. OpenSC is written by an international team and is licensed as Open Source software under the new BSD license. • EAP-TLS API: is a set of tools and libraries dealing and working with cryptographic capabilities of smartcards. It deals with the credential component that holds all the credentials required by EAP-TLS method (Certification Authority, Server Certificate and its associated private key). The main features of this API are: reset and start the EAP-TLS State Machine of smart card and initialize it with the appropriate identity and credentials.

• RACS: Remote APDU Call Secure protocol for Hardware Secure Module (HSM) Facilities, introduced in [START_REF] Urien | Remote APDU Call Secure[END_REF][15] [START_REF]PKCS#11[END_REF]. It provides all the features needed for the remote use of secure elements. The RACS protocol works over the TCP transport layer and is secured by the TLS protocol. One of the main targets of the RACS protocol is to efficiently push a set of ISO7816 requests [START_REF]EAP-Support in Smart card[END_REF] towards a secure element in order to perform cryptographic operations. In that case a RACS request typically comprises a prefix made with multiple ISO7816 requests and a suffix that collects the result of a cryptographic procedure.

• TLS server: secures connections and transactions between Xen server and VMs. EAP-TLS API launches a TLS server in Xen-side. A TLS Client is running on VM thanks to Remote TLS Client Application. A TLS session is opened between Xen and VM. Then, the Remote TLS Application may handle the remote Smart Card plugged on Xen server using the functionalities provided by EAP-TLS API and RACS. TLS sessions are opened. Upon success, a set of ephemeral cryptographic keys (named KeyBlock) and the list of negotiated algorithms (the CipherSuite parameter) are read from the smart card device.

B. Remote TLS Client Application

Remote TLS Client Application allows a secure connection between VM and Xen and transports ISO7816 commands. When the secure channel is established, the VM can send the APDU commands to smart card via a Secure Middleware API. Sessions with Xen are protected by the TLS protocol with mutual authentication based on certificates. VMs handle the smart card, thanks to RACS Client [START_REF] Urien | Remote APDU Call Secure[END_REF] and EAP-TLS API which hide the encoding of remote procedure call over the networks.

IV. THE FUNCTIONS SUPPORTED BY EAP-TLS API

The EAP-TLS API consists of number of functions: slot management and cryptographic instructions. This section indicates the main functions that are implemented in the Secure Middleware Library. The main procedures are the following:

• SCard_init_remote_reader: initialize the smart card reader, this function return the number of detected smart card readers, selects the javacard application identity (AID) in smart card device, sets identity and verify code pin, starts a session with a smart card identified by its index (first index is set to zero), setup the T-Protocol used by a smart card.

• SCard_start_remote_eap: resets the current EAP-TLS state machine and processes an EAP message.

• SCard_connect_remote_reader: opens a connection whit smart card reader.

• SCard_disconnect_remote_reader: closes a session with a smart card identified by its index.

• SCard_list_remote_reader: list all remote smart card readers.

• SCard_send_remote_apdu_cmd: send ISO7816 commands to a smart card.

• SCard_get_remote_cipher_suite: collects the negotiated CipherSuite parameter from smart card device, associated to an opened TLS session.

• SCard_get_remote_key_block: collects the ephemeral keys from smart card device, associated to an opened TLS session.

• SCard_parse_rw_message: read responses from smart card parses it and send messages to VM.

• SCard_process_eap_tls_packet: receives instructions from VMs and make an eap-tls packet to send to smart card.

The Remote TLS Client cans remotely inventory and exchange information with the secure elements. The RACS protocol and EAP-TLS API manage the GoSE or Smart Card and products ease authentication, encryption and digital signatures. This Secure Middleware allows to perform cryptographic operations and access to private key in the smart card associated to VM, and that access is granted only after the VM has been authenticated.

For example: The VM is running in the Domain U of hypervisor. Each VM holds a certificate (that instantiates its identity), but its private key is stored in a secure element plugged in a SIM array (GoSE). All cryptographic procedures dealing with secret keys (such as signing, symmetric encryption or decryption) are performed by the secure element associated to this VM. The protection of asymmetric cryptographic key is a critical issue for the SecFuNet Node.

The authentication between the VM and a Xen server requires the HV private key and VM-Auth-Token (VM Authentication Token). The VM-Auth-Token establishes the proof that the VM (identified by its certificate ID VM ) is attached to the HV (identified by its certificate ID HV ).

Conceptually the VM-Auth-Token is a signed piece of information (such as certificate) comprising ID VM , ID HV and a validity date; it works with HV public key. Therefore TLS sessions dealing with this certificate imply the use of the HV private key.

VM-Auth-Token = {ID VM ,ID HV ,Validity-Date}Sign-Priv HV .

V. THE EXPERIMENTAL RESULTS

We selected for the SecFuNet project the secure microcontrollers with timings are expressed in ms. Basic cryptographic performances are illustrated by "Fig. 5,": • Digest costs (T MD5 =0.49ms and T SHA1 =0.93ms). It is the time needed by hash functions to process a block of 64 bytes.

• Symmetric cryptographic costs (T 3xDES =2.10ms, T AES =2.60ms,T RC4 =0.50ms). It is the time required by symmetric encryption/decryption operations. 3xDES work with 8-byte blocks, while AES deals with 16-byte blocks. It should be noticed that RC4 is a stream cipher algorithm.

• Asymmetric cryptographic costs (T RSAPub =25ms and T RSAPriv =560ms). It is the time required for encryption and decryption with public and private RSA keys (whose size is 128-byte).

• T RW is the time measured by the docking host, to write and read data to/from the non-volatile memory (T RW # 0,17 ms/byte) Tests were carried out to assess the performance of the developed prototype and attest the infrastructure is viable.

The test configuration is a single server (Xeon Intel Dual Core 3.2 GHz, 12GB RAM), running a Xen hypervisor version 4.

Tow VMs (Fedora 20) were instanced on top of the hypervisor. We do not have access to the actual GoSE, so, the Authentication Server is executed in domain0 with two plugged smart cards in server.

The connection of VM to Xen server is done in three stages:

• Phase 1 is the time to establish secured connection with the Xen server, using tls session between VM and Xen server.

• Phase 2 is the Resume or Full TLS session: the time to establish a connection with the authentication server, using the smartcard associated to the VM. The measurement (T2) is the time between the moment the administrator confirms the PIN and the receipt of the finished message from AS, inside the TLS channel.

• Phase 3 is the time once the secure TLS channel has been established.

Our measures focus on the times of phase 2 and 3 that use transactions with the smart card:

• Phase 2 -Full TLS session: the opening of a full TLS session based on 1024 bits RSA cryptography, requests 2,0s. It exchanges 2500 bytes, whose transfer costs 0,40s. About 230 MD5 and SHA1 calculations are performed (dealing with 64-byte blocks), which consumes 0,32s. One RSA operation with the private key and two RSA procedures with public keys are computed in 0,60s. The sum of the previously cited operations is 1,32s; the remaining time (0,68s = 2,00 -1,32) is burnt by the java code execution.

• Phase 2 -Resume TLS session: the opening of a resume TLS session costs 0,50s. It requires the exchange of 250 bytes (0,04s), and the calculation of 75 MD5 and SHA1 that consumes 0,11 s. The remaining time (0,35s) is needed by the java code execution.

• Phase III : Once the secure TLS channel has been established, messages exchanged by the client and the server are encrypted with the selected algorithm and integrity is enforced by a HMAC procedure. This phase may either be computed by the EAP-TLS smartcard, or either exported in the docking terminal.

In the case where phase III is processed by the EAP-TLS device, we deduce according to "Fig. 5," that the computing time is mostly consumed by encryption and decryption operations (in other word hash calculations have a second order effect). As an illustration, for a classical RC4, HMAC-MD5 cipher suite, the transfer of 1000 bytes burnt 500ms for encryption, 16ms for MD5 hash, and 170ms for data transmission. Therefore it is mandatory to transfer TLS session from the secure microcontroller to the computer to which it is plugged, in order to obtain an operational data throughput.

VI. CONCLUSION

The architecture proposed and implemented in this deliverable presents a secure solution for VM connection, which aims to develop a highly secure identification scheme, for Cloud Computing, based on the secure elements (smart cards and grid of secure elements). This integration provides a higher level of security when compared to traditional password authentication, and establishes trust relationships among the resources, and eliminates the vulnerabilities like phishing attacks. Strong authentication is done directly between smart cards (owned by user or associated to VM) and a grid of secure processors arranged in a SecFuNet.

We have seen in this paper that the use of secure elements can make a real highlight in the certificate management and the application of a security policy. Adding our Secure Middleware, for large infrastructure, brings the real simplicity of management. In parallel, can be very interesting solution to access securely to a remote specific area of GoSE.

The experimental results of the platform developed for SecFuNet demonstrates that the scalability performances are compatible with today constraints. We are confident, that we will be able to achieve a platform whose authentication time will be reasonable enough to be massively deployed. Furthermore, the smart card shall be a great addition to the Xen architecture and will be as well as a key asset to securing Cloud Computing infrastructures.

Fig. 2 .

 2 Fig. 2. Xen Hypervisor Architecture

Fig. 3 .

 3 Fig. 3. Xen Hypervisor And Software Architecture

Fig. 4 .

 4 Fig. 4. Remote TLS Client and Secure Middleware

Fig. 5 .

 5 Fig. 5. Digest costs for the different elementary cryptographic operations.

Acknowledgment

This work has had financial support from CNPQ through process 590047/2011-6 (SecFuNet project) and also through processes 307588/2010-6 and 384858/2012-0. We also thank CAPES for the financial support with PhD scholarship.


Ting Wang
email: ting.wang@inria.fr

Claire Dune
email: claire.dune@univ-tln.fr

Jean-Pierre Merlet
email: jean-pierre.merlet@inria.fr

Philippe Gorce
email: philippe.gorce@univ-tln.fr

Guillaume Sacco

Philippe Robert

Jean-Michel Turpin

Bernard Teboul
email: bernard.teboul@libertysurf.fr

Audrey Marteu

Olivier Guerin
email: guerin.o@chu-nice.fr

A new application of smart walker for quantitative analysis of human walking

Keywords: Smart walker, gait analysis, step detection, turning, elderly

This paper presents a new nonintrusive device for everyday gait analysis and health monitoring. The system is a standard rollator equipped with encoders and inertial sensors. The assisted walking of 25 healthy elderly and 23 young adults are compared to develop walking quality index. The subjects were asked to walk on a straight trajectory and an L-shaped trajectory respectively. The walking trajectory, which is missing in other gait analysis methods, is calculated based on the encoder data. The obtained trajectory and steps are compared with the results of a motion capture system. The gait analysis results show that new index obtained by using the walker measurements, and not available otherwise, are very discriminating, e.g., the elderly have larger lateral motion and maneuver area, smaller angular velocity during turning, their walking accuracy is lower and turning ability is weaker although they have almost the same walking velocity as the young people.

Introduction

Ageing in society is a worldwide issue that especially impacts developed countries. In France, due to the high care cost and the limited number of rooms in care institution, the solution that has been chosen by care givers, frail people and their family is to maintain elderly at home and in the best conditions by giving them an adapted assistance. Among the possible assistance devices, the walkers have large number of users because of their simplicity while using the person's remaining locomotion capability.

Advances in robotics make it possible to develop more intelligent walkers by adding sensors and actuators. According to the user's needs, the functions of the proposed walkers are not restricted to their primary tasks, i.e. physical support and mobility assistance [START_REF] Alwan | Stability margin monitoring in steering-controlled intelligent walkers for the elderly[END_REF], [START_REF] Bostelman | Robotic patient transfer and rehabilitation device for patient care facilities or the home[END_REF], [START_REF] Glover | A robotically-augmented walker for older adults[END_REF]. There are other functions such as sensorial assistance, cognitive assistance and health monitoring [START_REF] Frizera | The smart walkers as geriatric assistive device. the simbiosis purpose[END_REF], [START_REF] Spenko | Robotic personal aids for mobility and monitoring for the elderly[END_REF]. Some walkers also focus on sit-to-stand transfer [START_REF] Chugo | A moving control of a robotic walker for standing, walking and seating assistance[END_REF], [START_REF] Médéric | Elderly people sit to stand transfer experimental analysis[END_REF], navigation help [START_REF] Grasse | Assisted navigation for persons with reduced mobility: path recognition through particle filtering (condensation algorithm)[END_REF], [START_REF] Kulyukin | iWalker: toward a rollator mounted wayfinding system for the elderly[END_REF], obstacle avoidance and fall detection [START_REF] Hirata | Fall prevention of passive intelligent walker based on human model[END_REF].

To study the extension of the functions of walkers we have developed our own family of walking aids based on the walker presented in [START_REF] Merlet | ANG, a family of multi-mode, low cost walking aid[END_REF]. Our walkers can perform multi-functions such as navigation, street mapping, fall detection/prevention and autonomous object recovery. In this work we use the simplest version (Fig. 1), which is based on a commercially available 3-wheels Rollator. It was instrumented with encoders at the two fixed rear wheels and a 3D accelerometer/gyrometer at the front with the purpose of determining the walker's trajectory on a 24/24 basis. A small, low energy consumption fit-pc computer manages the measurements and records all the data. Compared with the walkers proposed above, our walker is low cost, simple to be used at home and its functions can be easily extended. This paper will present how it can be used for medical monitoring of walking patterns.

Many studies have examined the effect of age on the walking by comparing younger with older adults [START_REF] Hollman | Age-related differences in spatiotemporal markers of gait stability during dual task walking[END_REF], [START_REF] Menz | Age-related differences in walking stability[END_REF], [START_REF] Woledge | The variable component of lateral body sway during walking in young and older humans[END_REF], [START_REF] Owings | Variability of step kinematics in young and older adults[END_REF], [START_REF] Grabiner | Age-related changes in spatial and temporal gait variables[END_REF], [START_REF] Hausdorff | Increased gait unsteadiness in community-dwelling elderly fallers[END_REF], [START_REF] Barbara | Gait characteristics of younger-old and older-old adults walking overground and on a compliant surface[END_REF], [START_REF] Su | Age-associated differences in the gait pattern changes of older adults during fast-speed and fatigue conditions: results from the baltimore longitudinal study of ageing[END_REF]. These studies calculated some gait parameters, such as step length, gait cycle, step width, cadence and gait speed [START_REF] Su | Age-associated differences in the gait pattern changes of older adults during fast-speed and fatigue conditions: results from the baltimore longitudinal study of ageing[END_REF], [START_REF] Callisaya | Sex modifies the relationship between age and gait: a population-based study of older adults[END_REF], [START_REF] Zijlstra | Assessment of spatio-temporal gait parameters from trunk accelerations during human walking[END_REF], [START_REF] Studenski | Gait speed and survival in older adults[END_REF]. Some studies have considered also gait variability [START_REF] Owings | Variability of step kinematics in young and older adults[END_REF], [START_REF] Callisaya | Ageing and gait variability-a population-based study of older people[END_REF], [START_REF] Grabiner | Age-related changes in spatial and temporal gait variables[END_REF], [START_REF] Hausdorff | Gait variability and fall risk in communityliving older adults: a 1-year prospective study[END_REF]. However, these usual gait parameters are not sufficient and sensitive enough to evaluate the health state of elderly people, and the results of previous studies are often inconsistent according to the conditions of experiment. Another drawback of most studies is that these measures are obtained on a reduced space with specialized laboratory equipment such as motion capture systems and instrumented walkways, which may not be available in many clinics and certainly not during daily activities. Thus the recent studies used the wearable inertial sensors to analyze the fall risk of elderly people [START_REF] Higashi | Quantitative evaluation of movement using the timed up-and-go test[END_REF], [START_REF] Zakaria | Quantitative analysis of fall risk using tug test[END_REF], [START_REF] Salarian | Analyzing 180 degrees turns using an inertial system reveals early signs of progression of parkinson's disease[END_REF]. Thanks to encoders and a 3D accelerometer/gyrometer, by using the smart walker we can not only calculate accurately all the classical gait parameters proposed above, but also avoid imposing constraints on the end-users. We also can obtain the trajectory of the walker and therefore compare it with the reference trajectory. This comparison will allow to establish several original gait parameters which have not been considered in previous studies. Overall, using our walker the gait characteristics can be described more comprehensively.

Methods

This section will explain how the walker can obtain the walking trajectory and determine the stride while the subject is walking along a straight line. In addition, in order to validate our method, the obtained gait features will be compared with the data given by a 6-camera motion capture system which is used for biomechanical analysis [START_REF] Richards | The measurement of human motion: A comparison of commercially available systems[END_REF].

Calculation of the trajectory

As shown in Fig. 2, the origin of the walker frame O 1 X 1 Y 1 is supposed to be the position of the middle point between the two rear wheels. The position of the walker in a reference frame OXY is described by [x, y, θ], where θ describes the walking direction of the rollator and represents the angle between the horizontal axis of two rear wheels and the X axis. In our experiment of 10m straight line walk test, as the reference trajectory was directed along the Y axis we have θ = 0 at the beginning of the walker's trajectory. The trajectory of the walker is determined by using the encoders. Assuming that at the (k + 1) th time sample moment the measurement of the encoders of two rear wheels are ∆ L and ∆ R , the displacement of the left and right wheel are obtained by using (1):

dL = 2πr 4C • 360 ∆ L , dR = 2πr 4C • 360 ∆ R ( 1
)
where r is the radius of the rear wheel and C is a constant parameter of the transformation between the value of encoder and the wheel radius. The change of the direction angle θ during the (k + 1) th sampling time can be estimated as:

dθ = (dL -dR)/D, (2)
where D is the distance between the two rear wheels.

According to the kinematic model shown in Fig. 2, the changes of the walker's position can be obtained as follows [START_REF] Lee | Accurate calibration of kinematic parameters for two wheel differential mobile robots[END_REF]:

dx = dL + dR 2 sin(θ k + dθ 2 ), dy = dL + dR 2 cos(θ k + dθ 2 ) (3)
Finally, the new position of the walker can be calculated by using:

x k+1 = x k + dx, y k+1 = y k + dy, θ k+1 = θ k + dθ. (4)
Using the above equations, the trajectory of the walker can be determined by using the encoders. The experiments have shown that after a straight line walking of 10 meters the estimated positioning has an absolute accuracy better than 1cm.

Step detection

Accelerometers are utilized to detect the walking stride in many studies [START_REF] Zijlstra | Assessment of spatio-temporal parameters during unconstrained walking[END_REF], [START_REF] Huang | An orientation free adaptive step detection algorithm using a smart phone in physical activity monitoring[END_REF]. Most of methods use the peak value of forward acceleration to detect the walking cycle. However, some steps do not lead to a high-peak forward acceleration, and hence they are not counted although there is displacement during these periods. A recent study [START_REF] Rebula | Measurement of foot placement and its variability with inertial sensors[END_REF] used thresholds on the magnitude of the gyroscope and accelerometer signals to identify the zero velocity instant and regarded it as the end of a step. Our walker also uses the gyrometer data to detect the walking stride. An interesting contribution is that it allows one to differentiate the right and left steps when the subject walks along a straight line. Indeed when the subject is on the left (right) support phase the walking aid rotates on the left (right). Hence the rotational velocity of the walker around the vertical axis, which can be easily obtained by the gyrometer, is used to detect the walking stride. Its zero value instant is regarded as the end of a step. An example of measured rotational velocity after a low pass filtering for an elderly people during walking along a 10m straight line is shown in Fig. 3. Since the position of the walker at every moment has been calculated by using the method presented in Section 2.1, the displacement of the walker during every step, which is regarded as the step length of the subject, can be easily calculated as soon as all the steps are detected, as shown in Fig. 4. Accordingly, all gait speed characteristics (such as mean value, minimum and maximum value) can be obtained for each step. With a sampling time of 1ms for the encoders and 4.8ms for the gyrometer, we may obtain a quiet reasonable accuracy on these parameters.

Validation

The proposed methods are tested on a standard 4-wheels rollator equipped with sensors and monitored by a motion capture system, as shown in Fig. 5. We attached 10 optical markers on the walker frame, wheel axis and handle. The motion capture gives the 3D trajectory of the markers at a framerate of 200Hz. The motion capture system is a high cost, whole body motion analysis system widely used for bio-mechanical analysis. Its accuracy is studied in [START_REF] Richards | The measurement of human motion: A comparison of commercially available systems[END_REF] among other commercially available systems for the measurement of human motion. The constructors give an accuracy of 0.6cm in a 3-m field of view. In [START_REF] Richards | The measurement of human motion: A comparison of commercially available systems[END_REF], the accuracy is evaluated by measuring the relative distance between two markers set on a rigid target of known size (9cm) that rotates around a vertical axis. The markers are seen by at most 3 cameras at each time frame. The results show a maximum distance error of 1.2cm. Yet, no result is given on the absolute distance error.

Validation of trajectory calculation Fig. 6 a. and Fig. 6 c. compare two walker trajectories (straight line and U-turn) that are estimated using the encoder values with the trajectories obtained using the motion capture system. The instantaneous velocity can also be estimated by mapping the time frames to these measurements (Fig. 6 b. and Fig. 6 d.).

The straight line experiment shows that the distance estimation given by the odometers and the motion capture system are similar. The encoders are able to capture small deviations from the straight line. The U-turn experiment also shows that the trajectory of the walker is well described, even when performing maneuver. Then, the encoders' values can be used to describe the walker trajectory and they give similar results as the motion capture system.

Validation of step detection The motion capture system tracks a marker set attached to the feet of the person. It gives the Cartesian position of the feet in a fixed Cartesian reference frame, which can be used to detect the steps. The velocity of the feet can also be used to determine the number of steps and their length. For example, the step starts when the velocity increases while it ends when the velocity is below an arbitrary threshold, as shown in Fig. 7.

The number of steps obtained using this method was compared with the result of the method proposed in Section. 2.2, as shown in Fig. 7. The first graph depicts the step detection by finding the zero value of the rotational velocity. The second graph represents the step detection by finding the sagital velocity changes. The third graph, at the bottom, compares the results of the two methods. 6. Walker's trajectory and velocity estimation using the encoder (odometry) and the motion capture system for two experiments: walking on a straight line (a. trajectory and b. velocity) and U-turn (c. trajectory and d. velocity).

It is worth mentioning that the number of steps shown in the first graph depends on the parameter of the low pass filter used for the angular velocity. If the parameter is well chosen, the number of steps is quite the same as measured with the motion capture system (see Fig. 7). If the low pass filter is not correctly tuned, the starting moment of the steps is not accurately found and that will impact the calculated step length and cadence.

In summary, our smart rollator gives results that are similar to motion capture data while being cheaper, more convenient to use and faster to set up. Indeed, the motion capture has to be calibrated for each experiment, the marker set up takes time and is sensible to skin motion, it is also dependent from the abilities of the system operator and it is restricted to a small area (3-m field of view for whole body tracking with a 6-camera system) with good illumination condition. On the opposite, the walker only needs to be turned on, no calibration is required, the measurements are totally independent from the system operator, the user does not need to wear any marker and it is insensitive to condition while it can be used for large displacements. Furthermore, the motion capture has a framerate of 200Hz, while the encoders and inertial sensors are 1000Hz. st ep Left > 0 Right < 0 Fig. 7. Comparison of the step numbers estimated with the gyroscope vs ground truth given by the motion capture system. The rotational velocity values are filtered with a Butterworth low pass filter with a coefficient of 0.02.

Gait analysis of a 10m walk test

Several classical tests used to assess the mobility of elderly people are the 10m walk test (10mWT) (measure: time duration) [START_REF] Salbach | Responsiveness and predictability of gait speed and other disability measures in acute stroke[END_REF], Timed Up and Go test (TUG) (measure: time duration) [START_REF] Podsiadlo | The timed 'up & go': a test of basic functional mobility for frail elderly persons[END_REF], Tinetti Test (TT) (analysis of gait parameters through a video) [START_REF] Tinetti | Performance-oriented assessment of mobility problems in elderly patients[END_REF]. Such tests are easy to implement but are basically global (the time for the 10mWT and the TUG may be identical for two subjects which have however very different walking patterns) or are subjective (for the TT [START_REF] Panella | Towards objective quantification of the tinetti test[END_REF]). Furthermore these tests are performed only during medical visits and consequently are not appropriate to detect abnormal events in the walking patterns.

Our objective was to examine if the measurements of our walking aid allow one to refine the output of the above walking tests.

For that purpose we have led a large scale experiment that was approved by the regional ethical committee (Comité de Protection des Personnes). Each subject was asked to walk along two specific trajectories while using the walker. One is a 10m straight line and another is an inverted L-shaped trajectory (see Fig. 12). The subjects are a group of young adults of INRIA (with age between 25 and 65 years, mean value 32) and a group of elderly people (age over 65 years) of Nice hospital. No subject has pathological walking diseases and all of them were asked to perform twice the trajectory with the walking aid.

Results

There were 23 young adults and 25 elderly people that participated to the 10mWT. The trajectories of all the subjects for the 10mWT are presented in Fig. 8. The reference trajectory is the horizontal axis and the vertical axis is scaled to illustrate the lateral deviations between the real and reference trajectories. Fig. 8 clearly shows that the elderly subjects have larger deviations than the young subjects. The maximum lateral deviation between the real and reference trajectory of the two groups is presented in Fig. 9. It illustrates that the results of the elderly have a significantly larger deviation than that of the young subjects. For example, the mean value of the maximum lateral deviation for the elderly is 11.048 ± 5.99 cm while that for the young people is only 3.963 ± 2.301 cm. In addition, several other walking indicators calculated from a trajectory, e.g., the area between the real and reference trajectory, the traveled Manhattan distance and the orientation angle of the walker also reveal that the lateral motion of the elderly is larger than the younger. It has been has shown that there exists significant group difference in the medio-lateral displacement of the center of mass between healthy elderly adults and elderly patients [START_REF] Hahn | Can motion of individual body segments identify dynamic instability in the elderly[END_REF] . It is consistent with our result and they reveal that walking accuracy can be regarded as a pertinent walking quality index.

Since the user's step can be detected by using our walker, most of gait parameters presented in the usual walking tests can be calculated or estimated, such as step length, step period, gait speed, cadence, forward acceleration and variability of these gait parameters. Although the step width cannot be calculated accurately, the analysis of the walker's lateral motion in the previous section can reflect the characteristic of the subjects' step width.

Usually gait speed or walking velocity is regarded as a very important indicator of health. Some studies claimed that older subjects exhibited significantly reduced gait speed compared to younger adults [START_REF] Hollman | Age-related differences in spatiotemporal markers of gait stability during dual task walking[END_REF], [START_REF] Menz | Age-related differences in walking stability[END_REF], [START_REF] Su | Age-associated differences in the gait pattern changes of older adults during fast-speed and fatigue conditions: results from the baltimore longitudinal study of ageing[END_REF] while other studies showed that there were little or no differences in that [START_REF] Grabiner | Age-related changes in spatial and temporal gait variables[END_REF], [START_REF] Hausdorff | Increased gait unsteadiness in community-dwelling elderly fallers[END_REF].

It has been shown that the gait speed depends on age, sex [START_REF] Callisaya | Sex modifies the relationship between age and gait: a population-based study of older adults[END_REF], use of mobility aids, chronic conditions, smoking history, blood pressure, body mass index, and hospitalization [START_REF] Studenski | Gait speed and survival in older adults[END_REF]. The instantaneous walking velocity can be derived from the encoder measurements and Fig. 10 gives their maximum values. It shows that the maximum walking velocity of the younger subjects is a little larger than that of the elderly. For the elderly the mean value of the maximum walking velocity is 117.969 ± 15.851 cm/s and for the young people this value is 119.967 ± 16.019 cm/s. Hence there is no obvious difference between the two groups, which is consistent with the result of [START_REF] Grabiner | Age-related changes in spatial and temporal gait variables[END_REF], [START_REF] Hausdorff | Increased gait unsteadiness in community-dwelling elderly fallers[END_REF].

We have also compared step period, step length and their variability for the two groups. Again, there is almost no difference between the two groups, which explains why the two groups have similar walking speed. The mean values of step period and step length of the elderly are 0.526 ± 0.1 s and 54.862 ± 11.643 cm, and 0.537 ± 0.095 s and 55.050 ± 8.605 cm for the young people.

Gait analysis of an inverted L-shaped trajectory test

An inverted L-shaped reference trajectory is described in Fig. 12. During the test, the subject was asked to walk along a 5m straight line, turn right (or left), walk along a 5m straight line again, turn 180 • and return to the starting point. Hence the total traveled distance is 20m. The inverted L-shaped trajectory test is considered because walking with turning is a common activity in our daily life. Turning frequently causes loss of balance in all gait disorders [START_REF] Joseph | Gait in the elderly[END_REF], and it is associated with increased risk of falling in elderly adults [START_REF] Tinetti | Risk factors for falls among elderly persons living in the community[END_REF], [START_REF] Tinetti | Fall risk index for elderly patients based on number of chronic disabilities[END_REF].

Many analysis of turning are considered in TUG test or figure-of-eight walking test. These studies examined the usual gait parameters such as used time [START_REF] Beauchet | Imagined timed up & go test: a new tool to assess higher-level gait and balance disorders in older adults?[END_REF], gait speed, stride length, cadence, double-limb support duration during the turning [START_REF] Shkuratova | Effects of age on balance control during walking[END_REF]. Other methods assess accelerometer-derived parameters [START_REF] Weiss | Can an accelerometer enhance the utility of the timed up & go test when evaluating patients with parkinson's disease?[END_REF], [START_REF] Zakaria | Quantitative analysis of fall risk using tug test[END_REF] and gyroscope-derived parameters [START_REF] Greene | Quantitative falls risk assessment using the timed up and go test[END_REF]. Our objective is to find new quality index during turning phase by using the walker instrumented with an accelerometer/gyrometer, which may help doctors to evaluate the fall risk of the elderly.

Detection of distinct phases

Since the inverted L-shaped trajectory test consists of a sequence of activities such as walking, 90 • turning and 180 • turning, it is necessary to evaluate the user's ability in distinct phases or activities to make the test more sensitive. Several studies use wearable inertial sensors to detect the phases in TUG test by using the Euler angles or angular velocity signals that had been processed by filters [START_REF] Higashi | Quantitative evaluation of movement using the timed up-and-go test[END_REF], [START_REF] Zakaria | Quantitative analysis of fall risk using tug test[END_REF]. Some specific algorithms used to filter out the noise of the sensors for the detection are presented in [START_REF] Salarian | Analyzing 180 degrees turns using an inertial system reveals early signs of progression of parkinson's disease[END_REF], [START_REF] Al-Jawad | Using multi-dimensional dynamic time warping for tug test instrumentation with inertial sensors[END_REF]. Here we use filtered yaw angle and yaw angular velocity to detect each activity and the test was divided into seven phases in sequence: walk 1, turn 1 (walking with 90 • turning), walk 2, turn 2 (walking with 180 • turning), walk 3, turn 3 (walking with 90 • turning) and walk 4, as shown in Fig. 11 and Fig. 12. For example, Fig. 11 presents both yaw angle and yaw angular velocity signals of an elderly subject during the test.

In order to follow the reference trajectory shown in Fig. 12, the yaw angle of the walker changes from 0 • to a nearby value of -90 • , 90 • and 180 • respectively at the end of turn 1, turn 2 and turn 3. When the yaw angular velocity passes zero and the yaw angle signal starts to increase or decrease suddenly, the turn phase starts. From this moment when the yaw angular velocity passes zero again, the turn phase ends.

As soon as the onset and offset of the three turn phases are detected, the duration of all the seven phases are obtained. During the four walk phases, the yaw angular velocity passes zero several times although the yaw angle almost keeps a constant value, and these zero-pass moments can be used to detect the user's steps, which has been presented in the previous section.

Results

There were 25 young adults (between 25 and 65 years, mean value 32) and 12 elderly people (> 69 years) that participated to the inverted L-shaped trajectory test. The walking trajectories of the two groups using the walker are presented in Fig. 13, which clearly shows that most of the elder subjects cannot follow exactly the reference trajectory. This is consistent with the result of 10mWT (see Fig. 8 and Fig. 9), meaning that the deviation between the reference trajectory and the followed trajectory can be always regarded as a walking indicator.

Turning duration is always an important index to be examined in the TUG test. It has been claimed that the elderly or people that has higher fall risk use more time during turning [START_REF] Higashi | Quantitative evaluation of movement using the timed up-and-go test[END_REF], [START_REF] Zakaria | Quantitative analysis of fall risk using tug test[END_REF], [START_REF] Salarian | Analyzing 180 degrees turns using an inertial system reveals early signs of progression of parkinson's disease[END_REF]. Fig. 14 presents the duration of each phase for the two groups. It clearly shows that the duration of all turning phases for the elderly is longer than that of the younger group, which is consistent with the result of previous studies. However, during the walking phases without turning, the elderly use less time. One reason which could explain this is shorter walking distance of the elderly (see Fig. 13) , with similar walking velocity of two groups, as presented in the previous section.

Using only the time parameter to evaluate the fall risk of the subject is obviously not sufficient. To further investigate this issue, other gait parameters derived from the walker will be assessed. From the walking trajectories shown in Fig. 13, we can see that the elderly need larger space to complete a turn. Hence it is interesting to compare the maneuver area, which is surrounded by a turning trajectory and a straight line that connects the start and end of the turning, as shown in green shadings of Fig. 12. Fig. 15 indicates that the maneuver area of the elderly is significantly larger. Hence it can be regarded as an important indicator to quantify the maneuvering ability of subjects. Root mean square (RMS) is a statistical measure of the magnitude of a varying quantity and is especially useful for sinusoids. We compared the RMS value of yaw angular velocity during turning phases between the two groups (Fig. 16). It clearly shows that during turning phases the young people has a higher RMS value of yaw angular velocity. The value of every subject during the turn around phase (turn 2) is presented in Fig. 17. It indicates that almost all members of the elderly group has lower RMS value of yaw angular velocity except for one member (92.54 • /s). In addition, the trajectory of this subject, presented in Fig. 18, shows that the maneuver area during turn 2 phase is very small, and almost equal to the one of the young people. In order to further investigate this issue, Fig. 19 presents the maneuver area and yaw angular velocity during turn 2 phase for all subjects. It shows that during the turning phase the young people has high RMS value of yaw angular velocity and a small maneuver area. This indicates that the two groups can almost be discriminated only by using these two indicators. Moreover, for the one subject that has the maximum RMS value of yaw angular velocity (92.54 • /s) in the elderly group, the maneuver area is only 857.4 cm 2 . Thus in Fig. 19 this point can be classified into the young group, which can explain the abnormality of Fig. 17 during each phase is compared for the two groups in Fig. 20. It reveals:

all the subjects exhibit a larger acceleration during turning phases than during walking phases, -for the elderly, the difference between the acceleration during turning phase and walking phase is significantly smaller, -the elderly have smaller acceleration than the young people during the turn around phase while they exhibit a larger acceleration during walking phases. In summary, by using the smart walker instrumented with inertial sensors, the subject's gait can be assessed by the above indicators. Besides, we also compared other parameters for the two groups. Derived indicators during the 10 m straight line test and the turn around phase of inverted L-shaped trajectory test are summarized in Table 1.

Conclusions

This paper proposed a gait analysis method an instrumented walker. The results of a 10 m straight line test and an inverted L-shaped trajectory test for a group of young adults and a group of elderly people were compared comprehensively. The gait of the two groups were analyzed during both walking phases and turning phases. Several indicators that exhibit significant difference between the two groups were obtained, such as the maximum lateral deviations, the orientation angle and the traveled Manhattan distance. For the elderly these indicators are much larger than that of the young people, which indicates that the walking accuracy of the elderly is lower. Moreover, we found there is no obvious difference in step length, step period and walking speed between two groups, which is consistent with some previous studies. It reveals that these classical gait indicators are neither sufficient nor sensitive enough to evaluate the fall risk of the elderly.

For the inverted L-shaped trajectory test we found that indicators derived form the walker's sensors can be used to classify the gait quality of subjects. For example, during the turning phases for the elderly the maneuver area is larger and the RMS value of angular velocity is smaller than that of the young people. In addition, for the elderly during turning phases the RMS value of acceleration is lower than that of the young people while during the walking phases it is larger than that of young people. Surprisingly this is inconsistent with previous studies [START_REF] Higashi | Quantitative evaluation of movement using the timed up-and-go test[END_REF], [START_REF] Zakaria | Quantitative analysis of fall risk using tug test[END_REF], which claimed that during both walking phases and turning phases the RMS values of acceleration for the group with high fall risk are smaller than that of the group with lower fall risk. However, during the 10mWT we also found that the elderly have larger forward acceleration and the walking velocity of the two groups are similar. Our interpretation is that elderly people are more familiar with walking aids and have walking patterns that benefit from such an aid, while younger people have a more dynamic pattern that is jeopardized by the aid. An open issue is to deduce the gait indicators for natural gait without the walker from the measurements obtained with the walker.

FITFig. 1 . 1 Fig. 2 .

 112 Fig. 1. The instrumented walking aid.

Fig. 3 .

 3 Fig. 3. Yaw angular velocity of the walker when the subject walks along a straight reference trajectory with 10m. One step is finished when it passes zero.

Fig. 4 .

 4 Fig. 4. Displacement of the walker during every step. The results of left steps and right steps were put together and they appeared alternately.

Fig. 5 .

 5 Fig.5. Walker and motion capture optical markers set-up[START_REF] Wu | {ISB} recommendations for standardization in the reporting of kinematic data[END_REF]

Fig. 8 .

 8 Fig. 8. Trajectory of the subjects in the xy plane, where the blue color denotes the elderly subjects and the red one denotes the young subjects. The reference trajectory is the horizontal line: y = 0.

Fig. 9 .

 9 Fig. 9. Maximum lateral deviation between the real and reference trajectory. The results are sorted in ascending order.

Fig. 10 .

 10 Fig. 10. Maximum instantaneous walking velocity of the two groups. The results are sorted in ascending order.

Fig. 11 .

 11 Fig. 11. Yaw angle and yaw angular velocity of an elderly subject during the inverted L-shaped trajectory test.

Fig. 12 .

 12 Fig. 12. Inverted L-shaped reference trajectory.

Fig. 13 .

 13 Fig. 13. Trajectory of the elderly (blue color) and young people (red color) in the xy plane.

Fig. 14 .

 14 Fig. 14. Comparison of the duration of each phase for young and elderly people.

Fig. 15 .

 15 Fig. 15. Comparison of the maneuver area of each turning phase for the two groups.

Fig. 16 .

 16 Fig. 16. Comparison of the RMS value of yaw angular velocity for the two groups.

Fig. 17 .Fig. 18 .

 1718 Fig. 17. RMS value of yaw angular velocity during the turn 2 phase for the two groups. The results are sorted in ascending order.

2 Fig. 19 .

 219 Fig. 19. RMS value of yaw angular velocity as a function of the maneuver area during the turn around phase (turn 2).

Fig. 20 .

 20 Fig. 20. RMS value of acceleration during each phase for young and elderly people.

Table 1 .

 1 Gait indicators obtained during 10mWT and turn around phase of the inverted L-shaped trajectory test by using the walker, where S means small and L means large.

	10m straight trajectory			Inverted L-shaped trajectory
	Indicators	Old Young	Indicators	Old Young
	Lateral deviation	L	S	Turning duration	L	S
	Manhattan distance	L	S	Maneuver area	L	S
	Orientation angle	L	S	RMS value of angular velocity	S	L
	Forward acceleration	L	S	Amplitude of angular velocity	S	L
	Percentage of forward support	S	L	RMS value of acceleration	S	L


Alexandru Amȃrioarei
email: alexandru.amarioarei@inria.fr

Cristian Preda
email: cristian.preda@polytech-lille.fr

Approximation for the Distribution of Three-dimensional Discrete Scan Statistic

Introduction

Let T 1 , T 2 , T 3 be positive integers, R = [0, T 1 ] × [0, T 2 ] × [0, T 3 ] be a rectangular region and {X ijk |1 ≤ i ≤ T 1 , 1 ≤ j ≤ T 2 , 1 ≤ k ≤ T 3 } be a family of independent and identically distributed integer valued random variables from a specified distribution. In practice, X ijk can be interpreted as the number of events that occur in the elementary subregion

r ijk = [i-1, i]×[j -1, j]×[k -1, k].
For each j ∈ {1, 2, 3}, consider the positive integers m j such that 2 ≤ m j ≤ T j -1, and define the random variables

Y i1i2i3 = i1+m1-1 i=i1 i2+m2-1 j=i2 i3+m3-1 k=i3 X ijk , 1 ≤ i j ≤ T j -m j + 1, (1.1)
as the number of events occurring in the rectangular region

R(i 1 , i 2 , i 3 ) = [i 1 -1, i 1 + m 1 -1] × [i 2 -1, i 2 + m 2 -1] × [i 3 -1, i 3 + m 3 -1].
The three dimensional discrete scan statistic is defined as the maximum number of events in any rectangle R(i 1 , i 2 , i 3 ) within the region R, S m1,m2,m3 (T 1 , T 2 , T 3 ) = max 1≤ij ≤Tj -mj +1 j∈{1,2,3} Y i1i2i3 .

(1.

2)

The distribution of scan statistics, P (S m1,m2,m3 (T 1 , T 2 , T 3 ) ≤ n) , n ∈ {1, 2, . . . , m 1 m 2 m 3 } is used with success in astronomy [START_REF] Darling | Approximations for three dimensional scan statistic[END_REF]), image analysis and reliability theory [START_REF] Boutsikas | Reliability approximations for Markov chain imbeddable systems[END_REF]) and many other domains.

For an overview of the potential application of scan statistics one can refer to the monographs of [START_REF] Glaz | Scan statistic[END_REF] and more recently the one of [START_REF] Glaz | Scan statistic: Methods and Applications[END_REF].

From a statistical point of view, the scan statistic S m1,m2,m3 (T 1 , T 2 , T 3 ) is used for testing the null hypothesis of randomness that X ijk 's are independent and identically distributed according to some specified distribution. Under the alternative hypothesis there exists one cluster location where the X ijk 's have a larger mean than outside the cluster. As an example, in the Poisson model, the null hypothesis, H 0 , assumes that X ijk 's are i.i.d. with X ijk ∼ P ois(λ) whereas the alternative hypothesis of clustering, H 1 , assumes the existence of a rectangular subregion R(i 0 , j 0 , k 0 ) such that for any

i 0 ≤ i ≤ i 0 +m 1 -1, j 0 ≤ j ≤ j 0 +m 2 -1 and k 0 ≤ k ≤ k 0 +m 3 -1, X ijk are i.i.d.
Poisson random variables with parameter λ ′ > λ. Outside the region R(i 0 , j 0 , k 0 ), X ijk are i.i.d. distributed according to the distribution specified by the null hypothesis. The generalized likelihood ratio test rejects H 0 in favor of the local change alternative H 1 , whenever S m1,m2,m3 (T 1 , T 2 , T 3 ) exceeds the threshold τ determined from P (S m1,m2,m3 (T 1 , T 2 , T 3 ) ≥ τ |H 0 ) = α and where α represents 1 the significance level of the testing procedure (Glaz, Naus and Wallenstein [2001, Chapter 13]). Since there are no exact formulas available for the distribution of three dimensional scan statistics, approximation methods are necessary. For the Bernoulli model, [START_REF] Glaz | Approximations for three dimensional scan statistic[END_REF] propose four approximation formulas: one Markov like product type approximation and three Poisson type approximations that extends the special case described by [START_REF] Darling | Approximations for three dimensional scan statistic[END_REF] when n = m 1 m 2 m 3 . The advantage of the method described in this paper is that it can be used for any distribution of the random field and provides accurate approximations and sharp error bounds. The methodology used to obtain the approximation and the error bounds is presented in Section 2. In Section 3 we describe adapt the importance sampling algorithm developed by [START_REF] Naiman | Computing Scan Statistic p Values Using Importance Sampling, with Applications to Genetics and Medical Image Analysis[END_REF] to estimate the simulation error. A simulation study is conducted in Section 4 for considered Bernoulli, binomial and Poisson models. Concluding remarks are given in Section 5.

Methodology

In order to approximate the distribution of S m1,m2,m3 (T 1 , T 2 , T 3 ) we use a similar approach as in [START_REF] Haiman | Estimation for the distribution of two-dimensional scan statistics[END_REF]. The key idea is to observe that we can write the scan statistic random variable as a maximum of 1-dependent stationary sequence of random variables. A sequence (Z k ) k≥1 is m-dependent, m ≥ 1, if for any h ≥ 1 the σ-fields generated by {Z 1 , . . . , Z h } and {Z h+m+1 , . . . } are independent.

The method is based on the following result developed in Haiman [1999, Theorem 4] and improved in Amarioarei [2012, Theorem 2.6]: Let (Z k ) k≥1 be a strictly stationary 1-dependent sequence of random variables and for x < sup{u|P(Z 1 ≤ u) < 1}, let q m = q m (x) = P(max(Z 1 , . . . , Z m ) ≤ x).

(2.1)

Theorem 2.1. For all x such that q 1 (x) ≥ 1-α ≥ 0.9, the following approximation formula holds:

q m - 2q 1 -q 2 [1 + q 1 -q 2 + 2(q 1 -q 2 ) 2 ] m ≤ mF (α, m)(1 -q 1 ) 2 (2.2) with F (α, m) = 1 + 3 m + Γ(α) m + K(α) (1 -q 1 ) (2.3) where Γ(α) = L(α) + E(α), K(α) = 11-3α (1-α) 2 + 2l(1 + 3α) 2+3lα-α(2-lα)(1+lα) 2 [1-α(1+lα) 2 ] 3 1 -2α(1+lα) [1-α(1+lα) 2 ] 2
(2.4)

L(α) = 3K(α)(1 + α + 3α 2 )[1 + α + 3α 2 + K(α)α 3 ] + α 6 K 3 (α) + 9α(4 + 3α + 3α 2 ) + 55.1 (2.5) E(α) = η 5 [1 + (1 -2α)η] 4 [1 + α(η -2)] 1 + η + (1 -3α)η 2 2(1 -αη 2 ) 4 [(1 -αη 2 ) 2 -αη 2 (1 + η -2αη) 2 ] (2.6)
and where η = 1 + lα with l = l(α) > t 3 2 (α) and t 2 (α) the second root in magnitude of the equation αt 3t + 1 = 0.

In this section we obtain an approximation formula for the distribution of scan statistic defined by Eq.(1.2) in three steps as follows.

Let assume that L j = Tj mj -1 , j ∈ {1, 2, 3}, are positive integers and define for each k ∈ {1, 2, . . . , L 3 -1} the random variables

Z k = max 1≤i1≤(L1-1)(m1-1) 1≤i2≤(L2-1)(m2-1) (k-1)(m3-1)+1≤i3≤k(m3-1) Y i1i2i3 .
(2.7)

The set of random variables {Z 1 , . . . , Z L3-1 } forms a 1-dependent stationary sequence. Indeed, from Eq.(2.7) and the independence of X ijl we observe that for

any k ≥ 1, σ(• • • , Z k ) and σ(Z k+2 , • • • ) are included in σ({X ijl |1 ≤ i ≤ T 1 , 1 ≤ j ≤ T 2 , 1 ≤ l ≤ (k+1)(m 3 -1)}) and σ({X ijl |1 ≤ i ≤ T 1 , 1 ≤ j ≤ T 2 , (k+1)(m 3 -1)+1 ≤ l})
, respectively, which are independent (see Fig. 1).

T 2 T 1 T 3 1 2(m3 -1) m3
3(m3 -1) 2m3 -1 4(m3 -1) Z3 Z2 Z1 m1 m2 m3 Figure 1. Illustration of Z k emphasizing the 1-dependence
Notice that from Eq.(1.2) and Eq.(2.7) we have

S(L 1 , L 2 , L 3 ) = S m1,m2,m3 (T 1 , T 2 , T 3 ) = max 1≤k≤L3-1 Z k .
(2.8) Take for s ∈ {2, 3},

Q s = Q s (n) = P s-1 k=1 {Z k ≤ n} = P      max 1≤i1≤(L1-1)(m1-1) 1≤i2≤(L2-1)(m2-1) 1≤i3≤(s-1)(m3-1) Y i1i2i3 ≤ n      . (2.9)
Notice that in the notation of Eq.(2.1) we have

Q s = q s-1 . For n such that Q 2 (n) ≥ 1 -α 1 ≥ 0.
9 we apply the result in Theorem 2.1 to obtain the first step approximation

P (S(L 1 , L 2 , L 3 ) ≤ n) ≈ 2Q 2 -Q 3 [1 + Q 2 -Q 3 + 2(Q 2 -Q 3 ) 2 ] (L3-1) , (2.10) with an error bound of (L 3 -1)F (α 1 , L 3 -1)(1 -Q 2 ) 2 .
Observe that Q 2 and Q 3 represents the distribution of the scan statistics over the rectangular subregions 1, 2(m 3 -1)] and[1, T 1 ] × [1, T 2 ] × [1, 3(m 3 -1)], respectively (see also Fig 1). To simplify the results of the presentation, in what follows we abbreviate the approximation formula by

[1, T 1 ] × [1, T 2 ] × [
H(x, y, m) = 2x -y [1 + x -y + 2(x -y) 2 ] m-1 .
(2.11)

In order to evaluate the approximation in Eq.(2.10) it is necessary to find approximations for Q 2 and Q 3 . Thus, the second step consists in applying Theorem 2.1 for each Q s . We define, as in Eq.(2.7), for s ∈ {2, 3} and l ∈ {1, 2, . . . , L 2 -1} the sequences Z

(s) l = max 1≤i1≤(L1-1)(m1-1) (l-1)(m2-1)+1≤i2≤l(m2-1) 1≤i3≤(s-1)(m3-1) Y i1i2i3 , (2.12)
which are strictly stationary, 1-dependent and satisfy

Q s = P(S(L 1 , L 2 , s) ≤ n) = P max 1≤l≤L2-1 Z s l ≤ n . (2.13) Set for t, s ∈ {2, 3}, Q ts = Q ts (n) = P t-1 l=1 {Z (s) l ≤ n} = P      max 1≤i1≤(L1-1)(m1-1) 1≤i2≤(t-1)(m2-1) 1≤i3≤(s-1)(m3-1) Y i1i2i3 ≤ n      . (2.14) If the condition Q 2s (n) ≥ 1 -α 2s ≥ 0.
9 is fulfilled, then using Theorem 2.1, we find, for s ∈ {2, 3}, the approximation

|Q s -H (Q 2s , Q 3s , L 2 )| ≤ (L 2 -1)F (α 2s , L 2 -1)(1 -Q 2s ) 2 .
(2.15)

The last step involves the evaluation of Q ts in Eq.(2.15). For s, t ∈ {2, 3} and j ∈ {1, 2, . . . , L 1 -1} let consider the following random sequences

Z (ts) j = max (j-1)(m1-1)+1≤i1≤j(m1-1) 1≤i2≤(t-1)(m2-1) 1≤i3≤(s-1)(m3-1) Y i1i2i3 .
(2.16)

We observe that Z (ts) j j≥1

forms 1-dependent stationary sequences and

Q ts = P (S(L 1 , t, s) ≤ n) = P max 1≤j≤L1-1 Z (ts) j ≤ n .
(2.17)

Put for r, t, s ∈ {2, 3}

Q rts = Q rts (n) = P   r-1 j=1 {Z ts j ≤ n}   = P      max 1≤i1≤(r-1)(m1-1) 1≤.i2≤(t-1)(m2-1) 1≤i3≤(s-1)(m3-1) Y i1i2i3 ≤ n      (2.18)
Then, under supplementary condition that Q 2ts ≥ 1α 3ts ≥ 0.9, we apply the result in Theorem 2.1 to obtain

|Q ts -H (Q 2ts , Q 3ts , L 1 )| ≤ (L 1 -1)F (α 3ts , L 1 -1)(1 -Q 2ts ) 2 . ( 2
P(S(L 1 + 1, L 2 + 1, L 3 + 1) ≤ n) ≤ P ≤ P(S(L 1 , L 2 , L 3 ) ≤ n),
(2.20)

we can approximate P = P(S m1,m2,m3 (T 1 , T 2 , T 3 ) ≤ n) by linear interpolation (see Table 3).

2.1.

Computing the approximation error. To simplify the presentation and the derivation of the approximation formulae, it is convenient to introduce the following notations for s, t ∈ {2, 3}:

α 3 = 1 -Q 3 , α 23 = 1 -Q 23 , α 233 = 1 -Q 233 , γ ts = H(Q 2ts , Q 3ts , L 1 ), γ s = H(γ 2s , γ 3s , L 2 ), F 1 = F (α 3 , L 3 -1), F 2 = F (α 23 , L 2 -1), F 3 = F (α 233 , L 1 -1). It is not hard to see that Q 3 ≤ Q 2 , Q 23 ≤ Q 22 and Q 233 ≤ Q 2ts
, so that the choice for the thresholds α 3 , α 23 and α 233 becomes natural. Based on the mean value theorem in two dimensions, one can easily verify that for m ≥ 6 and y i ≤ x i , i ∈ {1, 2} we have the inequality:

|H(x 1 , y 1 , m) -H(x 2 , y 2 , m)| ≤ (m -2) [|x 1 -x 2 | + |y 1 -y 2 |] . (2.21)
In what follows we use the result from Eq.(2.21) without restrictions. This is in agreement with the numerical values considered in Section 4. We begin by observing that applying Eq.(2.21) into Eq.(2.10) we obtain

|P -H (γ 2 , γ 3 , L 3 )| ≤ |P -H (Q 2 , Q 3 , L 3 )| + |H (Q 2 , Q 3 , L 3 ) -H (γ 2 , γ 3 , L 3 )| ≤ (L 3 -1)F 1 (1 -Q 2 ) 2 + (L 3 -2) [|Q 2 -γ 2 | + |Q 3 -γ 3 |] , (2.22)
where for simplicity we used the notation P = P (S(L 1 , L 2 , L 3 ) ≤ n). In the same manner, one can see that for s ∈ {2, 3} we have

|Q s -γ s | ≤ |Q s -H (Q 2s , Q 3s , L 2 )| + |H (Q 2s , Q 3s , L 2 ) -H (γ 2s , γ 3s , L 2 )| ≤ (L 2 -1)F 2 (1 -Q 2s ) 2 + (L 2 -2) [|Q 2s -γ 2s | + |Q 3s -γ 3s |] . (2.23)
We notice that Eq.(2.19) can be rewritten as

|Q ts -γ ts | ≤ (L 1 -1)F 3 (1 -Q 2ts ) 2 , s, t ∈ {2, 3}. (2.24)
Finally, in order to find the approximation error it is sufficient to determine bounds for 1 -Q 2 and 1 -Q 2s . It can be easily checked that

1 -Q 2s ≤ 1 -γ 2s + |Q 2s -γ 2s | ≤ δ 2s (2.25)
where

δ 2s = 1 -γ 2s + (L 1 -1)F 3 (1 -Q 22s ) 2 .
(2.26)

Similarly, we can write

1 -Q 2 ≤ 1 -γ 2 + |Q 2 -γ 2 | ≤ δ 2 , (2.27) with δ 2 = 1-γ 2 +(L 2 -1)F 2 δ 22 +(L 2 -2)(L 1 -1)F 3 (1 -Q 222 ) 2 + (1 -Q 232 ) 2 . (2.
E app = (L 3 -1)F 1 δ 2 2 + (L 3 -2)(L 2 -1)F 2 δ 2 22 + δ 2 23 + + (L 3 -2)(L 2 -2)(L 1 -1)F 3   t,s∈{2,3} (1 -Q 2ts ) 2   . ( 2
H (γ 2 , γ 3 , L 3 ) -H Q2 , Q3 , L 3 ≤ (L 3 -2) γ 2 -Q2 + γ 3 -Q3 . (2.31)
Observe that the differences in the right hand term in Eq.(2.31) can be bounded by

γ s -Qs = H (γ 2s , γ 3s , L 2 ) -H Q2s , Q3s , L 2 ≤ (L 2 -2) γ 2s -Q2s + γ 3s -Q3s .
(2.32)

In the same way we can write for t, s ∈ {2, 3}

γ ts -Qts = H (γ 2ts , γ 3ts , L 1 ) -H Q2ts , Q3ts , L 1 ≤ (L 1 -2) γ 2ts -Q2ts + γ 3ts -Q3ts ≤ (L 1 -2) [β 2ts + β 3ts ] .
(2.33)

Combining Eqs.(2.33), (2.32) and (2.31) we get the simulation error corresponding to the approximation formula

E sf = (L 1 -2)(L 2 -2)(L 3 -2)   r,t,s∈{2,3} β rts   .
(2.34)

In order to obtain the simulation error corresponding to the approximation error bound in Eq.(2.29) we follow the lines of Section 2.1. With the following notations

u rts = 1 -qrts + β rts , u ts = 1 -qts + (L 1 -2)(β 2ts + β 3ts ), u s = 1 -qs + (L 1 -2)(L 2 -2)(β 22s + β 32s + β 23s + β 33s ), δ2s = u 2s + (L 1 -1)F 3 u 2 22s , δ2 = u 2 + (L 2 -1)F 2 δ22 + (L 2 -2)(L 1 -1)F 3 (u 2 222 + u 2 232 ),
the error can be expressed as

E sapp = (L 3 -1)F 1 δ2 2 + (L 3 -2)(L 2 -1)F 2 δ2 22 + δ2 23 + + (L 3 -2)(L 2 -2)(L 1 -1)F 3   t,s∈{2,3} u 2 2ts   . (2.35)
The total simulation error is obtained by adding the two terms from Eq.(2.34) and Eq.(2.35)

E sim = E sf + E sapp .
(2.36)

To evaluate Eq.( 2.36), one needs to find suitable values for the bounds β rts . If IT ER is the number of iterations used in the Monte Carlo simulation algorithm for the estimation of Q rts then, one can consider, for example, the naive bound provided by the Central Limit Theorem with a 95% confidence level

β rts = 1.96 Qrts (1 -Qrts )
IT ER .

(2.37)

This bound has been used with some success for the two dimensional case (see [START_REF] Haiman | Estimation for the distribution of two-dimensional scan statistics[END_REF]). As the authors pointed out, the main contribution to the total error is due to the simulation error, especially for small sizes of the window scan with respect to the scanning region. Our numerical study shows that Eq.(2.37) is not feasible for the three dimensional case, the simulation error being to large with respect to the approximation error. Thus, for the simulation of Qrts , we use an importance sampling technique introduced in Naiman and Priebe [2001]. Next section illustrates how to adapt theirs algorithm to our problem.

Simulation by importance sampling

In this section we present a simulation method for Q rts , which gives an unbiased estimate whose variance is typically smaller then that of the naive hit or miss Monte Carlo approach. The method is an adaptation of the importance sampling algorithm developed in [START_REF] Naiman | Computing Scan Statistic p Values Using Importance Sampling, with Applications to Genetics and Medical Image Analysis[END_REF] to our problem. The main idea behind is to express the tail of the scan distribution as a Bonferroni upper bound (B) with some correction factor (ρ). Let define for 1 ≤ i j ≤ N j , j ∈ {1, 2, 3} the events

A i1i2i3 = {Y i1i2i3 ≥ τ }. Then P (S m1,m2,m3 (T 1 , T 2 , T 3 ) ≥ τ ) = P T1-m1+1 i1=1 T2-m2+1 i2=1 T3-m3+1 i3=1 A i1i2i3 = B T1-m1+1 i1=1 T2-m2+1 i2=1 T3-m3+1 i3=1 p i1i2i3 I(i 1 , i 2 , i 3 ) = Bρ, (3.1)
where

ρ = T1-m1+1 i1=1 T2-m2+1 i2=1 T3-m3+1 i3=1 p i1i2i3 I(i 1 , i 2 , i 3 ). (3.2)
Under the null hypothesis (H 0 ), B is the Bonferroni upper bound given by

B = T1-m1+1 i1=1 T2-m2+1 i2=1 T3-m3+1 i3=1 P(A i1i2i3 ) = (T 1 -m 1 + 1)(T 2 -m 2 + 1)(T 3 -m 3 + 1)P(A 111 ), (3.3) p i1i2i3 defines an uniform probability distribution over {1, . . . , T 1 -m 1 + 1} × {1, . . . , T 2 -m 2 + 1} × {1, . . . , T 3 -m 3 + 1}, p i1i2i3 = P(A i1i2i3 ) T1-m1+1 s1=1 T2-m2+1 s2=1 T3-m3+1 s3=1 P(A s1s2s3 ) = 1 (T 1 -m 1 + 1)(T 2 -m 2 + 1)(T 3 -m 3 + 1) , (3.4)
and

I(i 1 , i 2 , i 3 ) = 1 C(Y ) 1 Ai 1 i 2 i 3 P(A i1i2i3 )
dP where C(Y ) represents the number of triples (i 1 , i 2 , i 3 ) such that Y i1i2i3 exceeds the threshold τ , that is

C(Y ) = T1-m1+1 i1=1 T2-m2+1 i2=1 T3-m3+1 i3=1 1 Ai 1 i 2 i 3 . (3.5)
Based on these identities the simulation algorithm (similar with the one described in Naiman and Priebe [2001, page 303]) can be written as follows:

Begin Repeat for each k from 1 to IT ER (iterations number)

Step 1 Generate T ∈ {τ, . . . } according to the probabilities

p T (t) = P(Y 111 = t) s≥τ P(Y 111 = s) , t ≥ τ.
Step 2 Conditionally, given T = t, generate the triple (J 1 , J 2 , J 3 ) in the set {1, . . . , T 1 -m 1 +1}×{1, . . . , T 2 -m 2 +1}×{1, . . . , T 3 -m 3 +1} uniformly.

Step 3 Conditionally, given T and (J 1 , J 2 , J 3 ), generate the set of random variables { Ỹi1i2i3 |J s ≤ i s ≤ J s + m s -1, s ∈ {1, 2, 3}}, uniformly from the set of all the vectors of length m 1 × m 2 × m 3 over the set of values taken by Y i1i2i3 and whose sum is equal with T . Take the remaining Ỹi1i2i3 to be i.i.d. and distributed according to the null hypothesis distribution.

Step 4 Take c k = C( Ỹk ) the number of all triples (i 1 , i 2 , i 3 ) such that Ỹi1i2i3 ≥ T and put ρk =

1 c k . End Repeat Return ρ = 1 IT ER IT ER k=1 ρk .

End

Clearly, ρ is an unbiased estimator for ρ with estimated variance

V ar(ρ) ≈ 1 IT ER -1 IT ER k=1 ρk - 1 IT ER IT ER k=1 ρk 2 .
(3.6)

For IT ER sufficiently large, as a consequence of CLT the error between the true and the estimated value of the tail P (S m1,m2,m3 (T 1 , T 2 , T 3 ) ≥ τ ), corresponding to a 95% confidence level, is given by

β = 1.96B V ar(ρ) IT ER . (3.7)
Notice that for the simulation of Q rts , we substitute T 1 , T 2 and T 3 in the above relations with r(m 1 -1), t(m 2 -1) and s(m 3 -1) respectively. Therefore, we obtain the corresponding values for β rts as described by Eq.(3.7).

Numerical values for Binomial and Poisson models

In this section, for selected values of the parameters of the binomial and Poisson distributions, we evaluate the approximation introduced in Section 2 and provide the corresponding error bounds. We show the contributions of the approximation (Eq.(2.29)) and simulation (Eq.(2.36)) errors in the overall error. For all our simulations we used the importance sampling algorithm with IT ER = 10 5 replications. We compare our results with those existing in literature, see [START_REF] Glaz | Approximations for three dimensional scan statistic[END_REF] for Bernoulli model, and with the simulated value of the scan statistics obtained by scanning the whole region R, denoted by P(S ≤ n).

The scanning of R being more time consuming than the scanning of the subregions corresponding to Q rst , we used 10 3 repetitions of the algorithm.

In Table 1, we compare the results obtained by our approximation with the product type approximation presented by [START_REF] Glaz | Approximations for three dimensional scan statistic[END_REF]. We observe that our approximation is very sharp. In order to compare the binomial and Poisson models, in Table 4, we have evaluated the distribution of the scan statistics over a region of size 84 × 84 × 84 scanned with a 4 × 4 × 4 cubic window, in the two situations. In the first case we have a binomial random field with parameters m and p, that is X ijk ∼ B(m, p), while in the second we considered that X ijk ∼ P (λ), with λ = mp.

Notice that the contribution of the approximation error (E app ) to the total error is almost negligible in most of the cases with respect to the simulation error (E sim ). Thus, the precision of the method will depend mostly on the number of iterations (IT ER) used to estimate Q rts .

= m2 = m3 = 10, T1 = T2 = T3 = 185, L1 = L2 = L3 = 20, IT ER = 10 5 n P (S(L 1 + 1, L 2 + 1, L 3 + 1) ≤ n) P(S ≤ n) P (S(L 1 , L 2 , L 3 ) ≤ n) 4
The time required for the computations presented in this section was about two hours for each table on a computer of medium performances. The programs are written in MATLAB and are available from the authors.

Conclusions

In this article we derived an approximation for the three dimensional discrete scan statistic viewed as the maximum of a 1-dependent stationary sequence of random variables. We also provide the corresponding theoretical and simulation error bounds. In the three dimensional scan statistics framework, it is essential to reduce the variance of simulated values. For this purpose we used an importance sampling method. A simulation study for the binomial and Poisson models shows the accuracy as well as the limit of our method.

  In the next section we show how to choose these values.Remark 2.2. If T 1 , T 2 and T 3 are not multiples of m 1 -1, m 2 -1 and m 3 -1, respectively, then let consider L j =

	Tj mj -1	for j ∈ {1, 2, 3}. Based on the
	inequalities

.19)

Combining the Eqs.(2.10), (2.15) and (2.19) we obtain an approximation formula for the distribution of the scan statistic depending on the eight quantities Q rts , that we propose to evaluated by simulation. Note that in the above approximations, at each step we consider different values for α.

  28) Substituting Eqs.(2.23), (2.24), (2.25) and (2.27) in Eq.(2.22) we derive the formula for the approximation error

  .29) 2.2. Computing the simulation errors. Since, from our knowledge, there are no exact formulas available for the computation of Q rts we propose to evaluate them by simulation. It is obvious that the simulation error appears from two terms: first, from the approximation formula in Eq.(2.22) and second, from the error bound in Eq.(2.29). Usually, between the true and the estimated value we have a relation of the form Q

rts -Qrts ≤ β rts , r, t, s ∈ {2, 3} (2.30) where Qrts are the simulated values corresponding to Q rts . Provided a simulation error bound β rts as in Eq.(2.30), let denote the simulated values by Qts = H( Q2ts , Q3ts , L 1 ), Qs = H( Q2s , Q3s , L 2 ). From Eq.(2.21) one obtains

Table 1 .

 1 Approximation for P(S ≤ n) in Bernoulli case: m1 = m2 = m3 = 5, T1 = T2 = T3 = 60, IT ER = 10 5

	n P(S ≤ n)	Glaz et al.	Our	Eapp	E sim	Total
		Product type Approximation Eq.(2.29) Eq.(2.36)	Error
			p = 0.00005
	1 0.841806	0.841424	0.851076	0.011849	0.064889 0.076738
	2 0.999119	0.999142	0.999192	0.000000	0.000170 0.000170
	3 0.999997	0.999998	0.999997	0.000000 3 × 10 -7 3 × 10 -7
			p = 0.0001
	2 0.993294	0.993241	0.993192	0.000010	0.001367 0.001377
	3 0.999963	0.999964	0.999963	0.000000	0.000005 0.000005
	4 0.999999	0.999999	0.999999	0.000000 2 × 10 -9 2 × 10 -9

  Table2presents the numerical results obtained by scanning the region R of size 60 × 60 × 60 with two windows of the same volume but different sizes, first a cubic window of size 4 × 4 × 4 and second a rectangular region of size 8 × 4 × 2. We observe that the results are closely related, but significantly different.In Table3we have included numerical values emphasizing the situation described by Remark 2.2. We consider the Bernoulli model of parameter p = 0.0001 over the region R of size 185 × 185 × 185 and scan it with a cubic window of length 10. The second and forth columns gives the values corresponding to the bounds described in Eq.(2.20), while in the third column we presented the simulated values for P (S 10,10,10 (185, 185, 185) ≤ n).

Table 2 .

 2 Approximation for P(S ≤ n) over the region R with windows of the same volume by different sizes: T1 = T2 = T3 = 60, p = 0.0025, IT ER = 10 5

	n P(S ≤ n)	Our	Eapp	E sim	Total
		Approximation Eq.(2.29) Eq.(2.36)	Error
		m 1 = m 2 = m 3 = 4
	5 0.961691	0.963506	0.000038	0.003622 0.003660
	6 0.999006	0.999023	0.000000	0.000071 0.000071
	7 0.999980	0.999980	0.000000	0.000001 0.000001
	8 0.999999	0.999999	0.000000 2 × 10 -9 2 × 10 -9
		m 1 = 8, m 2 = 4, m 3 = 2
	5 0.969189	0.969110	0.000007	0.003387 0.003395
	6 0.999297	0.999228	0.000000	0.000071 0.000071
	7 0.999984	0.999984	0.000000	0.000001 0.000001
	8 0.999999	0.999999	0.000000 2 × 10 -9 2 × 10 -9

Table 3 .

 3 Approximation for P(S ≤ n) based on Eq.(2.20): m1

Table 4 .

 4 Approximation for P(S ≤ n) in Binomial and Poisson cases: m1 = m2 = m3 = 4, T1 = T2 = T3 = 84, IT ER = 10 5

	n P(S ≤ n)	Our	Eapp	E sim	Total
		Approximation	Eq.(2.29)	Eq.(2.36)	Error
		Binomial :	m = 10, p = 0.0025
	10 0.726386	0.723224	0.007763	0.032197 0.039960
	11 0.954605	0.955417	0.000123	0.003079 0.003202
	12 0.993938	0.993906	0.000001	0.000331 0.000333
	13 0.999289	0.999284	0.000000	0.000033 0.000033
	14 0.999923	0.999921	0.000000	0.000003 0.000003
	15 0.999992	0.999992	0.000000	3 × 10 -7 3 × 10 -7
		P oisson :	λ = 0.025
	10 0.713184	0.708481	0.009211	0.035294 0.044506
	11 0.950947	0.950197	0.000143	0.003345 0.003488
	12 0.993624	0.993452	0.000002	0.000365 0.000367
	13 0.999218	0.999210	0.000000	0.000038 0.000038
	14 0.999912	0.999911	0.000000	0.000003 0.000003
	15 0.999990	0.999990	0.000000	3 × 10 -7 3 × 10 -7