Umeå University's logo

In this thesis, we study numerical analysis for random processes and fields. We investigate the behavior of the approximation accuracy for specific linear methods based on a finite number of observations. Furthermore, we propose techniques for optimizing performance of the methods for particular classes of random functions. The thesis consists of an introductory survey of the subject and related theory and four papers (A-D).

In paper A, we study a Hermite spline approximation of quadratic mean continuous and differentiable random processes with an isolated point singularity. We consider a piecewise polynomial approximation combining two different Hermite interpolation splines for the interval adjacent to the singularity point and for the remaining part. For locally stationary random processes, sequences of sampling designs eliminating asymptotically the effect of the singularity are constructed.

In Paper B, we focus on approximation of quadratic mean continuous real-valued random fields by a multivariate piecewise linear interpolator based on a finite number of observations placed on a hyperrectangular grid. We extend the concept of local stationarity to random fields and for the fields from this class, we provide an exact asymptotics for the approximation accuracy. Some asymptotic optimization results are also provided.

In Paper C, we investigate numerical approximation of integrals (quadrature) of random functions over the unit hypercube. We study the asymptotics of a stratified Monte Carlo quadrature based on a finite number of randomly chosen observations in strata generated by a hyperrectangular grid. For the locally stationary random fields (introduced in Paper B), we derive exact asymptotic results together with some optimization methods. Moreover, for a certain class of random functions with an isolated singularity, we construct a sequence of designs eliminating the effect of the singularity.

In Paper D, we consider a Monte Carlo pricing method for arithmetic Asian options. An estimator is constructed using a piecewise constant approximation of an underlying asset price process. For a wide class of Lévy market models, we provide upper bounds for the discretization error and the variance of the estimator. We construct an algorithm for accurate simulations with controlled discretization and Monte Carlo errors, andobtain the estimates of the option price with a predetermined accuracy at a given confidence level. Additionally, for the Black-Scholes model, we optimize the performance of the estimator by using a suitable variance reduction technique.

We present a procedure for performing an ANOVA test on functional data, including pairwise group comparisons. in a Scheff´e-like perspective. The test is based on the Interval Testing Procedure, and it selects intervals where the groups significantly differ. The procedure is applied on the 3D kinematic motion of the knee joint collected during a functional task (one leg hop) performed by three groups of individuals.

Motivated by the analysis of the dependence of knee movement patterns during functional tasks on subject-specific covariates, we introduce a distribution-free procedure for testing a functional-on-scalar linear model with fixed effects. The procedure does not only test the global hypothesis on the entire domain but also selects the intervals where statistically significant effects are detected. We prove that the proposed tests are provided with an asymptotic control of the intervalwise error rate, that is, the probability of falsely rejecting any interval of true null hypotheses. The procedure is applied to one-leg hop data from a study on anterior cruciate ligament injury. We compare knee kinematics of three groups of individuals (two injured groups with different treatments and one group of healthy controls), taking individual-specific covariates into account.

Functional data are smooth, often continuous, random curves, which can be seen as an extreme case of multivariate data with infinite dimensionality. Just as component-wise inference for multivariate data naturally performs feature selection, subset-wise inference for functional data performs domain selection. In this paper, we present a unified testing framework for domain selection on populations of functional data. In detail, p-values of hypothesis tests performed on point-wise evaluations of functional data are suitably adjusted for providing a control of the family-wise error rate (FWER) over a family of subsets of the domain. We show that several state-of-the-art domain selection methods fit within this framework and differ from each other by the choice of the family over which the control of the FWER is provided. In the existing literature, these families are always defined a priori. In this work, we also propose a novel approach, coined threshold-wise testing, in which the family of subsets is instead built in a data-driven fashion. The method seamlessly generalizes to multidimensional domains in contrast to methods based on a-priori defined families. We provide theoretical results with respect to consistency and control of the FWER for the methods within the unified framework. We illustrate the performance of the methods within the unified framework on simulated and real data examples, and compare their performance with other existing methods.

We propose a new nonparametric clustering method for dependent functional data, the double clustering bagging Voronoi method. It consists of two levels of clustering. Given a spatial lattice of points, a function is observed at each grid point. In the first‐level clustering, features of the functional data are clustered. The second‐level clustering takes dependence into account, by grouping local representatives, built from the resulting first‐level clusters, using a bagging Voronoi strategy. Depending on the distance measure used, features of the functions may be included in the second‐step clustering, making the method flexible and general. Combined with the clustering method, a multiresolution approach is proposed that searches for stable clusters at different spatial scales, aiming to capture latent structures. This provides a powerful and computationally efficient tool to cluster dependent functional data at different spatial scales, here illustrated by a simulation study. The introduced methodology is applied to varved lake sediment data, aiming to reconstruct winter climate regimes in northern Sweden at different time resolutions over the past 6,000 years.

We consider a multivariate piecewise linear interpolation of a continuous random field on a-dimensional cube. The approximation performance is measured by the integrated mean square error. Multivariate piecewise linear interpolator is defined by N field observations on a locations grid (or design). We investigate the class of locally stationary random fields whose local behavior is like a fractional Brownian field in mean square sense and find the asymptotic approximation accuracy for a sequence of designs for large N. Moreover, for certain classes of continuous and continuously differentiable fields we provide the upper bound for the approximation accuracy in the uniform mean square norm.

We consider a Monte Carlo method to price a continuous arithmetic Asian option with a given precision. Piecewise constant approximation and plain simulation are used for a wide class of models based on L\'{e}vy processes. We give bounds of the possible discretization and simulation errors. The sufficient numbers of discretization points and simulations to obtain requested accuracy are derived. To demonstrate the general approach, the Black-Scholes model is studied in more detail. We undertake the case of continuous averaging and starting time zero,  but the obtained results can be applied to the discrete case  and generalized for any time before an execution date. Some numerical experiments and comparison to the PDE based method are also presented.

We consider a piecewise-multilinear interpolation of a continuous random field on a d-dimensional cube. The approximation performance is measured using the integrated mean square error. Piecewise-multilinear interpolator is defined by N-field observations on a locations grid (or design). We investigate the class of locally stationary random fields whose local behavior is like a fractional Brownian field, in the mean square sense, and find the asymptotic approximation accuracy for a sequence of designs for large N. Moreover, for certain classes of continuous and continuously differentiable fields, we provide the upper bound for the approximation accuracy in the uniform mean square norm.

Let a continuous random process X defined on [0,1] be (m+β)-smooth, 0m, 0<β1, in quadratic mean for all t>0 and have an isolated singularity point at t=0. In addition, let X be locally like a m-fold integrated β-fractional Brownian motion for all nonsingular points. We consider approximation of X by piecewise Hermite interpolation splines with n free knots (i.e., a sampling design, a mesh). The approximation performance is measured by mean errors (e.g., integrated or maximal quadratic mean errors). We construct a sequence of sampling designs with asymptotic approximation rate n^(m+β) for the whole interval.

We consider the problem of numerical approximation of integrals of random fields over a unit hypercube. We use a stratified Monte Carlo quadrature and measure the approximation performance by the mean squared error. The quadrature is defined by a finite number of stratified randomly chosen observations with the partition generated by a rectangular grid (or design). We study the class of locally stationary random fields whose local behavior is like a fractional Brownian field in the mean square sense and find the asymptotic approximation accuracy for a sequence of designs for large number of the observations. For the H¨older class of random functions, we provide an upper bound for the approximation error. Additionally, for a certain class of isotropic random functions with an isolated singularity at the origin, we construct a sequence of designs eliminating the effect of the singularity point.

Nonparametric bagging clustering methods are studied and compared to identify latent structures from a sequence of dependent categorical data observed along a one-dimensional (discrete) time domain. The frequency of the observed categories is assumed to be generated by a (slowly varying) latent signal, according to latent state-specific probability distributions. The bagging clustering methods use random tessellations (partitions) of the time domain and clustering of the category frequencies of the observed data in the tessellation cells to recover the latent signal, within a bagging framework. New and existing ways of generating the tessellations and clustering are discussed and combined into different bagging clustering methods. Edge tessellations and adaptive tessellations are the new proposed ways of forming partitions. Composite methods are also introduced, that are using (automated) decision rules based on entropy measures to choose among the proposed bagging clustering methods. The performance of all the methods is compared in a simulation study. From the simulation study it can be concluded that local and global entropy measures are powerful tools in improving the recovery of the latent signal, both via the adaptive tessellation strategies (local entropy) and in designing composite methods (global entropy). The composite methods are robust and overall improve performance, in particular the composite method using adaptive (edge) tessellations.

In this paper we introduce a novel functional clustering method, the Bagging Voronoi K-Medoid Aligment (BVKMA) algorithm, which simultaneously clusters and aligns spatially dependent curves. It is a nonparametric statistical method that does not rely on distributional or dependency structure assumptions. The method is motivated by and applied to varved (annually laminated) sediment data from lake Kassjön in northern Sweden, aiming to infer on past environmental and climate changes. The resulting clusters and their time dynamics show great potential for seasonal climate interpretation, in particular for winter climate changes.

Fiber (tracheid) length is an important trait targeted for genetic and silvicultural improvement. Such studies require large-scale non-destructive sampling, and accurate length determination. The standard procedure for non-destructive sampling is to collect increment cores, singularize their cells by maceration, measure them with optical analyzer and apply various corrections to suppress influence of non-fiber particles and cut fibers, as fibers are cut by the corer. The recently developed expectation-maximization method (EM) not only addresses the problem of non-fibers and cut fibers, but also corrects for the sampling bias. Here, the performance of the EM method has been evaluated by comparing it with length-weighing and squared length-weighing, both implemented in fiber analyzers, and with microscopy data for intact fibers, corrected for sampling bias, as the reference. This was done for 12-mm increment cores from 16 Norway spruce (Picea abies (L.) Karst) trees on fibers from rings 8-11 (counted from pith), representing juvenile wood of interest in breeding programs. The EM-estimates provided mean-fiber-lengths with bias of only +2.7% and low scatter. Length-weighing and length2-weighing gave biases of -7.3% and +9.3%, respectively, and larger scatter. The suggested EM approach constitutes a more accurate non-destructive method for fiber length (FL) determination, expected to be applicable also to other conifers.

BACKGROUND: Despite interventions, anterior cruciate ligament ruptures can cause long-term deficits. To assist in identifying and treating deficiencies, 3D-motion analysis is used for objectivizing data. Conventional statistics are commonly employed to analyze kinematics, reducing continuous data series to discrete variables. Conversely, functional data analysis considers the entire data series.

METHODS: Here, we employ functional data analysis to examine and compare the entire time-domain of knee-kinematic curves from one-leg hops between and within three groups. All subjects (n=95) were part of a long-term follow-up study involving anterior cruciate ligament ruptures treated ~20years ago conservatively with physiotherapy only or with reconstructive surgery and physiotherapy, and matched knee-healthy controls.

FINDINGS: Between-group differences (injured leg, treated groups; non-dominant leg, controls) were identified during the take-off and landing phases, and in the sagittal (flexion/extension) rather than coronal (abduction/adduction) and transverse (internal/external) planes. Overall, surgical and control groups demonstrated comparable knee-kinematic curves. However, compared to controls, the physiotherapy-only group exhibited less flexion during the take-off (0-55% of the normalized phase) and landing (44-73%) phase. Between-leg differences were absent in controls and the surgically treated group, but observed during the flight (4-22%, injured leg>flexion) and the landing (57-85%, injured leg<internal rotation) phases in the physiotherapy-only group.

INTERPRETATION: Functional data analysis identified specific functional knee-joint deviations from controls persisting 20years post anterior cruciate ligament rupture, especially when treated conservatively. This approach is suggested as a means for comprehensively analyzing complex movements, adding to previous analyses.

The recent sports science literature conveys a growing interest in robust statistical methods to analyze smooth, regularly-sampled functional data. This paper focuses on the inferential problem of identifying the parts of a functional domain where two population means differ. We considered four approaches recently used in sports science: interval-wise testing (IWT), statistical parametric mapping (SPM), statistical nonparametric mapping (SnPM) and the Benjamini-Hochberg (BH) procedure for false discovery control. We applied these procedures to both six representative sports science datasets, and also to systematically varied simulated datasets which replicated ten signal- and/or noise-relevant parameters that were identified in the experimental datasets. We observed generally higher IWT and BH sensitivity for five of the six experimental datasets. BH was the most sensitive procedure in simulation, but also had relatively high false positive rates (generally > 0.1) which increased sharply (> 0.3) in certain extreme simulation scenarios including highly rough data. SPM and SnPM were more sensitive than IWT in simulation except for (1) high roughness, (2) high nonstationarity, and (3) highly nonuniform smoothness. These results suggest that the optimum procedure is both signal and noise-dependent. We conclude that: (1) BH is most sensitive but also susceptible to high false positive rates, (2) IWT, SPM and SnPM appear to have relatively inconsequential differences in terms of domain identification sensitivity, except in cases of extreme signal/noise characteristics, where IWT appears to be superior at identifying a greater portion of the true signal.

Allocation of carbohydrates between competing organs is fundamental to plant development, growth and productivity. Carbohydrates are synthesized in mature leaves and distributed via the phloem vasculature to developing buds where they are consumed to produce new biomass. The distribution and mass-allocation processes within the plant remain poorly understood and may involve complex feedbacks between different plant functions, with implications for the emergent structure of the plant. Here, we investigate how the order in which dormant buds are flushed affects the development of tree size and reproductive output during the first 20 years of growth in full light and shaded canopy environments. We report the following findings: (i) Bud-flushing strategies strongly affect the temporal dynamics of height, mass and the size of reproduction pool, as well as the resulting architectures. (ii) Bud-flushing strategies affect tree growth by altering the rate of growth and final size of trees. (iii) No single bud-flushing strategy performs best when both the size and allocation for reproduction of the resulting trees are compared. However, we observe that the strategy that optimizes the net carbon gain for the entire tree architecture always results in a high reproduction output. (iv) Branch turnover and meristem regeneration enhance the performance of certain strategies with respect to the measured quantities. These results highlight the importance of employing generic models of architecture (i.e., non-species-specific) to identify general mechanisms of carbon allocation and the spatial distribution of newly formed biomass in growing trees.

One of the main concerns of silviculture and forest management focuses on finding fast, cost-efficient and non-destructive ways of measuring wood properties in standing trees. This paper presents an R package \verb+fiberLD+ that provides functions for estimating tree fiber length distributions in the standing tree based on increment core samples. The methods rely on increment core data measured by means of an optical fiber analyzer (OFA) or measured by microscopy. Increment core data analyzed by OFAs consist of the cell lengths of both cut and uncut fibers (tracheids) and fines (such as ray parenchyma cells) without being able to identify which cells are cut or if they are fines or fibers. The microscopy measured data consist of the observed lengths of the uncut fibers in the increment core. A censored version of a mixture of the fine and fiber length distributions is proposed to fit the OFA data, under distributional assumptions. Two choices for the assumptions of the underlying density functions of the true fiber (fine) lengths of those fibers (fines) that at least partially appear in the increment core are considered, such as the generalized gamma and the log normal densities. Maximum likelihood estimation is used for estimating the model parameters for both the OFA analyzed data and the microscopy measured data.

1. While theoretical models predict reproductive allocation (RA) should approach 100% of available energy as a plant ages, available empirical data suggest much lower RA values in perennial plants. In this study, we have two aims. First, we assess whether the discrepancy between theory and data arises from methodological differences in how growth and RA are calculated. Specifically, we hypothesize RA in older plants is large when compared to growth in leaf area, that is, after excluding turnover of stem and leaf tissues. Second, we hypothesize that species with cheap tissues or those that are shorter reach RA = 0.5 at a younger age.

2. We measured investment in leaf, stem and reproduction on individuals from 14 co-occurring woody perennial iteroparous species. A fire chronosequence allowed us to use a space-for-time substitution to estimate RA schedules for each species, simultaneously measuring reproductive and vegetative production on individuals differing in age.

3. For most (11 of 14) species, we found RA eventually reached 100% of available energy, with another two species reaching at least 80%. Increases in RA were associated with a decline in growth of leaf area. Comparing species, we found that species with cheap leaves reached RA = 0.5 sooner (they could be called fast-living), whereas delayed maturation and slower increases in RA were associated with greater maximum height.

4. Synthesis. Explicitly accounting for the cost of leaf replacement leads to the high estimates of reproductive allocation in perennial plants predicted by theoretical models, limiting or even halting leaf area expansion. For some species, so much energy is allocated to reproduction that leaf area declines year-upon-year for multiple growing seasons preceding death. Connecting lifetime reproductive allocation schedules with leaf area expansion, leaf life span, and plant maximum height demonstrates how reproductive allocation schedules synthesize a plant's life-history strategy, making them a valuable tool for connecting plant traits and demography.