Dynamic Duration of Load Models by Yongliang Zhai B.S., Zhejiang University, 2009 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF Master of Science in THE FACULTY OF GRADUATE STUDIES (Statistics) The University Of British Columbia (Vancouver) August 2011 c© Yongliang Zhai, 2011 Abstract The duration of load effect is a distinctive and important characteristic of wood strength. It refers to the fact that wood products can usually sustain a high load for a short time but the products may deteriorate and break in the long run. Modelling the duration of load effect and testing wood for specific properties of this effect are important in formulating wood construction standards. Damage accumulation models have been proposed by authors to model the duration of load effects. The models assume that damage is accumulated over time according to the load history, and once the accumulated damage reaches a threshold value, the board will break. Different authors have designed different experiments and proposed different methods for estimating the model parameters. In this work, we consider several damage accumulation models, with a focus on the U.S. model. We investigate the effects of the distributional assumptions for the models, and propose several methods to estimate parameters in the models. Our proposed methods are evaluated via simulation studies. Two real datasets are present for illustration. ii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvi 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2 Stochastic Duration of Load Models . . . . . . . . . . . . . . . . . . 7 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2 Madison Curve . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.3 The U.S. Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.3.1 Model Description . . . . . . . . . . . . . . . . . . . . . 11 2.3.2 Common Steps in Gerhards-Link and Foschi-Yao Analyses 12 2.3.3 Solutions for Ts and Tc from Gerhards and Link’s (1987) Approach . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.3.4 Solutions for Ts and Tc from Foschi and Yao’s (1986) Ap- proach . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.4 The Canadian Model and its Predecessors . . . . . . . . . . . . . 15 iii 2.4.1 Models Description . . . . . . . . . . . . . . . . . . . . . 15 2.4.2 Solutions for Ts and Tc . . . . . . . . . . . . . . . . . . . 17 2.4.3 Parameter Estimation . . . . . . . . . . . . . . . . . . . . 19 2.5 The Problem of Scale . . . . . . . . . . . . . . . . . . . . . . . . 21 3 Some Results Based on Simulated Data . . . . . . . . . . . . . . . . 25 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.2 Generating Data Based on the U.S. Model . . . . . . . . . . . . . 26 3.2.1 A Brief Summary of the U.S. Model . . . . . . . . . . . . 26 3.2.2 Simulation Setups . . . . . . . . . . . . . . . . . . . . . . 28 3.2.3 The Effect of the Type of Distribution of the Random Effects 28 3.2.4 The Effect of the Mean and the Standard Deviation of the Random Effects . . . . . . . . . . . . . . . . . . . . . . . 39 3.3 Generating Data Based on the Canadian model . . . . . . . . . . . 43 3.3.1 A Brief Summary of the Canadian Model . . . . . . . . . 43 3.3.2 Simulation Setups . . . . . . . . . . . . . . . . . . . . . . 45 3.4 Comparisons of the Generated Data with Foschi and Yao’s Data . 46 3.5 Fitting the Generated Breaking Times to Distributions . . . . . . . 48 4 Estimation Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.2 Models and Notation . . . . . . . . . . . . . . . . . . . . . . . . 54 4.3 Likelihood Method . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.3.1 The Likelihood Functions . . . . . . . . . . . . . . . . . . 56 4.3.2 Optimization Method . . . . . . . . . . . . . . . . . . . . 58 4.4 Approximations . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.4.1 Approximation Step 1 . . . . . . . . . . . . . . . . . . . . 62 4.4.2 Approximation Step 2 . . . . . . . . . . . . . . . . . . . . 62 4.4.3 Approximation Step 3: Linearization . . . . . . . . . . . . 64 4.4.4 Approximate Likelihood Functions . . . . . . . . . . . . . 69 4.5 Quantile Method . . . . . . . . . . . . . . . . . . . . . . . . . . 69 iv 4.5.1 Estimates of the Means . . . . . . . . . . . . . . . . . . . 71 4.5.2 Estimates of the Variances . . . . . . . . . . . . . . . . . 76 4.5.3 Approximation Errors . . . . . . . . . . . . . . . . . . . . 77 4.6 A Method that Combines the Maximum Likelihood and the Quan- tile Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5 Simulation Studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.2 General Setup for the Simulation Studies . . . . . . . . . . . . . . 80 5.3 Simulation Studies for the Approximate Maximum Likelihood Es- timates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.3.1 Simulation Setups . . . . . . . . . . . . . . . . . . . . . . 80 5.3.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 5.3.3 The Shapes of the Likelihood Functions . . . . . . . . . . 82 5.4 Simulation Studies for the Quantile Method . . . . . . . . . . . . 84 5.4.1 Simulation Setups . . . . . . . . . . . . . . . . . . . . . . 84 5.4.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.5 Simulation Studies for the Combined Method . . . . . . . . . . . 89 5.5.1 Simulation Setups . . . . . . . . . . . . . . . . . . . . . . 89 5.5.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 5.6 Summary of the Simulation Results . . . . . . . . . . . . . . . . . 94 6 Experiments and Data Analysis . . . . . . . . . . . . . . . . . . . . . 96 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 6.2 The FPInnovations Experiment and Data Analysis . . . . . . . . . 97 6.2.1 The FPInnovations Experiment . . . . . . . . . . . . . . . 97 6.2.2 Data Analysis for Experiments at FPInnovations . . . . . . 99 6.3 Foschi and Barrett’s Experiments and Data Analysis . . . . . . . . 101 6.3.1 Foschi and Barrett’s Experiments . . . . . . . . . . . . . . 102 6.3.2 Data Analysis for Foschi and Barrett’s Experiments . . . . 104 v 7 Conclusions and Future Work . . . . . . . . . . . . . . . . . . . . . 107 7.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 7.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 vi List of Tables Table 2.1 Commonly used notation for damage accumulation models. . . 8 Table 3.1 Summary of the setups of the simulation runs for studying the effect of the type of distribution of the random effects a and b when µa = 42 and µb = 50, and p = 0.2 for the constant load test. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 Table 3.2 Summary of the setups of the simulation runs for studying the effect of the mean and the standard deviation of the random effects a and b when a and b are both Normal, and p = 0.2 in the constant load test. . . . . . . . . . . . . . . . . . . . . . . 40 Table 3.3 Summary of the setups of the simulation runs for studying the distribution of the generated breaking times when µa = 42 and µb = 50, and p = 0.2 in the constant load test. . . . . . . . . . 49 Table 5.1 Summary of the simulation results in Chapter 5 when ns = nc = 200. In all simulation runs, (µa,µb) = (42,50) and in the con- stant load test p = 0.2. The values of (σa,σb) are shown in the table. The success rates denote the proportion of the converged simulation runs for the approximate maximum likelihood esti- mates and the combined method, and the rates of the simulation runs which produce positive estimates for both σ2a and σ2b for the quantile method. . . . . . . . . . . . . . . . . . . . . . . . 94 vii List of Figures Figure 1.1 An illustration of the ramp load test and the constant load test. In the ramp load test, the load increases linearly over time until the piece fails. In the constant load test, the load first increases linearly until it reaches a pre-determined load τa at time T0 and then stays constant thereafter until the piece fails or the experiment ends. . . . . . . . . . . . . . . . . . . . . 3 Figure 1.2 A picture of the ramp load test conducted at FPInnovations. The board is placed on the machine horizontally. It is sup- ported by two ends and the load is applied at the centre. The board will bend and finally break when the load increases. . . 4 Figure 1.3 A picture of the constant load test conducted at FPInnovations. The board is placed on the machine vertically. The board is supported by two ends and the load is applied to the centre through a pulley. . . . . . . . . . . . . . . . . . . . . . . . . 5 Figure 2.1 Data from constant load tests by Wood (1951), reproduced from Gerhards (1977). The x-axis is the breaking time T in hours and the y-axis is the applied load ratio σ in percentage. 9 Figure 2.2 Madison Curve with ramp loading test trends and constant loading test trends by Wood (1951), reproduced from Ger- hards (1977). The x-axis is the breaking time T in seconds and the y-axis is the applied load ratio σ in percentage. . . . . 10 viii Figure 3.1 Comparisons of the generated Ts in the U.S. model based on the three different distributional assumptions of a and b when µa = 42,µb = 50,σa = 1 and σb = 1. Upper left: the empir- ical cumulative distributions of log10(Ts) from the three sce- narios. Upper right: the quantile-quantile plot of the generated log10(Ts) from the assumption that a and b are both Normal and the assumption that a and b are both log-normal. Lower left: the quantile-quantile plot of the generated log10(Ts) from the assumption that a and b are both Normal and the assump- tion that a is Normal and b is log-normal. Lower right: the quantile-quantile plot of the generated log10(Ts) from the as- sumption that a and b are both log-normal and the assumption that a is Normal and b is log-normal. . . . . . . . . . . . . . 31 Figure 3.2 Comparisons of the generated Tc in the U.S. model based on the three different distributional assumptions of a and b when µa = 42,µb = 50,σa = 1, σb = 1 and p = 0.2. The panels are in the same arrangement as in Figure 3.1. . . . . . . . . . . . 32 Figure 3.3 Comparisons of the generated Ts in the U.S. model based on the three different distributional assumptions of a and b when µa = 42,µb = 50,σa = 10 and σb = 10. The panels are in the same arrangement as in Figure 3.1. . . . . . . . . . . . . . . 33 Figure 3.4 Comparisons of the generated Tc in the U.S. model based on the three different distributional assumptions of a and b when µa = 42,µb = 50,σa = 10, σb = 10 and p = 0.2. The panels are in the same arrangement as in Figure 3.1. . . . . . . . . . 34 Figure 3.5 Comparisons of the generated Ts in the U.S. model based on the three different distributional assumptions of a and b when µa = 42,µb = 50,σa = 10 and σb = 1. The panels are in the same arrangement as in Figure 3.1. . . . . . . . . . . . . . . 35 ix Figure 3.6 Comparisons of the generated Tc in the U.S. model based on the three different distributional assumptions of a and b when µa = 42,µb = 50,σa = 10, σb = 1 and p = 0.2. The panels are in the same arrangement as in Figure 3.1. . . . . . . . . . 36 Figure 3.7 Comparisons of the generated Ts in the U.S. model based on the three different distributional assumptions of a and b when µa = 42,µb = 50,σa = 1 and σb = 10. The panels are in the same arrangement as in Figure 3.1. . . . . . . . . . . . . . . 37 Figure 3.8 Comparisons of the generated Tc in the U.S. model based on the three different distributional assumptions of a and b when µa = 43,µb = 50,σa = 1, σb = 10 and p = 0.2. The panels are in the same arrangement as in Figure 3.1. . . . . . . . . . 38 Figure 3.9 Comparisons of the empirical cumulative distributions of the generated Ts and Tc from the U.S. model based the assumption that a and b are both Normal, σa = 1 and σb = 1. One mean is fixed at the value stated in the plot title. . . . . . . . . . . . 41 Figure 3.10 Comparisons of the empirical cumulative distributions of the generated Ts and Tc from the U.S. model based the assumption that a and b are both Normal, µa = 42 and µb = 50. One standard variation is fixed at the value stated in the plot title. . 42 Figure 3.11 Comparisons of the empirical cumulative distributions of break- ing times Tc’s from Foschi and Yao’s experiments and the gen- erated Tc based on the U.S. model (left panel), and the Cana- dian model (right panel). . . . . . . . . . . . . . . . . . . . . 47 Figure 3.12 The quantile-quantile plots of the logarithm of the generated breaking times Ts based on the U.S. model when (µa,µb,σa,σb) = (42,50, 0.4,0.4) and the logarithm of simulated data from the fitted distributions. The type of fitted distribution is shown in the title of the plot. . . . . . . . . . . . . . . . . . . . . . 50 x Figure 3.13 The quantile-quantile plots of the logarithm of the generated breaking times Ts based on the U.S. model when (µa,µb,σa,σb) = (42,50, 10,10) and the logarithm of simulated data from the fitted distributions. The type of fitted distribution is shown in the title of the plot. . . . . . . . . . . . . . . . . . . . . . . . 51 Figure 3.14 The quantile-quantile plots of the logarithm of the generated breaking times Tc based on the U.S. model when (µa,µb,σa,σb) = (42,50, 0.4,0.4) and p = 0.2, and the logarithm of simu- lated data from the fitted distributions. The fitted distribution is shown in the title of the plot. . . . . . . . . . . . . . . . . 51 Figure 3.15 The quantile-quantile plots of the logarithm of the generated breaking times Tc based on the U.S. model when (µa,µb,σa,σb) = (42,50, 10,10) and p = 0.2, and the logarithm of simulated data from the fitted distributions. The type of fitted distribu- tion is shown in the title of the plot. . . . . . . . . . . . . . . 52 Figure 3.16 The quantile-quantile plots of the logarithm of the generated breaking times Tc based on the Canadian model and the loga- rithm of simulated data from the fitted distributions. The type of fitted distribution is shown in the title of the plot. . . . . . 52 Figure 4.1 Plots of the approximation error versus the logarithm of the breaking time Tc when a and b are both Normal, µa = 42 and µb = 50, and p= 0.2 in the constant load test. The vertical line is log(T0). The values for σa and σb are shown in the title. Tc’s are the breaking times generated from the U.S. model without approximations and Tc,approx2’s are the approximated breaking times generated from the U.S. model with the approximation step 1 and 2. The approximation error is the difference of log(Tc,approx2) and log(Tc). . . . . . . . . . . . . . . . . . . . 63 xi Figure 4.2 Plots of the approximation error versus the logarithm of the breaking time Ts when a and b are both Normal, µa = 42 and µb = 50. The values for σa and σb are shown in the title. Ts’s are the breaking times generated from the U.S. model without approximations and Ts,approx3’s are the approximated breaking times generated from the U.S. model with the approximation steps 1, 2 and 3. The approximation error is the difference of log(Ts,approx3) and log(Ts). . . . . . . . . . . . . . . . . . . . 67 Figure 4.3 Plots of the approximation error versus the logarithm of the breaking time Tc when a and b are both Normal, µa = 42 and µb = 50, and p = 0.2 in the constant load test. The values for σa and σb are shown in the title. Tc’s are the breaking times generated from the U.S. model without approximations and Tc,approx3’s are the approximated breaking times generated from the U.S. model with the approximation steps 1, 2 and 3.The vertical line denotes the logarithm of the median of the breaking times Tc’s. The approximation error is the difference of log(Tc,approx3) and log(Tc). . . . . . . . . . . . . . . . . . 68 Figure 4.4 Plots for illustrating the quantile method when a and b are both Normal, µa = 42 and µb = 50, and p= 0.2 in the constant load test. The values for σa and σb are shown in the title. We paint all the simulated points grey, and then paint those points which satisfy Tc > mY black. These plots support our claim in equation (4.31). . . . . . . . . . . . . . . . . . . . . . . . . 75 Figure 5.1 Bar plots of the optimization codes for calculating the approx- imate maximum likelihoodestimates in nsim = 100 simulation runs, when θ = (µa,µb,σa,σb) = (42,50,0.4,0.4), ns = nc = 200 and p = 0.2 in the constant load test. . . . . . . . . . . . 81 xii Figure 5.2 Perspective plots of the log-likelihoods when σa = 0.4 and σb = 0.4 are fixed, and µa varies from 40 to 45 and µb varies from 47 to 52. . . . . . . . . . . . . . . . . . . . . . . . . . 83 Figure 5.3 Box lots of the estimates using the quantile method, when θ = (42,50,0.4,0.4) and p = 0.2. The sample size ns and nc used in the simulation runs are shown in plot titles. The grey line indicates the true value of the parameter. . . . . . . . . . . . 86 Figure 5.4 Box lots of the estimates using the quantile method, when θ = (µa,µb,σa,σb) = (42,50,5,5) and p = 0.2. The sample size ns = nc = 200. The grey line indicates the true value of the parameter. . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 Figure 5.5 Box lots of the estimates using the quantile method, when θ = (µa,µb,σa,σb) = (42,50,5,5) and p = 0.2. The sample size ns = nc = 200,000. No grey line is shown in the plots because the estimates are too far from the true values. . . . . . . . . . 88 Figure 5.6 Bar plots of the output of the optimization results in calcu- lating the approximate maximum likelihoodestimates for σa and σb in the combined method when θ = (µa,µb,σa,σb) = (42,50,0.4,0.4) and p = 0.2 in the constant load test. Fixed starting point means the starting point in the optimization step of the combined method is set to be the true value (σa,σb) = (0.4,0.4). Random starting point means the starting point in the optimization step of the combined method is set to be dif- ferent random numbers uniformly generated from the square 0 < σa < 1 and 0 < σb < 1 for different simulation runs. . . . 90 xiii Figure 5.7 Plots of the estimates for σa and σb when θ =(42,50,0.4,0.4) and p = 0.2 in the constant load test. Fixed starting point means the starting point in the optimization step of the com- bined method is set to be the true value (σa,σb) = (0.4,0.4). Random starting point means the starting point in the opti- mization step of the combined method is set to be different random numbers in the square 0 < σa < 1 and 0 < σb < 1 for different simulation runs. The grey line indicates the true value of the parameter. . . . . . . . . . . . . . . . . . . . . . 92 Figure 5.8 Plots of the estimates for σa and σb when θ = (42,50,5,5) and p = 0.2 in the constant load test. Fixed starting point means the starting point in the optimization step of the com- bined method is set to be the true value (σa,σb) = (5,5). Ran- dom starting point means the starting point in the optimiza- tion step of the combined method is set to be different random numbers in the square 0 < σa < 10 and 0 < σb < 10 for dif- ferent simulation runs. The grey line indicates the true value of the parameter. . . . . . . . . . . . . . . . . . . . . . . . . 93 Figure 6.1 The deflection history and the load history of one specimen in the ramp load test at FPInnovations. . . . . . . . . . . . . . . 99 Figure 6.2 Empirical cumulative distributions of log(Ts) from the ramp load test results at FPInnovations. In the plot, DF-High, DF- Low, LVL-High and LVL-Low denote Douglas-fir lumber with high moisture content, Douglas-fir lumber with low moisture content, laminated-veneer lumber with hight moisture content and laminated-veneer lumber with low moisture content re- spectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 Figure 6.3 The quantile-quantile plots of the breaking times Ts’s (in log- scale) and the fitted log-normal distributions for Ts’s in each group. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 xiv Figure 6.4 Empirical cumulative distribution functions of log(Tc) from the constant load test results in the FPInnovations. In the plot, DF-High and DF-Low denote Douglas-fir lumber with high moisture content and Douglas-fir lumber with low moisture content respectively. There are 20 boards in each group. . . . 103 Figure 6.5 Box plots of the 89 estimates of σa and σb. . . . . . . . . . . 105 xv Acknowledgements I would like to express my deepest gratitude to my supervisors, Dr. Nancy Heck- man and Dr. Lang Wu, for their support, guidance and patience. From them, I learned how to conduct a research project, and more importantly, how to become a successful statistician. I am also touched by the time they spent with me on revising this thesis, which helped me improving it greatly from the first version. I would like to thank Dr. James Zidek and the Forest Products Stochastic Modelling Group, for providing me the financial support and helpful comments on this project. I would like to thank our collaborators, Dr. Ciprian Pirvu and Dr. Conroy Lum, at FPInnovations. I also thank Dr. Ricardo Foschi for releasing the data from his experiments. I am grateful to the faculty members and the staff of the Statistics Department at UBC, in particular, Dr. Jiahua Chen, Dr. Harry Joe, and Dr. John Petkau, for the knowledge I learned from them and the nice time I had when taking their courses. Last, I leave the warmest part of my heart for my beloved parents, who gave birth to me, enlightened me and educated me with their unconditional support and continuous love. Without them, everything is impossible. xvi Chapter 1 Introduction 1.1 Motivation The duration of load effect is one of the most distinctive and important charac- teristics of wood strength. It refers to the fact that wood products can usually sustain a high load for a short time but the products may deteriorate and break in the long run. The strength of wood products is influenced by both the intensity and the duration of the applied load. Manufactured products, such as steel, do not manifest a duration of load effect. In other words, if a piece of steel can sustain a certain constant load for one hour, it is likely that this piece of steel can sustain the same load for one year. Modelling the duration of load effect and testing wood for specific properties of this effect are important in formulating wood construction standards and in ensuring that those standards are met. Models for the duration of load effect are usually formulated in terms of the accumulation of damage. It is difficult to propose a model for damage accumu- lation based on physical laws since our knowledge of material behavior at the microscopic level where the deterioration happens is generally incomplete. As an alternative, damage accumulation models have been proposed based instead on a combination of our incomplete understanding of the phenomena at the macro- scopic level and experimental data (Yao, 1987). In an accumulation of damage 1 1.1. MOTIVATION model, a board accumulates damage depending on some load τ that may vary over time. The damage accumulated by time t is denoted α(t) with, by conven- tion, α(0) = 0 and α(T ) = 1, where T is the breaking time of the board. Because wood products cannot “cure themselves”, α is a non-decreasing function of t. In practice, the accumulated damage process cannot be observed, but it may be inferred based on the observed breaking times. The advantage of damage accu- mulation models is that they may facilitate the prediction of damage produced by an arbitrary load sequence, which may be useful in reliability tests. All damage accumulation models considered in this thesis are based on the following differential equation: dα(t) dt = f (α(t),τ(t),θ), (1.1) where f is a known function, τ(t) is the known load at time t, and θ is a vector of parameters, usually unknown. The short term strength τs is often included as an argument of f in equation (1.1), in the form of σ(t) = τ(t)/τs. However, authors define τs in different ways. Gerhards and Link (1987) treat τs as a board dependent random parameter with an assumed distribution and do not define τs in terms of any breaking time or load pattern. Foschi and Yao (1986) also treat τs as a board dependent parameter, but define τs as the breaking load τ(Ts) in the ramp load test, when the loading rate k is set so that the mean breaking time is expected to be around one minute. Different definitions of τs lead to different damage accumulation models. We mostly focus on the model of Foschi and Yao (1986) in this thesis. In the literature, the parameter vector θ is often treated as a constant vector, depending on the type of lumber but constant among boards of the same type. However, this implies that all boards of the same type, when subjected to the same load, have the same breaking time T , since, for a fixed θ and load, at most one value of T can satisfy α(T ) = 1. Clearly this is not realistic, as breaking times do vary from board to board. A perhaps more realistic approach is to treat the parameters θ as random effects, which vary from board to board. This approach 2 1.1. MOTIVATION Time Lo ad l ramp load test constant load test (T0,τa) Figure 1.1: An illustration of the ramp load test and the constant load test. In the ramp load test, the load increases linearly over time until the piece fails. In the constant load test, the load first increases linearly until it reaches a pre-determined load τa at time T0 and then stays constant thereafter until the piece fails or the experiment ends. was taken by Foschi and Yao (1982, 1986), as well as Gerhards and Link (1987). There are typically two types of duration of load tests: the ramp load test and the constant load test. Figure 1.1 illustrates the load history in the ramp load test and in the constant load test. In the ramp load test with rate k, the applied load is linear in t, that is τ(t) = kt. In the ramp load test, the breaking time is Ts, and the breaking load is τ(Ts) = kTs. Figure 1.2 depicts the ramp load test conducted at an FPInnovations lab, located on the UBC campus in Vancouver. In the constant load test, the load first increases linearly at constant rate k until a predetermined time T0, similar to the initial period of the ramp load test, and then the load remains constant during the rest of time (see Figure 1.1). That is τ(t) = { kt for 0≤ t ≤ T0, kT0 for t > T0. 3 1.1. MOTIVATION Figure 1.2: A picture of the ramp load test conducted at FPInnovations. The board is placed on the machine horizontally. It is supported by two ends and the load is applied at the centre. The board will bend and finally break when the load increases. The pre-determined load level kT0 is denoted by τa, i.e., τa = kT0. The load level τa is usually set at a certain percentile of the empirical distribution of the short- term strength of the wood products tested during a ramp load test with load equal to kt, the same value of k as in the constant load test. The first part of the constant load test (i.e., τ(t) = kt when 0 ≤ t ≤ T0) is called the ramp loading part of the constant load test and the second part of the constant load test (i.e., τ(t) = τa when t > T0) is called the constant loading part of the constant load test (see Figure 1.1). Figure 1.3 contains a picture of the ramp load test conducted at FPInnovations. Authors have proposed various parameter estimation methods from breaking times Ts’s in the ramp load test and/or breaking times Tc’s in the constant load 4 1.2. OUTLINE Figure 1.3: A picture of the constant load test conducted at FPInnovations. The board is placed on the machine vertically. The board is supported by two ends and the load is applied to the centre through a pulley. test. However, the estimation of parameters has been done in an ad hoc way with no statistical principles, thereby leaving room for possible improvement. These authors have not discussed the distributions of the breaking times. Also, nor has anyone considered random effects in the U.S. model. 1.2 Outline In Chapter 2, we review several damage accumulation models in the literature: the Madison Curve model, the U.S. model, and the Canadian model. We derive the solutions or the approximate solutions they imply for the breaking time Ts in the ramp load test and the breaking time Tc in the constant load test for the U.S. model 5 1.2. OUTLINE and the Canadian model. We also review the parameter estimation methods found in the literature. Finally, we discuss the problem of scale for some models. In Chapter 3, we simulate breaking times from the U.S. model and Canadian model. We investigate the effects of the distributional assumptions on the breaking times Ts and Tc. We compare the results to those in Foschi and Yao’s (1986). We fit the simulated breaking times to some distributions. In Chapter 4, we propose a maximum likelihood based method, a quantile based method, and a method that combines the maximum likelihood and the quan- tile estimates for parameter estimation in the U.S. model. We propose some ap- proximations to the solutions of the breaking times Ts and Tc in the U.S. model, and investigate the accuracy of those approximations. In Chapter 5, we conduct several simulation studies to investigate the perfor- mances of the parameter estimation methods proposed in Chapter 4. In Chapter 6, we describe a trial experiment for the duration of load effects that we conducted at an FPInnovations laboratory. We analyze the empirical dis- tributions of the breaking times from the trial experiment. We also analyze the ex- periment described in Foschi and Barrett (1982). We revise the combined method proposed in Chapter 4 to estimate the parameters using the breaking times Tc’s in the constant load test from Barrett and Foschi’s experiments. In Chapter 7, we summarize the results and discuss possible future work. 6 Chapter 2 Stochastic Duration of Load Models 2.1 Introduction In this chapter, we review some popular damage accumulation models: the Madi- son Curve model, the U.S. model, and the Canadian model along with its prede- cessors. Note that, for damage accumulation models, the accumulated damage α(t) cannot be observed. Instead, we only observe the breaking time Ts and/or Tc. Therefore, we focus on the solutions for Ts and Tc implied by the damage accu- mulation models. We review the solutions for the breaking time Ts in the ramp load test and the breaking time Tc in the constant load test for the U.S. model. We discuss both Gerhards and Link (1987)’s approach and Foschi and Yao (1986)’s approach for the U.S. model. We review the approximate solutions of the breaking time Ts in the ramp load test and the breaking time Tc in the constant load test proposed by Foschi and Yao (1986) for the Canadian model. We discuss Foschi and Yao’s method for parameter estimation in the Canadian model. We discuss the problem of scale, which is a common problem in damage ac- cumulation models. We propose a method to revise the damage accumulation 7 2.2. MADISON CURVE models to eliminate the problem of scale. Table 2.1 summarizes commonly used notation and their definitions. Sections 2.2 through 2.4 present details of some damage accumulation mod- els. Section 2.2 is a brief review of a damage accumulation model based on the Madison Curve proposed by Hendrickson et al. (1987). Section 2.3 contains a re- view of an exponential damage accumulation model (the U.S. model) proposed by Gerhards (1979). Section 2.4 contains a review of a sequence of models includ- ing the Canadian model proposed by Foschi and his collaborators. In this thesis, we will focus on the U.S. model and the Canadian model. Section 2.5 discusses revisions of the damage accumulation models. Notations Definitions Comments α(t) damage accumulated by time t 0≤ α(t)≤ 1 Ts breaking time in the ramp load test α(Ts) = 1 Tc breaking time in the constant load test α(Tc) = 1 τ(t) applied load (applied stress) at time t - τs short term strength - σ(t) applied load ratio at time t σ(t) = τ(t)/τs σ0 threshold load ratio 0≤ σ0 ≤ 1 Table 2.1: Commonly used notation for damage accumulation models. 2.2 Madison Curve In North America, the duration of load effect was initially discussed in Wood’s (1951) work. Unfortunately, this work is unavailable and we rely on Gerhards (1977) for its contents. According to Gerhards, Wood conditioned small clear bending specimens of dried Douglas-fir to achieve 6% and 12% moisture contents and then subjected them to constant load levels ranging from 60% to 95% of their short-term strength. Wood’s data from the constant load test are shown in Fig- ure 2.1, but it is not clear whether that figure presents all data from both moisture contents. 8 2.2. MADISON CURVE Figure 2.1: Data from constant load tests by Wood (1951), reproduced from Gerhards (1977). The x-axis is the breaking time T in hours and the y-axis is the applied load ratio σ in percentage. According to Gerhards (1977), Wood proposed a hyperbolic empirical model: σ = a+bT c, (2.1) where σ is the applied load ratio, T is the breaking time, and a,b and c are un- known parameters. Based on the partial data of Wood’s constant load test and the results of a ramp load test by Markwardt and Liska (1948), Wood estimated the parameters as: â = 18.03, b̂ = 108.4 and ĉ =−0.04635. The equation (2.1) with Wood’s estimates is called the Madison Curve. Gerhards (1977) did not mention which loading rate was used to reach the constant load in Wood’s experiments 9 2.2. MADISON CURVE nor which moisture content level was used to obtain the Madison Curve. Wood’s original plot of the Madison Curve is shown in Figure 2.2. Figure 2.2: Madison Curve with ramp loading test trends and constant load- ing test trends by Wood (1951), reproduced from Gerhards (1977). The x-axis is the breaking time T in seconds and the y-axis is the ap- plied load ratio σ in percentage. According to Gerhards (1977), breaking times of various wood products have been modelled based on the Madison Curve. The validity of the Madison Curve has been questioned for a long time. However, due to its simplicity, the Madison Curve is still used today as the basis for the National Design Specification for Wood Construction (NDS) by the American Wood Council. Based on the Madison Curve, a damage accumulation model for the duration 10 2.3. THE U.S. MODEL of load effect was developed by Hendrickson et al. (1987): dα(t) dt = a{σ(t)−σ0}b+, (2.2) where a,b and σ0 are model parameters. Here, (σ −σ0)+ = max{(σ −σ0),0}, which means that the damage will only accumulate when the applied load ratio is larger than the threshold load ratio, i.e., σ(t)> σ0. It is easy to show that, if for all t, σ(t) = σa > σ0, then for (2.2), the solution for T with α(T ) = 1 satisfies σa = σ0+(1/a)1/bT−1/b, which has the same form of the Madison Curve. Wang (2009) called model (2.2) the model derived from the Madison Curve. We call model (2.2) the Madison Curve model in this thesis. 2.3 The U.S. Model 2.3.1 Model Description The U.S. model, also called the exponential damage rate model (EDRM), was proposed by Gerhards (1979). Based on the assumption that the accumulated damage is an exponential function of the applied load ratio, Gerhards proposed the following model: dα(t) dt = exp{−a+bσ(t)} , (2.3) where a and b are model parameters. Here, b > 0. Some authors consider the parameters are fixed while others consider the parameters are random. The U.S. model has been discussed in Gerhard and Link (1987) as well as Foschi and Yao (1986). Although these represent the U.S. model in the same form in their papers, they actually discuss two different models based on their 11 2.3. THE U.S. MODEL different definitions of the short term strength τs. Gerhards and Link treat the short term strength τs as a board dependent param- eter and assume that τs follows a log-normal distribution with median τm. They do not define τs in terms of any breaking time or load pattern. Foschi and Yao treat the U.S. model in a different way. They also consider the short term strength τs as a board dependent parameter, but they further define τs as the breaking load τ(Ts) of the ramp load test when the loading rate k is set so that the mean breaking time is expected to be around one minute. In both approaches, the breaking time is random since the short term strength τs is random. The Gerhards-Link and Foschi-Yao analyses had some common calculation steps, presented in Section 2.3.2. After that, differences arise in their analysis. Gerhard and Link’s analysis is summarized in Section 2.3.3, Foschi and Yao’s, in Section 2.3.4. 2.3.2 Common Steps in Gerhards-Link and Foschi-Yao Analyses Although Gerhards and Link, and Foschi and Yao define the U.S. model in differ- ent ways, they still go through some common steps to calculate the solution for the breaking time Ts in the ramp load test and the solution for the breaking time Tc in the constant load test. The integrations of the differential equation are exactly the same in their analysis, but the solutions appear in different forms due to their different definitions of the short term strength τs and their different notations for other parameters. In this section, we review the common steps of their calculations for solving the U.S. model for Ts and Tc. For the ramp load test, σ(t) = kt/τs, we can integrate (2.3) to get α(t): α(t) = ∫ t 0 exp(−a+bks/τs)ds = τsbk{exp(−a+bkt/τs)− exp(−a)}. (2.4) 12 2.3. THE U.S. MODEL Based on the fact that α(Ts) = 1, we get the equation: τs bk {exp(−a+bkTs/τs)− exp(−a)}= 1. (2.5) Gerhards and Link solve the above equation (2.5) for Ts in terms of a,b,k and τs. Foschi and Yao solve the above equation (2.5) for Ts with another condition that τs = kTs, so Foschi and Yao’s solution of Ts does not contain τs. The exact forms of the solutions are contained in Section 2.3.3 and 2.3.4. For the constant load test, as we discussed in the previous section, the applied load is τ(t) = kt in the ramp loading part for 0≤ t ≤ T0 and then τ(t) = kT0 = τa in the constant loading part for t > T0. If Ts ≤ T0 or equivalently, if α(T0) ≥ 1, then the board breaks during ramp loading. For the case that Ts > T0, the breaking time Tc will depend on the damage accumulated during the ramp loading part and during the constant loading part. The damage accumulated during the ramp loading part can be calculated from (2.4): α0 = α(T0) = τs bk {exp(−a+bkT0/τs)− exp(−a)} = τs bk {exp(−a+bτa/τs)− exp(−a)}. (2.6) In the constant loading part of the constant load test, τ(t) = τa, so we can integrate (2.3) from T0 to find the damage accumulated during the constant loading part, and then find the total damage accumulated by time t > T0: α(t) =α0+ ∫ t T0 exp(−a+bτa/τs)ds=α0+(t−T0)exp(−a+bτa/τs), for t > T0. (2.7) Setting α(Tc) in (2.7) equal to 1 and solving for Tc yields Tc = T0+(1−α0)/exp(−a+bτa/τs) (2.8) if Tc > T0. 13 2.3. THE U.S. MODEL The above steps for finding Tc are the same as in Gerhards and Link’s analysis as well as Foschi and Yao’s analysis. Gerhards and Link solve for Tc as in (2.8) in terms of a,b,k,T0 and τs. Foschi and Yao solve for τs in terms of a and b, and substitute the result into (2.8), so Foschi and Yao’s solution of Tc does not contain τs. The exact forms of the solutions are contained in Section 2.3.3 and 2.3.4. 2.3.3 Solutions for Ts and Tc from Gerhards and Link’s (1987) Approach Gerhards and Link (1987) assume that τs is log-normally distributed and expressed τs as follows: τs = τm exp(wR) where τm is the median short term strength, w is a (unitless) scale parameter and R is a standard normal random effect. To get their forms of Ts, α0 and Tc, let B = b/τs = b/{τm exp(wR)}. Then we can solve for the breaking time Ts in the ramp load test from (2.5): Ts = ln{Bk exp(a)+1} Bk . (2.9) Here, Ts 6= τs/k. Substituting B for b/τs in (2.6), we can write α0 as α0 = 1 Bk {exp(−a+Bτa)− exp(−a)}. (2.10) Substituting B for b/τs and substituting α0 from (2.10) in (2.8), we can write Tc in the U.S. model as: Tc = ln{Bk exp(a)+1} Bk , if ln{Bk exp(a)+1} Bk ≤ T0, T0− 1Bk + exp(−Bτa) { 1 Bk + exp(a) } , if ln{Bk exp(a)+1} Bk > T0. 14 2.4. THE CANADIAN MODEL AND ITS PREDECESSORS (2.11) 2.3.4 Solutions for Ts and Tc from Foschi and Yao’s (1986) Approach Foschi and Yao (1986) define the short term strength of a board as its breaking strength in the ramp load test with the ramp loading rate k set so that the mean breaking time is around one minute. In other words, they solve for τs from τs = τ(Ts) = kTs. For the ramp load test, replacing τs by kTs in (2.5), we can solve for the break- ing time Ts in the ramp load test: Ts = exp(a)b exp(b)−1 ≡ Ab (2.12) where A≡ exp(a)/{exp(b)−1}. Noticing the fact that τs = kTs, (2.12) is equivalent to: τs = exp(a)bk exp(b)−1 = Abk. (2.13) Substituting τs from (2.13) in (2.6) and (2.8), we can write the solution for Tc in terms of the model parameters A and b: Tc = { Ab, if Ab≤ T0, T0+A{exp(b−T0/A)−1} , if Ab > T0. 2.4 The Canadian Model and its Predecessors 2.4.1 Models Description Foschi and his collaborators proposed three models: Model I in Barrett and Foschi (1978), Model II in Barrett and Foschi (1978, 1982) and the Canadian Model in 15 2.4. THE CANADIAN MODEL AND ITS PREDECESSORS Foschi and Yao (1986). Model I is a generalization of the Madison Curve model, and is given by dα(t) dt = a{σ(t)−σ0}b+αc(t), (2.14) where a,b,c and σ0 are model parameters. In Model I, the damage will accumulate only when τ(t)/τs−σ0 > 0. For the ramp load test with τ(t) = kt, this means that damage can only accumulate at time t bigger than t0 = σ0τs/k. Barrett and Foschi’s model I takes the current damage status α(t) into consid- eration if c 6= 0. The current damage status appears in the model as a multiplier. If c = 0, Barrett and Foschi model I reduces exactly to the Madison Curve model. Model II uses an additive model to include the current damage status. The model is given by dα(t) dt = a{σ(t)−σ0}b++ cα(t), (2.15) where a,b,c and σ0 are random model parameters. Again, if c = 0, then model II reduces to the Madison Curve model. In Model II, the damage will accumulate if either τ(t)/τs−σ0 > 0 or α(t)> 0. If τ(t) is non-decreasing, as in the ramp load test and the constant load test, the damage will only accumulate when τ(t)/τs−σ0 > 0. Barrett and Foschi (1982) conducte an experiment on western hemlock lumber with an approximately 10% moisture content. From these data, they estimate the parameters in (2.15) as â = 0.721568× 10−14 hour−1; b̂ = 34.0; ĉ = 0.150× 10−2 hour−1 and σ̂0 = 0.5. Foschi and Yao (1986) find that inclusion of the second term of Barrett and Foschi’s model II lead to some unreasonable results. For example, if some dam- age was accumulated before some time t0 (e.g., if α(t0) = 0.5) and the load was quite small after t0 (e.g., if σ(t) = 0.01 for t ≥ t0), the damage would still increase exponentially according to model (2.15) since the second term in (2.15) is dom- inant. However, in practice the damage should not increase exponentially in this 16 2.4. THE CANADIAN MODEL AND ITS PREDECESSORS case. To correct this problem, Foschi and Yao (1986) propose a third model which is called the Foschi and Yao model or the Canadian model: dα(t) dt = a{τ(t)−σ0τs}b++ c{τ(t)−σ0τs}n+α(t), (2.16) where a,b,c,n and σ0 are model parameters. This model is called the Canadian model since it is the national standard of Canada. In the Canadian model, the damage will accumulate only when τ(t)/τs−σ0 > 0. Foschi and Yao (1986) assume that the parameters a,b,c,n,σ0 and τs in the Canadian model are all random effects, which means that a,b,c,n,σ0 and τs vary from board to board. There is a dependency among these random effects. In practice, a can be solved in terms of the others. For all boards of one type under the same conditions, the vector (a,b,c,n,σ0,τs) follows the same distribution, no matter what the load pattern is. But for boards of different types or boards of one type under different conditions such as differing moisture contents, the vector (a,b,c,n,σ0,τs) may follow a different distribution. 2.4.2 Solutions for Ts and Tc To get a closed form solution for the breaking time Ts in the ramp load test, Foschi and Yao (1986) disregarded the second term of the differential equation in (2.16) and obtained the following simplified model: dα(t) dt ≈ a{τ(t)−σ0τs}b+. (2.17) We can solve for α(t) explicitly in model (2.17) for the ramp load test with τ(t) = kt, and based on the fact that α(t0) = 0: α(t)≈ ∫ t t0 a(ks−σ0τs)b+ds = a k(b+1) (kt−σ0τs)b+1+ . (2.18) The model in (2.17) is over-parameterized, so Foschi and Yao use (2.18) to 17 2.4. THE CANADIAN MODEL AND ITS PREDECESSORS solve for a in terms of the other board-specific parameters. Specifically, they set α(Ts) = 1 and solve (2.18) for a, yielding a≈ k(b+1){τs(1−σ0)}b+1 . (2.19) By replacing τs with kTs, an equivalent expression of the approximation (2.19) yields: Ts ≈ {k(b+1)/a} 1/(b+1) k(1−σ0) . (2.20) In the constant load test, the applied load is τ(t) = kt in the ramp loading part for 0 ≤ t ≤ T0 and then τ(t) = τa in the constant loading part for t > T0. As we discussed in Section 2.3.2, if Ts ≤ T0 or equivalently, if α(T0)≥ 1, then the board breaks during ramp loading. For the case that Ts > T0, Foschi and Yao propose an approximate method to find Tc, the breaking time during constant loading. First, they use the integration (2.18) of the approximate model (2.17) to calcu- late α0, the approximate accumulated damage during the ramp loading part of the constant load test: α0 = α(T0)≈ ak(b+1)(kT0−σ0τs) b+1 = a k(b+1) (τa−σ0τs)b+1. (2.21) Substituting (2.19) into (2.21) gives an approximation for α0: α0 ≈ ( τa−σ0τs τs−σ0τs )b+1 . (2.22) Second, they solve for α(t)−α(T0), t > T0 by integrating the full Canadian model (2.16) for the constant loading part of the constant load test. The accumu- 18 2.4. THE CANADIAN MODEL AND ITS PREDECESSORS lated damage at time t ≥ T0 is given by α(t) = α0+ ∫ t T0 dα(s) = {α0+(a/c)(τa−σ0τs)b−n}exp{c(τa−σ0τs)n(t−T0)}− (a/c)(τa−σ0τs)b−n. They then find the breaking time Tc for the constant loading part of the constant load test by solving α(Tc) = 1, and obtained: Tc = T0+ 1 c(τa−σ0τs)n ln { 1+(a/c)(τa−σ0τs)b−n α0+(a/c)(τa−σ0τs)b−n } . (2.23) To sum up, Foschi and Yao give an approximate solution for the breaking time Tc in the Canadian model: Tc≈ {k(b+1)/a}1/(b+1) k(1−σ0) , if {k(b+1)/a}1/(b+1) k(1−σ0) ≤ T0, T0+ 1 c(τa−σ0τs)n ln { 1+(a/c)(τa−σ0τs)b−n α0+(a/c)(τa−σ0τs)b−n } , if {k(b+1)/a}1/(b+1) k(1−σ0) > T0, (2.24) with a approximated as in (2.19) and α0 approximated as in (2.22). 2.4.3 Parameter Estimation Foschi and Yao (1986) propose a criterion to estimate parameters in the Canadian model. They treated b,c,n and σ0 as independent random effects which follow log-normal distributions with means µb,µc,µn and µσ0 and standard deviations σb,σc,σn and σσ0 respectively. They treat τs as a random effect, which is inde- pendent of b,c,n and σ0, with a distribution equal to the empirical distribution of the short term strength of the boards from the one minute ramp load test. They treated a as another random effect, which can be solved from others. 19 2.4. THE CANADIAN MODEL AND ITS PREDECESSORS Let φ = (µb,µc,µn,µσ0 ,σb,σc,σn,σσ0), which is the vector of parameters to be estimated. Foschi and Yao use a numerical algorithm to find a value of φ to minimize a simulation-based function. The authors do not define a probability model based function to minimize, but we see that their method can be defined as minimizing an approximation to a function ψ . To define ψ , let Tc(1), Tc(2), · · · , Tc(N) be the order statistics of the observed breaking times in the constant load test and let F = F(·;φ) be the cumulative distribution function of the breaking times assuming that the Canadian model holds, with parameters equal to φ . Then ψ(φ) = N ∑ j=1 {1−F−1( j/N;φ)/Tc( j)}2. (2.25) Foschi and Yao minimize an approximation of (2.25) via a Newton Raph- son algorithm, an iterative method which requires calculation of F(·;φ) eval- uated at the current values of φ , along with calculation of all partial deriva- tives of F(·;φ) with respect to φ . To carry out this algorithm, they approx- imate F by simulating breaking times with distribution F(·;φ). The approxi- mation of F , which we denote by F̂ , is simply the empirical distribution func- tion of the simulated breaking times. For derivative approximation, they simply calculate the change in F̂ over the change in one element of φ when the oth- ers are fixed. For instance, the partial derivative of F at φ with respect to µb was approximated by {F̂(·;φ1)− F̂(·;φ2)}/(0.002µb), where φ1 = (µb+0.001µb, µc,µn,µσ0 ,σb,σc,σn,σσ0) and φ2 = (µb−0.001µb,µc,µn,µσ0 ,σb,σc,σn, σσ0). Foschi and Yao (1986) obtained the following estimates: the means of b,c,n and σ0 are estimated by 35.204,0.1559×10−6,1.429 and 0.578 respectively and the standard deviations of b,c,n and σ0 are estimated by 6.589,0.9621×10−7,0.139 and 0.163 respectively. They do not report the units of the parameters in their pa- per. We should be aware that these are estimates of the means and variances of the log-normal distributed random effects, but not the parameters we usually use to characterize log-normal distributions. Foschi and Yao’s approach to parameter estimation is problematic. The range 20 2.5. THE PROBLEM OF SCALE of σ0 under the log-normal assumption is not the same as the range of σ0 from the experimental perspective. From the experimental perspective, if σ(t) = τ(t)/τs = 1, i.e., the applied load equals the short-term strength, the board will break for sure. This implies (1−σ0)+ > 0 and yields σ0 < 1. However, the probability of getting a random sample larger than 1 from the log-normal distribution with mean 0.578 and standard deviation 0.163 is 0.017, which is not negligible. 2.5 The Problem of Scale In the previous studies, researchers measured the breaking times in different units of measurement and estimated the parameters in different units of measurement. Gerhards and Link (1987) measured the breaking times in minutes and estimated a and b with a unit of log(minute−1), while Foschi and Yao (1986) measured the breaking times in hours and estimated a and b with a unit of log(hour−1). However, the transformation between log(minute−1) and log(hour−1) is not clear, which makes their results incomparable. This problem is caused by the scale of the measurements. We call this problem the problem of scale in this thesis. To solve the problem of scale, we propose a general method to revise the dam- age accumulation models by adding reference levels to the models. We define all parameters in the revised damage accumulation models to be unitless, and we scale the breaking times to be unitless before estimating the model parameters. In this approach, the estimates will be the same regardless of the unit of the measure- ments used in the experiments. First, we consider the U.S. model: dα(t) dt = exp{−a+bσ(t)} . We notice that σ(t) is unitless by definition. If we define both a and b to be unitless, we need to add a unit quantity to the right side of the differential equation since the left side of the differential equation has a unit, which is the inverse of 21 2.5. THE PROBLEM OF SCALE the unit of time. Based on this idea, we revise the U.S. model as: dα(t) dt = λ exp{−a+bσ(t)}, (2.26) where λ = 1 hour−1, and a and b are unitless model parameters. By introduc- ing the unit quantity λ , the solutions for Ts and Tc in the revised U.S. model are changed. Here, we show the solutions for Ts and Tc in the revised Foschi and Yao’s version of the U.S. model. The idea can be applied to Gerhards and Link’s version of the U.S. model without difficulty. We can also choose λ = 1 minute−1, or the inverse of other units of measure- ment of time. The magnitude of the estimates will change accordingly. For the ramp load test, the breaking time Ts is now given by λTs = exp(a)b exp(b)−1 . (2.27) For the constant load test, the breaking time Tc is now given by λTc = { Ab, if Ab≤ λT0, λT0+A{exp(b−λT0/A)−1} , if Ab > λT0, (2.28) where A = exp(a)/{exp(b)−1}. We notice that by adding the unit quantity λ to the model, we scale the break- ing times Ts and Tc to be unitless. In the revised model, the parameter estimates will be consistent regardless of the unit of measurements of the breaking times since both a and b are unitless. Because λ is a unit quantity which is set as 1 hour−1 and does not need to be estimated, adding λ will not complicate the estimation of the parameters. The Madison Cuve model can be revised in the same way as the U.S. model. 22 2.5. THE PROBLEM OF SCALE Now we consider the Canadian model dα(t) dt = a{τ(t)−σ0τs}b++ c{τ(t)−σ0τs}n+α(t), (2.29) where the problem of scale also appears. First, as polynomial functions, the terms {τ(t)−σ0τs}b+ and {τ(t)−σ0τs}n+ in (2.29) need to be unitless for the sake of unit comparability. However, the stress τ(t) and the short term strenth τs have the unit of pressure while the threshold stress ratio σ0 does not have a unit by definition. So the terms {τ(t)−σ0τs}b+ and {τ(t)−σ0τs}n+ are not unitless in the Canadian model (2.29). One possible revision is to replace {τ(t)−σ0τs}b+ and {τ(t)−σ0τs}n+ in (2.29) by the unitless terms {σ(t)−σ0}b+ and {σ(t)−σ0}n+ respectively. Second, the two sides of the differential equation (2.29) should have the same units if the model parameters a,b,c,n and σ0 are considered as unitless. By intro- ducing two unit quantities λ1 and λ2, we revise the Canadian model as dα(t) dt = λ1a{σ(t)−σ0}b++λ2c{σ(t)−σ0}n+α(t), (2.30) where λ1 = λ2 = 1 hour−1 and a,b,c,n and σ0 are unitless parameters. For the ramp load test, the breaking time Ts is now approximated by λ1Ts ≈ b+1a(1−σ0)b+1 . (2.31) For the constant load test, the breaking time Tc is now approximated by Tc≈ λ−11 b+1 a(1−σ0)b+1 , if b+1 a(1−σ0)b+1 ≤ λ1T0, T0+λ−12 1 c(σa−σ0)n ln { 1+(a/c)(σa−σ0)b−n α0+(a/c)(σa−σ0)b−n } , if b+1 a(1−σ0)b+1 > λ1T0, (2.32) 23 2.5. THE PROBLEM OF SCALE where α0 ≈ {(1−σ0)/(σa−σ0)}b+1. In conclusion, the damage accumulation models can be revised by adding one or more unit quantities to solve the problem of scale. After the revising the mod- els, all parameters to be estimated are unitless, and the breaking times are scaled before the estimation procedures, so the estimation results will not depend on the units of the measurements used in the experiments. 24 Chapter 3 Some Results Based on Simulated Data 3.1 Introduction In this section, we perform simulation studies to generate the breaking times in the ramp load test and in the constant load test based on the U.S. model and the Canadian model. The main objective of this simulation study is to investigate the effects of the assumed distributions of the random effects on the breaking times. The simulation study will also provide generated datasets to validate our estimation methods in the next Chapter. For the simulations based on the U.S. model, we treat the model parameters a and b as random effects. We adopt the revised Foschi and Yao’s (1986) version of the U.S. model to get the solution for the breaking time Ts in the ramp load test and the solution for the breaking time Tc in the constant load test, as discussed in Section 2.5. We generate random effects that are Normal and/or log-normal, and study the effects of the mean, standard deviation, and the type of distribution on the distribution of the breaking times. For the simulation based on the Canadian model, we adopt Foschi and Yao’s approach to get the approximate solution for the breaking time Ts in the ramp load 25 3.2. GENERATING DATA BASED ON THE U.S. MODEL test and the solution for the breaking time Tc in the constant load test, as discussed in Section 2.4. We compare the simulated breaking times to the breaking times in Foschi and Yao’s experiments. We also fit the simulated breaking times to the Weibull distri- bution, the log-normal distribution, the Normal distribution, and the exponential distribution to investigate the distribution of generated breaking times. In this chapter, all comparisons are made on the logarithm of the breaking times. 3.2 Generating Data Based on the U.S. Model In this section, we generate the breaking time Ts in the ramp load test and breaking time Tc in the constant load test based on the revised Foschi and Yao’s version of the U.S. model, as shown in (2.26), and compare the distributions of the generated breaking times based on different assumptions. Based on the revised Foschi and Yao’s version of the U.S. model, we actually generate the scaled breaking times λTs and λTc in this study rather than Ts and Tc. However, for simplicity, we still refer to the generated values as Ts and Tc. Section 3.2.1 contains a review of the revised Foschi and Yao’s version of the U.S. model. Section 3.2.2 contains the basic setups for the simulation studies. Section 3.2.3 contains the comparisons of the generated Ts and Tc based on dif- ferent types of distributions of the random effects a and b. Section 3.2.4 contains the comparisons of the generated Ts and Tc based on different means and standard deviations of the random effects a and b. 3.2.1 A Brief Summary of the U.S. Model In Section 2.5, we discussed the revised Foschi and Yao’s version of the U.S. model. The revised Foschi and Yao’s version of the U.S. model is given by dα(t) dt = λ exp{−a+bσ(t)}. (3.1) 26 3.2. GENERATING DATA BASED ON THE U.S. MODEL where λ = 1 hour−1, a and b are unitless model parameters, and τs is the short term strength, which is defined as the breaking strength of the board in the ramp load test when the ramp loading rate k is set so that the mean breaking time is around one minute. The breaking time Ts in the ramp load test is given by: λTs = exp(a)b exp(b)−1 . (3.2) The breaking time Tc in the constant load test is given by: λTc = exp(a)b exp(b)−1 , if exp(a)b exp(b)−1 ≤ λT0, λT0+ exp(a) exp(b)−1 [ exp { b−λT0 exp(b)−1exp(a) } −1 ] , if exp(a)b exp(b)−1 > λT0 (3.3) where T0 is the loading time in the ramp loading part of the constant load test. To generate Ts, we need to specify the distributions of a and b, while to gen- erate Tc, we need to specify the loading time T0 in the loading part of the constant load test, along with the distributions of a and b. T0 is usually set to be the p-th percentile of the generated Ts in the ramp load test. We notice that, in the solution (3.2) for Ts, we do not need to specify the loading rate k, which implies that the distribution of random effects a and b should depend on the loading rate k. Otherwise, equation (3.2) would imply that the breaking times do not depend on the loading rate k, which contradicts the fact that, in experiments, the breaking times are shorter if the loading rate is large. We also notice that the constant load τa does not appear explicitly in the solu- tion (3.3) for Tc in the U.S. model. However, τa is actually included in the solution (3.3) since T0 = τa/k. 27 3.2. GENERATING DATA BASED ON THE U.S. MODEL 3.2.2 Simulation Setups We denote the means of a and b by µa and µb respectively and denote the standard deviations of a and b by σa and σb respectively. In all simulation studies, we assume that a and b are independent. To generate Ts, we generate ns = 1000 a’s and b’s from the assumed distri- butions with means µa and µb respectively and standard deviations σa and σb respectively, and then calculate Ts from (3.2) based on the pairs of a and b. To generate Tc, we choose T0 as the p-th percentile of the generated Ts’s, and then generate Tc’s based on another nc = 1000 a’s and b’s, which are independent of the a’s and b’s we generated for Ts. 3.2.3 The Effect of the Type of Distribution of the Random Effects In this section, we will investigate the effects of the choice of random effects distribution on the generated breaking times Ts and Tc. So far, to our knowl- edge, random effects have not previously been used in the U.S. model. However, researchers have used random effects in other damage accumulation models, as- suming that the random effects follow log-normal distributions, as in Foschi and Yao (1986). In this study, we compare three different types of distributional assumptions on the random effects a and b: (1) a and b are both Normal, (2) a and b are both log-normal, and (3) a is Normal and b is log-normal. We compare the differences of the generated breaking times from the three scenarios when the means and the standard deviations of a and b are the same. For the means µa and µb of the random effects a and b, we use µa = 42 and µb = 50. We choose these values based on our estimates of the real data from Foschi and Yao’s experiments (details of the estimation procedure are in Chapter 6). These values are similar to those in the literature: Gerhards and Link (1986) considered a and b as fixed across all boards and estimated them as â= 43.17 and b̂ = 49.75. 28 3.2. GENERATING DATA BASED ON THE U.S. MODEL a b Breaking times Distribution σa Distribution σb Figure for Ts (Tc) Ts(Tc) Normal 1 Normal 1 Figure 3.1 (3.2) Ts(Tc) Normal 1 log-normal 1 Figure 3.1 (3.2) Ts(Tc) log-normal 1 log-normal 1 Figure 3.1 (3.2) Ts(Tc) Normal 10 Normal 10 Figure 3.3 (3.4) Ts(Tc) Normal 10 log-normal 10 Figure 3.3 (3.4) Ts(Tc) log-normal 10 log-normal 10 Figure 3.3 (3.4) Ts(Tc) Normal 10 Normal 1 Figure 3.5 (3.6) Ts(Tc) Normal 10 log-normal 1 Figure 3.5 (3.6) Ts(Tc) log-normal 10 log-normal 1 Figure 3.5 (3.6) Ts(Tc) Normal 1 Normal 10 Figure 3.7 (3.8) Ts(Tc) Normal 1 log-normal 10 Figure 3.7 (3.8) Ts(Tc) log-normal 1 log-normal 10 Figure 3.7 (3.8) Table 3.1: Summary of the setups of the simulation runs for studying the effect of the type of distribution of the random effects a and b when µa = 42 and µb = 50, and p = 0.2 for the constant load test. Table 3.1 summarizes the setups of the simulation runs for studying the effect of the type of distribution of the random effects a and b as displayed in Figure 3.1 to Figure 3.8. Figure 3.1 to Figure 3.8 are generated when (µa,µb) = (42,50) and p = 0.2 for the constant load test. We use these figures to compare the distributions of Ts and Tc under the three different types of distribution of the random effects. Figure 3.1 and Figure 3.2 show that the distributions of the generated log10(Ts) and log10(Tc) are similar under the three different distributional assumptions when the standard deviations σa and σb are both relatively small. This can be explained by the fact that when σa and σb are both relatively small, the Normal distribution and the log-normal distribution are similar and, as a result, the distributions of the generated breaking times Ts and Tc are similar. Figure 3.3 shows results similar to those in Figure 3.1. The distributions of the generated log10(Ts) are roughly the same under the three different distribu- 29 3.2. GENERATING DATA BASED ON THE U.S. MODEL tional assumptions, although there are observable differences in the lower tail and the upper tail of the distribution of the generated log10(Ts). The differences can be explained by the fact that the log-normal distribution is heavy-tailed, so the distribution of the generated log10(Ts) is heavy-tailed when σa and σb are large. Figure 3.4 shows that the distributions of the generated log10(Tc) are different under the three distributional assumptions, especially in the centre of the distribu- tion. In other words, those generated breaking times Tc which are slightly larger than the loading time T0 are more sensitive to the distributional assumptions of the random effects a and b when σa and σb are large. Figure 3.5 and Figure 3.6 show that the distributions of the generated log10(Ts) and log10(Tc) are mainly affected by the distributional assumptions on a when σa is relatively large and σb is relatively small. The lower left panels of Figure 3.5 and Figure 3.6 show that the distribution of the generated log10(Ts) and the distribution of the generated log10(Tc) both stay roughly the same when the distributional assumption on b changes from Normal to log-normal. However, from the lower right panels of Figure 3.5 and Figure 3.6, the distributions of the breaking times are quite different when the distributional assumption on a changes from Normal to log-normal, which means the distributional assumption on a has an effect on the distribution of log10(Ts) and log10(Tc). Figure 3.7 and Figure 3.8 show that the distributions of the generated log10(Ts) and log10(Tc) are mainly affected by the distributional assumption on b when σa is relatively small and σb is relatively large. The lower right panels of Figure 3.7 and Figure 3.8 show that the distribution of the generated log10(Ts) and the distribution of the generated log10(Tc) both stay roughly the same when the distributional assumption on a changes from Normal to log-normal. However, from the lower left panels of Figure 3.7 and Figure 3.8, the distributions of the generated breaking times are different when the distributional assumption on b changes from Normal to log-normal, which means the distributional assumption on b has an effect on the distribution of log10(Ts) and log10(Tc). In conclusion, the distributions of Ts and Tc are not greatly affected by the 30 3.2. GENERATING DATA BASED ON THE U.S. MODEL −3 −2 −1 0 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 ECDF of log10(Ts) log10(Ts) Pe rc e n ta ge s a and b are Normal a and b are log−normal a is Normal and b is log−normal lll lll lll lll llll lll llll llll llll llll llll llll llll llll llll llll lll llll lll llll lll l ll −3 −2 −1 0 − 3 − 2 − 1 0 QQ plot for log10(Ts) log10(Ts): a and b are Normal lo g 1 0(T s): a a n d b ar e lo g− no rm a l ll lll ll llll lll llll llll llll llll llll llll llll llll llll llll llll llll llll llll llll ll ll −3 −2 −1 0 − 3 − 2 − 1 0 QQ plot for log10(Ts) log10(Ts): a and b are Normal lo g 1 0(T s): a is N or m a l a nd b is lo g− no rm a l lll lll ll lll lll lll llll llll llll llll llll llll llll llll llll llll llll llll llll llll lll l ll −3 −2 −1 0 − 3 − 2 − 1 0 QQ plot for log10(Ts) log10(Ts): a and b are log−normal lo g 1 0(T s): a is N or m a l a nd b is lo g− no rm a l Figure 3.1: Comparisons of the generated Ts in the U.S. model based on the three different distributional assumptions of a and b when µa = 42,µb = 50,σa = 1 and σb = 1. Upper left: the empirical cumulative distributions of log10(Ts) from the three scenarios. Upper right: the quantile-quantile plot of the generated log10(Ts) from the assumption that a and b are both Normal and the assumption that a and b are both log-normal. Lower left: the quantile-quantile plot of the generated log10(Ts) from the assumption that a and b are both Normal and the assumption that a is Normal and b is log-normal. Lower right: the quantile-quantile plot of the generated log10(Ts) from the assumption that a and b are both log-normal and the assumption that a is Normal and b is log-normal. 31 3.2. GENERATING DATA BASED ON THE U.S. MODEL 0 5 10 15 20 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 ECDF of log10(Tc) log10(Tc) Pe rc e n ta ge s a and b are Normal a and b are Log−normal a is Normal and b is Log−normal llll llll llll llll lll lllll llll llll llll lll llll lll lll lll lll llll llll llll llll llll llll llll llll lll lll 0 5 10 15 20 0 5 10 15 20 QQ plot for log10(Tc) log10(Tc): a and b are Normal lo g 1 0(T c): a a n d b ar e lo g− no rm a l llll llll lll llll lll llll llll llll lll lll llll llll lll llll llll llll llll llll llll llll llll llll llll llll l 0 5 10 15 20 0 5 10 15 20 QQ plot for log10(Tc) log10(Tc): a and b are Normal lo g 1 0(T c): a is N or m a l a nd b is lo g− no rm a l lll llll llll llll lll lll llll lll llll lll llll llll llll llll llll lll llll llll llll llll llll llll llll llll lll 0 5 10 15 20 0 5 10 15 20 QQ plot for log10(Tc) log10(Tc): a and b are log−normal lo g 1 0(T c): a is N or m a l a nd b is lo g− no rm a l Figure 3.2: Comparisons of the generated Tc in the U.S. model based on the three different distributional assumptions of a and b when µa = 42,µb = 50,σa = 1, σb = 1 and p = 0.2. The panels are in the same arrangement as in Figure 3.1. 32 3.2. GENERATING DATA BASED ON THE U.S. MODEL −20 −10 0 10 20 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 ECDF of log10(Ts) log10(Ts) Pe rc e n ta ge s a and b are Normal a and b are log−normal a is Normal and b is log−normal l lll lllll llll lll llll lll lll lllll llll llll llll llll llll llll lll llll lll llll llll llll ll ll l l −20 −10 0 10 20 − 20 − 10 0 10 20 QQ plot for log10(Ts) log10(Ts): a and b are Normal lo g 1 0(T s): a a n d b ar e lo g− no rm a l ll ll lll lllll llll lll llll llll lll llll lll llll llll llll llll llll lllll llll llll llll lllll ll l ll −20 −10 0 10 20 − 20 − 10 0 10 QQ plot for log10(Ts) log10(Ts): a and b are Normal lo g 1 0(T s): a is N or m a l a nd b is lo g− no rm a l l l lll lll llll lll llll lll llll llll llll lll llll llll lll llll llll llll llll lll lllll lllll llll ll l l −20 −10 0 10 20 − 20 − 10 0 10 QQ plot for log10(Ts) log10(Ts): a and b are log−normal lo g 1 0(T s): a is N or m a l a nd b is lo g− no rm a l Figure 3.3: Comparisons of the generated Ts in the U.S. model based on the three different distributional assumptions of a and b when µa = 42,µb = 50,σa = 10 and σb = 10. The panels are in the same arrange- ment as in Figure 3.1. 33 3.2. GENERATING DATA BASED ON THE U.S. MODEL −20 −10 0 10 20 30 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 ECDF of log10(Tc) log10(Tc) Pe rc e n ta ge s a and b are Normal a and b are Log−normal a is Normal and b is Log−normal ll lll lllll llll lll llll llllll ll ll ll l ll ll lll lll lll lll lllll llll llll llll lll lll lll lll ll ll −20 −10 0 10 20 30 − 20 0 10 20 30 QQ plot for log10(Tc) log10(Tc): a and b are Normal lo g 1 0(T c): a a n d b ar e lo g− no rm a l lll llll llll lllll lll llll lllll ll llll l ll ll llll llll lll llll llll llll llll llll llll llll llll ll −20 −10 0 10 20 30 − 20 0 10 20 30 QQ plot for log10(Tc) log10(Tc): a and b are Normal lo g 1 0(T c): a is N or m a l a nd b is lo g− no rm a l lll llll lll lll llll llll lll l lll ll ll ll llllll lll llll lllll llll llll llll llll llll lllll lllll ll −20 −10 0 10 20 30 − 20 0 10 20 30 QQ plot for log10(Tc) log10(Tc): a and b are log−normal lo g 1 0(T c): a is N or m a l a nd b is lo g− no rm a l Figure 3.4: Comparisons of the generated Tc in the U.S. model based on the three different distributional assumptions of a and b when µa = 42,µb = 50,σa = 10, σb = 10 and p= 0.2. The panels are in the same arrangement as in Figure 3.1. 34 3.2. GENERATING DATA BASED ON THE U.S. MODEL −15 −10 −5 0 5 10 15 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 ECDF of log10(Ts) log10(Ts) Pe rc e n ta ge s a and b are Normal a and b are log−normal a is Normal and b is log−normal llll llllll llllll llllll lllll lllll lllll llll llll llll lll llll lll lll lll lll lll lll ll l l l −15 −10 −5 0 5 10 15 − 10 0 5 10 20 QQ plot for log10(Ts) log10(Ts): a and b are Normal lo g 1 0(T s): a a n d b ar e lo g− no rm a l ll ll llll lll lll llll llll llll llll llll llll llll llll llll llll llll llll lll lll lll ll l l −15 −10 −5 0 5 10 15 − 15 − 5 0 5 10 15 QQ plot for log10(Ts) log10(Ts): a and b are Normal lo g 1 0(T s): a is N or m a l a nd b is lo g− no rm a l ll ll lll lll lll lll lll lll lll lll lll lll lll llll llll llll llll llll lllll lll l lll l l −10 −5 0 5 10 15 20 − 15 − 5 0 5 10 15 QQ plot for log10(Ts) log10(Ts): a and b are log−normal lo g 1 0(T s): a is N or m a l a nd b is lo g− no rm a l Figure 3.5: Comparisons of the generated Ts in the U.S. model based on the three different distributional assumptions of a and b when µa = 42,µb = 50,σa = 10 and σb = 1. The panels are in the same arrange- ment as in Figure 3.1. 35 3.2. GENERATING DATA BASED ON THE U.S. MODEL −10 0 10 20 30 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 ECDF of log10(Tc) log10(Tc) Pe rc e n ta ge s a and b are Normal a and b are Log−normal a is Normal and b is Log−normal llll lllll lllll lllll lll l l l l lll llll lll lll lll lll lll llll llll llll lll lll lll lll lll ll l l −10 0 10 20 30 − 10 0 10 20 30 QQ plot for log10(Tc) log10(Tc): a and b are Normal lo g 1 0(T c): a a n d b ar e lo g− no rm a l ll llll lll llll llll ll lll l ll lll lll ll lll llll lll llll llll llll llll llll llll lll lll −10 0 10 20 30 − 10 0 10 20 30 QQ plot for log10(Tc) log10(Tc): a and b are Normal lo g 1 0(T c): a is N or m a l a nd b is lo g− no rm a l ll lll llll lll llll ll lll ll llll ll lll ll l lllll lllll lllll llll llll llll lllll lllll lllll lllll l l −10 0 10 20 30 − 10 0 10 20 30 QQ plot for log10(Tc) log10(Tc): a and b are log−normal lo g 1 0(T c): a is N or m a l a nd b is lo g− no rm a l Figure 3.6: Comparisons of the generated Tc in the U.S. model based on the three different distributional assumptions of a and b when µa = 42,µb = 50,σa = 10, σb = 1 and p = 0.2. The panels are in the same arrangement as in Figure 3.1. 36 3.2. GENERATING DATA BASED ON THE U.S. MODEL −15 −10 −5 0 5 10 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 ECDF of log10(Ts) log10(Ts) Pe rc e n ta ge s a and b are Normal a and b are log−normal a is Normal and b is log−normal l l ll lll lll lll lll lll lll llll llll llll llll llll lllll lllll lllll lllll lllll lllll llll ll l −15 −10 −5 0 5 10 − 20 − 10 0 5 10 QQ plot for log10(Ts) log10(Ts): a and b are Normal lo g 1 0(T s): a a n d b ar e lo g− no rm a l l l ll lll lll lll lll lll lll llll llll llll llll llll lllll lllll lllll lllll lllll lllll llll ll l −15 −10 −5 0 5 10 − 20 − 10 0 5 10 QQ plot for log10(Ts) log10(Ts): a and b are Normal lo g 1 0(T s): a is N or m a l a nd b is lo g− no rm a l l l l ll lll lll lll lll lll lll llll llll llll llll llll llll llll llll llll llll llll l l −20 −15 −10 −5 0 5 10 − 20 − 10 0 5 10 QQ plot for log10(Ts) log10(Ts): a and b are log−normal lo g 1 0(T s): a is N or m a l a nd b is lo g− no rm a l Figure 3.7: Comparisons of the generated Ts in the U.S. model based on the three different distributional assumptions of a and b when µa = 42,µb = 50,σa = 1 and σb = 10. The panels are in the same arrange- ment as in Figure 3.1. 37 3.2. GENERATING DATA BASED ON THE U.S. MODEL −15 −5 0 5 10 15 20 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 ECDF of log10(Tc) log10(Tc) Pe rc e n ta ge s a and b are Normal a and b are Log−normal a is Normal and b is Log−normal l lll lll lll llll llll lll llll l ll lll ll lll l l l l ll ll lllll ll ll lllll lllll lllll ll −15 −5 0 5 10 15 20 − 20 − 10 0 10 20 QQ plot for log10(Tc) log10(Tc): a and b are Normal lo g 1 0(T c): a a n d b ar e lo g− no rm a l l lll lll lll llll llll lll llll l ll lll l llll l l l ll ll llllll ll lll llllll lllll llll −15 −5 0 5 10 15 20 − 20 − 10 0 10 20 QQ plot for log10(Tc) log10(Tc): a and b are Normal lo g 1 0(T c): a is N or m a l a nd b is lo g− no rm a l l ll lll lll llll llll llll llll lll ll ll lll ll l l ll llll lll lll llll llll l −20 −10 0 10 20 − 20 − 10 0 10 20 QQ plot for log10(Tc) log10(Tc): a and b are log−normal lo g 1 0(T c): a is N or m a l a nd b is lo g− no rm a l Figure 3.8: Comparisons of the generated Tc in the U.S. model based on the three different distributional assumptions of a and b when µa = 43,µb = 50,σa = 1, σb = 10 and p = 0.2. The panels are in the same arrangement as in Figure 3.1. 38 3.2. GENERATING DATA BASED ON THE U.S. MODEL distributional assumptions if the standard deviations of both random effects a and b are relatively small. The distributions of Ts and Tc are mainly affected by the distributional assumption on the random effect with a large standard deviation, and not greatly affected by the distributional assumption on the random effect with a small standard deviation. 3.2.4 The Effect of the Mean and the Standard Deviation of the Random Effects In this section, we assume that both a and b are Normal, and investigate the effects of the means and the standard deviations of a and b on the breaking times Ts and Tc. We investigate the effects of the means of the random effects a and b on the breaking times Ts and Tc by fixing the standard deviations of a and b. First, we compare the empirical cumulative distributions of the generated Ts’s and Tc’s as a function of µa with σa,σb and µb fixed. Second, we compare the empirical cumulative distributions of the generated Ts’s and Tc’s as a function of µb with σa,σb and µa fixed. Similarly, we investigate the effects of the standard deviations of the random effects of a and b on the breaking times Ts and Tc by fixing the means of a and b. Table 3.2 summarizes the setups of the simulation runs for studying the effect of the mean and the standard deviation of the random effects a and b when a and b are both Normal, as displayed in Figure 3.9 and Figure 3.10. Figure 3.9 compares the empirical cumulative distributions of the generated log10(Ts) and the empirical cumulative distributions of the generated log10(Tc) based on the assumption that a and b are both Normal, σa = 1 and σb = 1 are fixed, and either µa or µb is not fixed. Figure 3.10 shows the comparisons of the empirical cumulative distributions of the generated log10(Ts) and the empirical cumulative distributions of the generated log10(Tc) based on the assumption that a and b are both Normal, µa = 42 and µb = 50 are fixed, and either σa or σb is not fixed. 39 3.2. GENERATING DATA BASED ON THE U.S. MODEL a b T µa σa µb σb Figure Ts and Tc 37 1 50 1 Figure 3.9 Ts and Tc 42 1 50 1 Figure 3.9 Ts and Tc 47 1 50 1 Figure 3.9 Ts and Tc 42 1 45 1 Figure 3.9 Ts and Tc 42 1 50 1 Figure 3.9 Ts and Tc 42 1 55 1 Figure 3.9 Ts and Tc 42 1 50 1 Figure 3.10 Ts and Tc 42 3 50 1 Figure 3.10 Ts and Tc 42 5 50 1 Figure 3.10 Ts and Tc 42 1 50 1 Figure 3.10 Ts and Tc 42 1 50 3 Figure 3.10 Ts and Tc 42 1 50 5 Figure 3.10 Table 3.2: Summary of the setups of the simulation runs for studying the effect of the mean and the standard deviation of the random effects a and b when a and b are both Normal, and p = 0.2 in the constant load test. Figure 3.9 shows that the means of the random effects a and b mainly affect the locations of the distributions of the breaking times. In the upper left panel, upper right panel and the lower left panel, the location shifts are uniform over the range of the generated breaking times. In the lower right panel, the location shift effect is significant in the lower part of the distribution of the generated Tc, but is small in the upper part of the distribution of the generated Tc. This panel shows that the magnitude of µb affects the breaking time more for small values of Tc, and does not affect the maximum breaking time Tc in the constant load test. The two upper panels also show that the generated Ts’s and Tc’s increase when µa increases, and the two lower panels show that the generated Ts’s and Tc’s decrease when µb increases since log(Ts) = a+ log(b)− log{exp(b)− 1}, the effect of µa on Ts’s distribution is clear. Since exp(b) dominates b, the effect of µb is also clear. For Tc, this cannot be explained easily since the solution (3.3) is complicated. From 40 3.2. GENERATING DATA BASED ON THE U.S. MODEL −8 −6 −4 −2 0 2 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 ECDF of log10(Ts), µb=50 log10(Ts) Pe rc e n ta ge s µa=37 µa=42 µa=47 −5 0 5 10 15 20 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 ECDF of log10(Tc), µb=50 log10(Tc) Pe rc e n ta ge s µa= 37 µa= 42 µa= 47 −8 −6 −4 −2 0 2 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 ECDF of log10(Ts), µa=42 log10(Ts) Pe rc e n ta ge s µb= 45 µb= 50 µb= 55 −5 0 5 10 15 20 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 ECDF of log10(Tc), µa=42 log10(Tc) Pe rc e n ta ge s µb= 45 µb= 50 µb= 55 Figure 3.9: Comparisons of the empirical cumulative distributions of the generated Ts and Tc from the U.S. model based the assumption that a and b are both Normal, σa = 1 and σb = 1. One mean is fixed at the value stated in the plot title. Figure 3.9, we also notice that in order to generate reasonable breaking times Ts, i.e., with the mean of Ts’s around 1 minute, µa−µb should be between 5 and 10. Figure 3.10 shows that the standard deviations of the random effects a and b mainly affect the shapes of the distributions of the breaking times. The variation among the generated Ts’s and the variation among the generated Tc’s increase 41 3.2. GENERATING DATA BASED ON THE U.S. MODEL −8 −6 −4 −2 0 2 4 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 ECDF of log10(Ts), σb=1 log10(Ts) Pe rc e n ta ge s σa= 1 σa= 3 σa= 5 −10 −5 0 5 10 15 20 25 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 ECDF of log10(Tc), σb=1 log10(Tc) Pe rc e n ta ge s σa= 1 σa= 3 σa= 5 −8 −6 −4 −2 0 2 4 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 ECDF of log10(Ts), σa=1 log10(Ts) Pe rc e n ta ge s σb= 1 σb= 3 σb= 5 −10 −5 0 5 10 15 20 25 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 ECDF of log10(Tc), σb=1 log10(Tc) Pe rc e n ta ge s σb= 1 σb= 3 σb= 5 Figure 3.10: Comparisons of the empirical cumulative distributions of the generated Ts and Tc from the U.S. model based the assumption that a and b are both Normal, µa = 42 and µb = 50. One standard variation is fixed at the value stated in the plot title. 42 3.3. GENERATING DATA BASED ON THE CANADIAN MODEL when either σa or σb increases. The breaking times Tc’s in the constant load test concentrate more in the upper tail when either σa or σb is large. In other words, if σa and σb are small, we would expect more boards to break during the early stage of the constant load test, and if either σa or σb is large, we would expect few boards to break during the early stage of the constant load test and more boards to break at almost the same time at the late stage in the constant load test. In conclusion, the locations of the distributions of the generated Ts and Tc are mainly affected by the means of the random effects a and b, and the shapes of the distributions of the generated Ts and Tc are mainly affected by the standard deviations of the random effects a and b. 3.3 Generating Data Based on the Canadian model In this section, we generate the breaking time Tc in the constant load test based on the Canadian model using the parameter estimates from Foschi and Yao (1986). In Section 2.4 and Section 2.5, we discussed the Canadian model and the revised Canadian model. The forms of the solutions for the breaking times of the Canadian model and those of the revised Canadian model are different, unlike the solutions for the breaking times of the U.S. model and those of the revised U.S. model, which only differ by a unit quantity λ . To use the parameter estimates from Foschi and Yao (1986), we adopt the Canadian model in this section for generating breaking times because Foschi and Yao estimated the parameters using the Canadian model. Section 3.3.1 contains a review of the Canadian model. Section 3.3.2 contains the basic setups for the simulation studies. 3.3.1 A Brief Summary of the Canadian Model The Canadian model is given by dα(t) dt = a{τ(t)−σ0τs}b++ c{τ(t)−σ0τs}n+α(t), (3.4) 43 3.3. GENERATING DATA BASED ON THE CANADIAN MODEL where a,b,c,n and σ0 are model parameters. For the ramp load test, the breaking time Ts is approximated by Ts ≈ {k(b+1)/a} 1/(b+1) k(1−σ0) . (3.5) For the constant load test, the breaking time Tc is approximated by Tc≈ {k(b+1)/a}1/(b+1) k(1−σ0) , if {k(b+1)/a}1/(b+1) k(1−σ0) ≤ T0, T0+ 1 c(τa−σ0τs)n ln { 1+(a/c)(τa−σ0τs)b−n α0+(a/c)(τa−σ0τs)b−n } , if {k(b+1)/a}1/(b+1) k(1−σ0) > T0, (3.6) with a approximated by a≈ k(b+1){τs(1−σ0)}b+1 , (3.7) and α0 approximated by α0 ≈ ( τa−σ0τs τs−σ0τs )b+1 . (3.8) We discussed Foschi and Yao’s method of parameter estimation in Section 2.4.3. They treated b,c,n and σ0 as independent random effects which follow log-normal distributions with means µb,µc,µn and µσ0 and standard deviations σb,σc,σn and σσ0 respectively. They treated τs as a random effect, which is in- dependent of b,c,n and σ0, with distribution equal to the empirical distribution of the short term strength of the boards from the one minute ramp load test. They treated a as another random effect, but did not need to specify its distribution since a can be expressed in terms of the other random effects. They obtained the following estimates: the means of b,c,n and σ0 are es- timated by 35.204,0.1559× 10−6,1.429 and 0.578 respectively and the standard 44 3.3. GENERATING DATA BASED ON THE CANADIAN MODEL deviations of b,c,n and σ0 are estimated by 6.589,0.9621×10−7,0.139 and 0.163 respectively. We should be aware that these are estimates of the means and vari- ances of the log-normal distributed random effects, but not the parameters we usually use to characterize log-normal distributions. 3.3.2 Simulation Setups In this section, we describe the setups of generating breaking times based on the Canadian model. As discussed in the previous section, we have the estimates of the distributions of the random effects b,c,n and σ0 from Foschi and Yao (1986). We also have the loading rate k = 6474×60 psi/hour from Foschi and Barrett (1982). We need the distribution of either τs or a to generate the breaking time Tc in the constant load test. These distributions are not available to us. There are two alternatives for generating the breaking time Tc. We can gen- erate τs using the U.S. model and then calculate a from (3.7) given τs, or we can generate a based on some distribution and then calculate τs from (3.7) given a. We adopt the first approach in this section since we do not have much information on the distribution of a. The procedures for generating Tc based on the Canadian model are described as follows. First, we generate ns = 1000 Ts’s from the U.S. model when (µa,µb,σa,σb)= (42,50,0.4,0.4) and we calculate τs using the fact that Ts = τs/k in Foschi and Yao’s version of the U.S. model. Second, we generate nc = 1000 b’s, c’s, n’s and σ0’s from the log-normal distributions with the means µb = 35.204, µc = 0.1559×10−6, µn = 1.429 and µσ0 = 0.578 respectively and the standard devia- tions σb = 6.589, σc = 0.9621×10−7, σn = 0.139 and σσ0 = 0.163 respectively. Third, we calculate a from (3.7) and α0 from (3.8) for the vector (τs,b,c,n,σ0), and then calculate the breaking time Tc from (3.6). The breaking times Tc’s generated based on the Canadian model by this setup are discussed in Section 3.4 and 3.5. 45 3.4. COMPARISONS OF THE GENERATED DATAWITH FOSCHI AND YAO’S DATA 3.4 Comparisons of the Generated Data with Foschi and Yao’s Data In this section, we compare the generated breaking times based on the U.S. model and the Canadian model from the previous sections with the breaking times from Foschi and Yao’s experiments. We do not discuss parameter estimation in this section. The empirical distribution of the logarithm of the breaking times Tc’s from Foschi and Yao’s experiments are shown in Figure 3.11. According to Foschi and Yao (1986), the constant load was set to be the 20-th percentile of the short term strength in their ramp load test. We do not have the data from their ramp load test. By trial and error, we find that the empirical distribution of the generated breaking time Tc is close to the empirical distribution of the breaking time Tc from Foschi and Yao’s experiments when (µa,µb,σa,σb) = (42,50,0.4,0.4) in the U.S. model, as shown in the left panel of Figure 3.11. The empirical distribution of the generated breaking time Tc based on the Canadian model is shown in the right panel of Figure 3.11 for comparison. The left panel of Figure 3.11 shows that the empirical distribution of the gen- erated breaking time Tc based on the U.S. model and the empirical distribution of the breaking time Tc from Foschi and Yao’s experiments do not seem consistent overall. However, the two curves are in close agreement over the ramp loading part, i.e., when Tc ≤ T0, but differ in the constant load part, i.e., when Tc > T0. In particular, in the dataset, the percentage of boards breaking increases rapidly for the times with log10(Tc)> 2. The right panel of Figure 3.11 also shows that the empirical distribution of the generated breaking times Tc’s based on the Canadian model differs from the empirical distribution of the breaking times Tc’s from Foschi-Yao experiments. Again, the two curves are close during the ramp loading part, but different during the constant load part. The percentage of the generated boards that break shows the same rapid increase as does the percentage of the breaking boards from the Foschi-Yao experiments. However, the rapid increase occurs later in the generated 46 3.4. COMPARISONS OF THE GENERATED DATAWITH FOSCHI AND YAO’S DATA −2 0 2 4 6 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 log10(Tc) Pe rc e n ta ge s Foschi and Yao' data Generated data (U.S. model) −2 0 2 4 6 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 log10(Tc) Pe rc e n ta ge s Foschi and Yao' data Generated data (Canadian model) Figure 3.11: Comparisons of the empirical cumulative distributions of breaking times Tc’s from Foschi and Yao’s experiments and the gener- ated Tc based on the U.S. model (left panel), and the Canadian model (right panel). boards. The generated Tc’s based on the revised U.S. model or the Canadian model do not fit the data from the Foschi-Yao experiments well. In an attempt to improve the fit, we tried various values for µa, µb, σa and σb in the U.S. model and var- ious values for k in the Canadian model. For the values we tried, the empirical distribution function of the breaking times Tc’s generated from the U.S. model did not show the rapid increase at the same stage as the data from the Foschi-Yao experiments. The empirical distribution function of the breaking times Tc’s gen- erated from the Canadian model is closer to that of the data from the Foschi-Yao experiments when k is large. 47 3.5. FITTING THE GENERATED BREAKING TIMES TO DISTRIBUTIONS 3.5 Fitting the Generated Breaking Times to Distributions In this section, we investigate the type of distribution of the generated breaking times based on the U.S. model and the Canadian model. We simulate breaking times from each model as described in Sections 3.2.2 and 3.3.2. We then fit the Weibull distribution, the log-normal distribution, the Normal distribution and the exponential distribution to the simulated data. In this section, we study the breaking times T = Ts in the ramp load test, and T = Tc−T0 for Tc > T0 in the constant load test. We do not separately consider the distribution of Tc for Tc ≤ T0 because by definition, P(Tc ≤ t) = P(Ts ≤ t) whenever t ≤ T0. We assess the fits of the parametric distributions as follows. First, we cal- culate the maximum likelihood estimates of the parameters of the distribution to be fitted using parametric models and the R-function survreg. For the breaking times generated from the U.S. model, there is no censoring, but for the breaking times generated from the Canadian model, there is censoring since some boards do not break during the constant load test. To study how well the fitted distribution matches that of the generated breaking times, we simulate n = 1000 Tf it’s from the fitted distribution. Then we draw the quantile-quantile plot of the logarithm of the generated breaking times T ’s from the damage accumulation models and the logarithm of the simulated Tf it’s from the fitted distribution. The quantile-quantile plots are shown in Figure 3.12 to Figure 3.16. Table 3.3 summarizes the setups of the simulation runs in Figure 3.12 to Fig- ure 3.16. In all scenarios, the Normal distribution and the exponential distribution fit the generated breaking times poorly, so we only show one set of quantile-quantile plots for the Normal distribution and the exponential distribution. See Figure 3.12. The fitting results for the Weibull distribution and the log-normal distribution are shown in all figures in this section. Figure 3.12 shows that the log-normal distribution fits generated breaking 48 3.5. FITTING THE GENERATED BREAKING TIMES TO DISTRIBUTIONS T Model Fitted Distribution σa σb Figure Ts U.S. Weibull, log-normal, Normal and exponential 0.4 0.4 Figure 3.12 Ts U.S. Weibull and log-normal 10 10 Figure 3.13 Tc U.S. Weibull and log-normal 0.4 0.4 Figure 3.14 Tc U.S. Weibull and log-normal 10 10 Figure 3.15 Tc Canadian Weibull and log-normal 0.4 0.4 Figure 3.16 Table 3.3: Summary of the setups of the simulation runs for studying the distribution of the generated breaking times when µa = 42 and µb = 50, and p = 0.2 in the constant load test. times Ts’s very well although there are slight differences in the lower tail and the upper tail. The other types of distribution do not provide reasonable fits for the generated breaking times. Figure 3.13 shows similar results as in Figure 3.12. The log-normal distribu- tion fits the generated breaking times Ts’s very well when the standard deviations of the random effects a and b are large, but the Weibull distribution fits the gener- ated breaking times Ts’s poorly. Figure 3.14 shows that neither the Weibull distribution nor the log-normal distribution fits the generated breaking times (Tc− T0)’s for those Tc > T0 well as a whole. The left panel of Figure 3.14 shows that the Weibull distribution is heavier than the distribution of (Tc−T0) in the lower tail, but provides a reasonable fit for other (Tc−T0)’s. The right panel of Figure 3.14 shows that the log-normal distribution fits the distribution of (Tc−T0) well in the centre of the distribution, but is heavier than the distribution of (Tc−T0) in the lower tail and the upper tail. Figure 3.15, where the standard deviations of the random effects a and b are large, shows similar results as in Figure 3.14. The left panel of Figure 3.15 shows that the Weibull distribution fits the distribution of (Tc−T0) reasonably well al- though there are still minor differences in the lower tail and the centre. The right panel of Figure 3.15 shows that the fitted log-normal distribution is lighter than the distribution of (Tc−T0) in the lower tail and heavier in the upper tail. Figure 3.16 shows similar results as in Figure 3.14 and Figure 3.15 when the 49 3.5. FITTING THE GENERATED BREAKING TIMES TO DISTRIBUTIONS breaking times Tc’s are generated from the Canadian model. The Weibull distri- bution fits the data (Tc−T0)’s well except for the lower tail. In conclusion, if the breaking times Ts’s are generated from the U.S. model under the conditions described, the Ts’s seem to follow a log-normal distribution and if the breaking times Tc’s are generated from the U.S. model or the Canadian model under the conditions described, all (Tc− T0)’s which satisfy Tc > T0 and Tc−T0 not close to 0 seem to follow a Weibull distribution. l l l lll lll llll llll llll llll lll llll lllll lllll lllllll lllllll llllllll lllllllll lllllll ll ll −5.5 −4.5 −3.5 −2.5 − 9 − 8 − 7 − 6 − 5 − 4 − 3 Weibull log(Ts) fit te d di st rib u tio n lll llll llll lll lll llll lll llll llll lllll lllll llll llll lllll llll lll llll lllll llll lll ll l l −5.5 −4.5 −3.5 −2.5 − 5 − 4 − 3 − 2 Log−normal log(Ts) fit te d di st rib u tio n l l ll ll lll llll lll lll ll llll lll llll lllll llllll llllll llllllll llllllllll lllllllllll lllllll l l −5.5 −4.5 −3.5 −2.5 − 9 − 7 − 5 − 3 Normal log(Ts) fit te d di st rib u tio n ll ll lll ll ll lll lll lll lll lll lll lll llll lll llll llll lllll llllll llllll llllll llllll llllll lll l l −5.5 −4.5 −3.5 −2.5 − 10 − 8 − 6 − 4 − 2 Exponential log(Ts) fit te d di st rib u tio n Figure 3.12: The quantile-quantile plots of the logarithm of the generated breaking times Ts based on the U.S. model when (µa,µb,σa,σb) = (42,50, 0.4,0.4) and the logarithm of simulated data from the fitted distributions. The type of fitted distribution is shown in the title of the plot. 50 3.5. FITTING THE GENERATED BREAKING TIMES TO DISTRIBUTIONS l l ll lll llll llll llll lll lllll lllll llllllll llllllll llllllllll llll ll −40 −20 0 20 40 − 10 0 − 50 0 Weibull log(Ts) fit te d di st rib u tio n lll lllll llll llll llll llll lllll llll llll lllll llll llll llllll llll lll l l −40 −20 0 20 40 − 40 0 20 Log−normal log(Ts) fit te d di st rib u tio n Figure 3.13: The quantile-quantile plots of the logarithm of the generated breaking times Ts based on the U.S. model when (µa,µb,σa,σb) = (42,50, 10,10) and the logarithm of simulated data from the fitted distributions. The type of fitted distribution is shown in the title of the plot. l l ll lll lll llll lllll lllllll lllllllll llllllll lllllllllll lllllll llllllllll ll −10 0 10 20 30 40 − 60 − 20 20 Weibull log(Tc − T0) fit te d di st rib u tio n ll lll lll lll lllll lllllll lllllll llllll llllll llllll llllll lll llll lll lllll ll l −10 0 10 20 30 40 − 20 0 20 40 Log−normal log(Tc − T0) fit te d di st rib u tio n Figure 3.14: The quantile-quantile plots of the logarithm of the generated breaking times Tc based on the U.S. model when (µa,µb,σa,σb) = (42,50, 0.4,0.4) and p = 0.2, and the logarithm of simulated data from the fitted distributions. The fitted distribution is shown in the title of the plot. 51 3.5. FITTING THE GENERATED BREAKING TIMES TO DISTRIBUTIONS l l l ll llll llll lllllll lllll llll llll llll llllll llllll llllll lll −20 0 20 40 60 − 40 0 40 Weibull log(Tc − T0) fit te d di st rib u tio n lll l llllll lllllll lllll llll llll lll lll llll lll llll lll lll llll lll ll l −20 0 20 40 60 20 40 60 80 Log−normal log(Tc − T0) fit te d di st rib u tio n Figure 3.15: The quantile-quantile plots of the logarithm of the generated breaking times Tc based on the U.S. model when (µa,µb,σa,σb) = (42,50, 10,10) and p= 0.2, and the logarithm of simulated data from the fitted distributions. The type of fitted distribution is shown in the title of the plot. l l l ll ll llll lllll llllll llllll llllll lllll lllll llllll lll l −10 0 5 10 15 − 30 − 10 10 Weibull log(Tc − T0) fit te d di st rib u tio n ll lll lllll llll lllllll llllll lllll lll lll llll lll lll llll lll ll l −10 0 5 10 15 − 10 0 10 20 Log−normal log(Tc − T0) fit te d di st rib u tio n Figure 3.16: The quantile-quantile plots of the logarithm of the generated breaking times Tc based on the Canadian model and the logarithm of simulated data from the fitted distributions. The type of fitted distri- bution is shown in the title of the plot. 52 Chapter 4 Estimation Methods 4.1 Introduction In this section, we propose several methods to estimate parameters in the U.S. model for the ramp load test and for the constant load test. In the ramp load test, the loading rate k is set so that the mean breaking time is around one minute. In the constant load test, the constant load is set to be the p-th (p < 0.5) percentile of the breaking load in the ramp load test. We choose p = 0.2 as an illustration in this chapter. We propose several methods to estimate parameters in the U.S. model. The first one is the maximum likelihood method. However, the likelihood functions seem complicated. The other methods rely on approximations to the solutions of the breaking times Ts and Tc. We also consider approximate likelihoods based on the approximate solutions. We also propose a quantile method based on the approximate solutions of the breaking times. Finally, we consider a method that combines approximate maximum likeli- hoods and the approximate quantile method. Section 4.2 gives some general notation. Section 4.3 contains a description of the maximum likelihood method. Section 4.4 discusses approximations to the 53 4.2. MODELS AND NOTATION solutions for the breaking times. Section 4.5 describes the quantile method. Sec- tion 4.6 describes a method that combines the approximate maximum likelihood method and the quantile method. 4.2 Models and Notation In Section 2.5, we discussed the revised Foschi and Yao’s version of the U.S. model, which is given by dα(t) dt = λ exp{−a+bσ(t)}. where λ = 1 hour−1, a and b are unitless model parameters and b > 0. The breaking time Ts in the ramp load test is given by: λTs = exp(a)b exp(b)−1 . The breaking time Tc in the constant load test is given by: λTc = exp(a)b exp(b)−1 , if exp(a)b exp(b)−1 ≤ λT0, λT0+ exp(a) exp(b)−1 [ exp { b−λT0 exp(b)−1exp(a) } −1 ] , if exp(a)b exp(b)−1 > λT0 where T0 is the loading time in the ramp loading part of the constant load test. In this chapter, we write Ts and Tc instead of λTs and λTc for simplicity since λ = 1 hour−1. The observed data are T̃s = (Ts,1,Ts,2, · · · ,Ts,ns) from the ramp load test, and T̃c = (Tc,1,Tc,2, · · · ,Tc,nc) from the constant load test. We assume that a and b are two random variables following some distribu- tions with means µa and µb, and variance σ2a and σ2b respectively. Let θ = (µa,µb,σ2a ,σ2b ) be the parameters to be estimated. We assume that a and b are independent in this thesis. 54 4.3. LIKELIHOODMETHOD We define the random variables X and Y as follows: X = exp(a)b exp(b)−1 , (4.1) Y = T0+ exp(a) exp(b)−1 [ exp { b−T0 exp(b)−1exp(a) } −1 ] = T0+ X b { exp ( b−T0 bX ) −1 } . (4.2) Then the random variable Ts can be written as Ts = X , (4.3) and the random variable Tc can be written as Tc = { X , if X ≤ T0, Y, if X > T0. (4.4) 4.3 Likelihood Method In this section, we propose a parameter estimation method based on the likelihood method. We consider of collection the data in three different ways. One way is just the ramp load test, another, the constant load test with T0 pre-determined independently from the constant load test, and the third, the ramp load test and the constant load test with T0 depending on the ramp load test. Let fTs and fTc be the density functions of the random variable Ts and Tc re- spectively. The likelihood function for the ramp load test is given by Ls(θ |T̃s) = P(T̃s|θ) = ns ∏ i=1 fTs(Ts,i|θ). (4.5) 55 4.3. LIKELIHOODMETHOD The likelihood function for the constant load test (when T0 is determined in- dependently of T̃c) is given by Lc(θ |T̃c,T0) = P(T̃c|T0,θ) = nc ∏ i=1 fTc(Tc,i|T0,θ). (4.6) The joint likelihood function for both the ramp load test and the constant load test is given by Lb(θ |T̃s, T̃c,T0) = P(T̃s, T̃c,T0|θ) = P(T̃s|θ)P(T0|T̃s,θ)P(T̃c|T0, T̃s,θ) = P(T̃s|θ)P(T̃c|T0,θ) = ns ∏ i=1 fTs(Ts,i|θ) nc ∏ j=1 fTc(Tc, j|T0,θ). The maximum likelihood estimates can be obtained by maximizing the likeli- hood functions Ls(θ |T̃s), Lc(θ |T̃c,T0) and Lb(θ |T̃s, T̃c,T0). Because the likelihood functions are complicated, the quasi-Newton method can be applied to calculate the maximum likelihood estimates. 4.3.1 The Likelihood Functions In this section, we derive the likelihood functions Ls(θ |T̃s), Lc(θ |T̃c,T0) and Lb(θ |T̃s, T̃c,T0). In order to calculate the density functions of Ts and Tc, we need to calculate the density functions of X and Y denoted fX and fY respectively. The density function of Ts can be written as fTs(t|θ) = fX(t|θ), and the density function of Tc can be written as fTc(t|θ) = fX(t|θ)I(t ≤ T0)+ fY (t|θ)I(t > T0), 56 4.3. LIKELIHOODMETHOD where I is the indicator function. To calculate the density function of X , we define the variable transformation as X = exp(a)b exp(b)−1 , V1 = b. Then, a an b can be solved bya = log X{exp(V1)−1} V1 , b =V1. The Jacobian matrix Js of this variable transformation is given by Js = ( ∂a ∂X ∂a ∂V1 ∂b ∂X ∂b ∂V1 ) = ( 1 X exp(V1) exp(V1)−1 − 1 V1 0 1 ) The joint density function of X and V1 is fX ,V1(x,v) = fa,b ( log x{exp(v)−1} v ,v ) |Js| = fa ( log x{exp(v)−1} v ) fb(v) 1 |x| . (4.7) The density function of X can be integrated from (4.7): fX(x) = ∫ fX ,V1(x,v)dv. To calculate the density function of Y , we define the variable transformation as Y = T0+ exp(a) exp(b)−1 [ exp { b−T0 exp(b)−1exp(a) } −1 ] V2 = exp(a) exp(b)−1 . 57 4.3. LIKELIHOODMETHOD Then a and b can be solved by a = log { (Y −T0+V2)exp ( T0 V2 ) −V2 } b = log ( Y −T0 V2 +1 ) + T0 V2 . The Jacobian matrix Jc of this variable transformation is given by Jc = ( ∂a ∂Y ∂a ∂V2 ∂b ∂Y ∂b ∂V2 ) = exp(T0/V2)(Y−T0+V2)exp(T0/V2)−V2 exp(T0/V2)−(T0/V2)exp(T0/V2)−1(Y−T0+V2)exp(T0/V2)−V2 1 Y−T0+V2 −Y+T0 (Y−T0)V2+V 22 − T0 V 22 The joint density function of Y and V2 is fY,V2(y,v) = fa,b ( log { (y−T0+ v)exp ( T0 v ) − v } , log ( y−T0 v +1 ) + T0 v ) |Jc| = fa ( log { (y−T0+ v)exp ( T0 v ) − v }) × fb ( log ( y−T0 v +1 ) + T0 v ) |Jc| (4.8) The density function of Y can be integrated from (4.8): fY (y) = ∫ fY,V2(y,v)dv. We can also calculate the density functions of the logarithm of the breaking times Ts and Tc and estimate the parameters in the logarithm scale of the breaking times. The density functions are simpler in the logarithm scale of the breaking times. 4.3.2 Optimization Method We use the quasi-Newton method to optimize the likelihood functions. We choose a starting point for θ as θ (0), and then update θ (0) to θ (1) by calculating the gra- 58 4.4. APPROXIMATIONS dient of the likelihood functions with respect to θ at the value θ (0). This recursive optimization routine is performed until convergence. We use nlm in R for opti- mization. For the choice of the starting point θ (0), we use both the true values and ran- dom values in the simulation studies. The simulation results are discussed in the next chapter. For real data analysis, we do not have a good starting point θ (0). So we use random starting points in the analysis of data from Foschi and Barrett’s experiments in Section 6.2. Although the maximum likelihood estimates seem promising in theory, it does not work well in practice because the integration limits are hard to determine. To solve this problem, we propose some approximations to the solutions for the breaking times in the U.S. model, in the next section. 4.4 Approximations In this section, we discuss three approximations of the solutions for the breaking times in the U.S. model. The approximate solutions yield simple likelihoods. In addition, the approximations are required for the quantile method we propose in the next section. Details of the approximations are given in Sections 4.4.1, 4.4.2 and 4.4.3. We first approximate exp(b)− 1 by exp(b) in the expressions of X and Y , 59 4.4. APPROXIMATIONS yielding X = exp(a)b exp(b)−1 1≈ exp(a)b exp(b) = bexp(a−b), Y = T0+ exp(a) exp(b)−1 { exp ( b−T0 exp(b)−1exp(a) ) −1 } 1≈ T0+ exp(a)exp(b) { exp ( b−T0 exp(b)exp(a) ) −1 } , which we approximate further as Y 2≈ T0+ exp(a)exp(b) exp { b−T0 exp(b)exp(a) } = T0+ exp{a−T0 exp(b−a)}. We then approximate the logarithm of X and the logarithm of Y−T0 via Taylor series. For X , we expand log(b) about µb = E(b), and for Y − T0, we expand exp(a−b) about µb−µa = E(b)−E(a), yielding log(X) ≈ a−b+ log(b) 3≈ a−b+ log(µb)+ 1µb (b−µb), (4.9) log(Y −T0) ≈ a−T0 exp(b−a) 3≈ a−T0 exp(µb−µa)(b−a−µb+µa+1). (4.10) The approximations of the solutions for Ts and Tc can be written in terms of the approximations of X and Y . 60 4.4. APPROXIMATIONS For Ts, the first approximation yields Ts 1≈ Ts,approx1 ≡ bexp(a−b), (4.11) and according to (4.9), the Taylor series expansion of Ts after the first approxima- tion yields log(Ts) 3≈ log(Ts,approx3)≡ a−b+ log(µb)+ 1µb (b−µb). (4.12) For Tc, the first approximation yields Tc 1≈Tc,approx1≡ bexp(a−b), if exp(a)b exp(b)−1 ≤ λT0, T0+ exp(a−b) [exp{b−T0 exp(b−a)}−1] , if exp(a)bexp(b)−1 > λT0. (4.13) The second approximation yields Tc 2≈ Tc,approx2 ≡ bexp(a−b), if exp(a)b exp(b)−1 ≤ λT0, T0+ exp{a−T0 exp(b−a)}, if exp(a)bexp(b)−1 > λT0, (4.14) and according to (4.9) and (4.10), the Taylor series expansions after the first two approximations yields log(Tc,approx3)≡ a−b+ log(µb)+ 1µb (b−µb), if exp(a)b exp(b)−1 ≤ λT0, log(Tc,approx3−T0)≡ a−T0 exp(µb−µa)(b−a−µb+µa+1), if exp(a)bexp(b)−1 > λT0. (4.15) 61 4.4. APPROXIMATIONS We study the accuracy of our approximations by comparing Ts simulated from the original expressions (4.3) with Ts simulated from the approximations (4.11) and (4.12), and by comparing Tc simulated from the original expressions (4.4) with Tc simulated from the approximations (4.13), (4.14) and (4.15). 4.4.1 Approximation Step 1 The first approximation step is approximating exp(b)− 1 by exp(b) in the solu- tions for the breaking times Ts and Tc. It is easy to see that if b is large, 1 is negligible compared to exp(b). We can write the ratio of X and the approximate X as X Xapprox1 = exp(b) exp(b)−1 , which is approximately 1 when b is large. As discussed in Section 3.2, b is positive and the value of b should not be too small in order to generate reasonable breaking times. We conduct a simulation study to investigate the accuracy of the approxima- tion step 1 when µb = 50. The approximation error is negligible for both Ts and Tc when σb = 0.4 or σb = 5. 4.4.2 Approximation Step 2 The approximation step 2 is to approximate exp{b−T0 exp(b−a)}−1 by exp{b−T0 exp(b−a)} in the solution of Tc after the first approximation. In Section 3.2, we noticed that in order to keep the mean breaking times in the ramp load test around one minute, µb should be bigger than µa and the difference should be between 6 and 8. Since T0 is around 0.01, then T0 exp(µb− µa) < 30. We can show that b−T0 exp(b− a) > 10 for most pairs of the random variables (a,b). Thus 1 is negligible compared to exp{b−T0 exp(b−a)}. We perform a simulation study to investigate the accuracy of this approxi- mation. In all simulation studies in Section 4.4.2, we assume that a and b are 62 4.4. APPROXIMATIONS lll l llll lll ll lllll ll ll l l l l lll lll l l l l l l l ll ll l l l l l l llll l l ll l l l l lll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l ll l ll 0 10 20 30 40 0. 00 0 0. 01 0 0. 02 0 σa= 0.4 and σb= 0.4 log(Tc) a pp ro xi m at io n er ro r lll l lll ll ll l l llll lll lllll l ll ll lll ll l l ll l ll ll lll l l l l ll ll l l ll lll l l l l l l l −20 −10 0 10 20 30 40 0. 00 0 0. 00 5 0. 01 0 0. 01 5 σa= 1 and σb= 5 log(Tc) a pp ro xi m at io n er ro r ll ll lll lll l lll lll ll ll l lll ll ll ll llll ll l ll l lll l ll l l l l ll l l l l ll l lllll l l ll l −20 0 20 40 0. 00 0 0. 00 5 0. 01 0 0. 01 5 0. 02 0 σa= 5 and σb= 1 log(Tc) a pp ro xi m at io n er ro r ll ll l llll ll lll llll llll l lll l ll ll l ll llll l ll l l l l l l l lll lll lll ll ll l ll ll l l −20 0 20 40 0. 00 0 0. 00 5 0. 01 0 0. 01 5 σa= 5 and σb= 5 log(Tc) a pp ro xi m at io n er ro r Figure 4.1: Plots of the approximation error versus the logarithm of the breaking time Tc when a and b are both Normal, µa = 42 and µb = 50, and p = 0.2 in the constant load test. The vertical line is log(T0). The values for σa and σb are shown in the title. Tc’s are the break- ing times generated from the U.S. model without approximations and Tc,approx2’s are the approximated breaking times generated from the U.S. model with the approximation step 1 and 2. The approximation error is the difference of log(Tc,approx2) and log(Tc). 63 4.4. APPROXIMATIONS both Normal, µa = 42 and µb = 50, and p = 0.2 in the constant load test. We set (σa,σb) = (0.4,0.4), (1,5), (5,1) and (5,5). We simulate 1000 Tc’s from the U.S. model with and without approximation steps 1 and 2, and then compare the differences of the Tc and the approximate Tc. Figure 4.1 contains plots of the approximation error versus the logarithm of the breaking time Tc. The logarithm of T0 is indicated as the vertical line. From Figure 4.1, we notice that the approximation step 2 is very accurate for most Tc’s except for some values in the lower middle. We recall that the approx- imation step 2 is only for the random variable Y , and Tc can be written in terms of X and Y as in (4.4). For Tc ≤ T0, we have that Tc = X and the approximation step 2 is not used, so the only error is caused by the approximation step 1 and this error is small. For Tc > T0, if we have that Tc−T0 is close to 0, the approximation step 2 is not very accurate. However, the magnitude of the approximation error is still small. For Tc > T0 and Tc much greater than T0, the approximation is very accurate. From Figure 4.1, we see that, the smaller σa and σb, the more inaccuracies in the approximation when Tc > T0 but close to T0. This is consistent with the results from Section 3.2.4, where showed that we would expect to see more boards break during the early stage of the constant load test when σa and σb are small. In conclusion, the approximation step 2 is very accurate for most Tc’s except for those in the starting stage of the constant loading part of the constant load test. However, the magnitude of the difference caused by the approximation is still small. 4.4.3 Approximation Step 3: Linearization In this section, we discuss the linearization of the approximate random variables Ts and Tc in the logarithm scale after the first two approximation steps. We can also linearize the random variable Ts and Tc without the first two approximation steps. However, the results are too complicated and not useful. The linearization of the logarithm of the approximated random variable X after 64 4.4. APPROXIMATIONS the approximation step 1 is given by the Taylor series about µb = E(b): log(X) ≈ a−b+ log(b) ≈ a−b+ log(µb)+ 1µb (b−µb). (4.16) Using the Lagrange form of the remainder of the Taylor series, the approxi- mation error caused by the Taylor series is equal to (ξ −µb)2 2µ2b , where ξ is between b and µb. If we assume that b is Normal, by the empirical rule, it can be shown that P { (ξ −µb)2 2µ2b > 9σ2b 2µ2b } ≤ P { (ξ −µb)2 2µ2b > 9σ2b 2µ2b } = P { (ξ −µb)2 > 9σ2b } ≤ P(|b−µb|> 3σb) < 0.01. The last inequality is from the empirical rule (or so-called three-σ rule). We choose the value 9σ2b/2µ 2 b to use the empirical rule. The above result means when σb is small and µb is large, the approximation error caused by linearization is negligible with a probability close to 1. We can write the approximation error caused by the Taylor series expansion of the logarithm of X as log(X)− log(Xapprox) = log ( b µb ) + 1 µb (b−µb), (4.17) 65 4.4. APPROXIMATIONS which also shows that the linearization is more accurate when b is closer to µb. The linearization of the logarithm of the random variable Y −T0 after the ap- proximation steps 1 and 2 is given by the Taylor series expansion of exp(b− a) about µb−µa: log(Y −T0) ≈ a−T0 exp(b−a) ≈ a−T0 exp(µb−µa)(b−a−µb+µa+1). (4.18) Using the Lagrange form of the remainder of the Taylor series, the approxi- mation error caused by the Taylor series expansion is equal to 1 2 T0 exp(µb−µa){ξ − (µb−µa)}2, where ξ is between b−a and µb−µa. However, unlike the linearization of X , the approximation error caused by the linearization of Y −T0 is only negligible when b−a is very close to µb−µa. We perform a simulation study to investigate the accuracy of the linearization. The simulation setups are the same as in Section 4.4.2. The approximate breaking times Ts’s and Tc’s are calculated after all three approximations. Figure 4.2 plots of the approximation error versus the logarithm of the break- ing time Ts. Figure 4.3 plots of the approximation error versus the logarithm of the breaking time Tc. Figure 4.2 suggests the Taylor series expansion for Ts is reasonably accurate. In all panels, the magnitude of the approximation errors is relatively small com- pared to the magnitude of Ts for all panels. In the upper left panel of Figure 4.2, the magnitude of the approximation errors is much smaller than in the other pan- els, which means the Taylor series expansion is more accurate when σb is small, as shown in (4.17). However, in the other three panels, the magnitude of the ap- proximation errors is still relatively small compared to the magnitude of Ts’s for most values. In the upper right panel of Figure 4.2, we see a pattern of the change of the approximations errors when Ts increases, which does not appear in the other 66 4.4. APPROXIMATIONS l l l ll l l l l l ll l l ll l l l l l lll l l l l l ll l l l l l ll ll l l l l l l l l l l l l l l l l l l l ll l l l l l ll l l l l ll l l l l l l ll l ll l l ll l l ll l l l l l l l l l l ll l l ll l l ll l l l l l l l l l l l ll l l l l l l l l l l l l l l l l ll l l l l l l l l ll l l l ll ll l l l l l ll l l l l ll l l l l l ll l l l l l l l l l l l ll l l l l ll l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l ll l l l l l l l l l l l l l l l ll l l l l l l l l ll l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l ll l l l l l l ll l l l l l l l l l l l l l ll ll l l l l l l l l l l l l l l l l l ll l ll l l l l l l l l l l l −5.5 −4.5 −3.5 −2.5 0e +0 0 2e −0 4 4e −0 4 σa= 0.4 and σb= 0.4 log(Ts) a pp ro xi m at io n er ro r ll l ll l l l l l ll l l l l l l lll l l l l l l l l l l l l ll l l l l ll l l l l l l l l l l l l l l l l l l l ll ll l l l l l l l l l l l l l l l ll l l l l l l l l l ll l l l l l l l l l l l l l l l l l ll l l l ll l l l ll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l lll l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l ll l l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l ll l −20 −10 −5 0 5 10 0. 00 0. 02 0. 04 0. 06 σa= 1 and σb= 5 log(Ts) a pp ro xi m at io n er ro r l l l ll l l l l l ll l l lll l l l l lll l l l l l l ll l ll l l ll ll l l l ll l l l l l l l l l l l l l l l l ll l l l ll l ll ll l ll l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l ll l l l ll l l ll l ll l ll l l ll l l l ll l l l l l l l l l l l ll l l l l l l l l l l l l l l l l l ll l l ll ll l l l l l l l l ll l lll l l l l l l l l l l l l l l l l ll ll ll l l ll l l l l l l ll l l l l l l l l l l l l l l l ll l l l l l l l l l l l ll l l l l l l l l l ll l l l l l l l l l ll ll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l ll l l ll ll l l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l ll l l l ll l l l l l ll l l ll ll l l ll ll l ll l l l l l l l l l l −20 −10 −5 0 5 10 0. 00 00 0. 00 10 0. 00 20 σa= 5 and σb= 1 log(Ts) a pp ro xi m at io n er ro r l l l ll l l l l l ll l l ll l l l l ll ll l l l l l lll l l l l ll ll l l l l l l l l l l l l l l l l l l l ll l l l l l ll l l l l ll l l l l l ll l ll l l ll l l ll l l l l l l l l l l l l l l ll l ll ll l l l l l l l l l l l ll l l l l l l l l l l l l l l l l l ll l l l l l l l l ll l l l ll ll l l l l ll l l l l l ll l l l l l ll l l l l l l l l l ll l l l l ll l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l ll l l l l l l l l ll l l l l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l l ll l ll l l l l l l ll l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l ll l l l ll l ll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l −20 −10 0 10 20 0. 00 0. 02 0. 04 0. 06 σa= 5 and σb= 5 log(Ts) a pp ro xi m at io n er ro r Figure 4.2: Plots of the approximation error versus the logarithm of the breaking time Ts when a and b are both Normal, µa = 42 and µb = 50. The values for σa and σb are shown in the title. Ts’s are the break- ing times generated from the U.S. model without approximations and Ts,approx3’s are the approximated breaking times generated from the U.S. model with the approximation steps 1, 2 and 3. The approxima- tion error is the difference of log(Ts,approx3) and log(Ts). 67 4.4. APPROXIMATIONS lll l ll l l l l l l l l l ll l l l l l l l l ll l l l ll lll ll ll l ll l l l l l l l l ll l l l l l l l l lll l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l 0 10 20 30 40 0 10 20 30 40 σa= 0.4 and σb= 0.4 log(Tc) a pp ro xi m at io n er ro r lll l lll ll ll l l l l l l l l l lll l ll l lll ll ll l ll l l l l l l l ll l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l ll l l l l l ll l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l l l l −20 −10 0 10 20 30 40 0 10 20 30 40 σa= 1 and σb= 5 log(Tc) a pp ro xi m at io n er ro r ll ll lll lll l l ll ll l l l l ll l l l l ll l l l ll l ll l l l l l l l l l l l l l ll l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l ll l l l l l l l l l l l l l l l l l l −20 0 20 40 0 10 20 30 40 σa= 5 and σb= 1 log(Tc) a pp ro xi m at io n er ro r ll ll l llll ll lll lll l llll l lll l l l l lll l l l l l l l l l l l l l l lll l l l ll l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l l l l ll l l l l l ll l l l l l l l l l ll l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l −20 0 20 40 0 10 20 30 40 50 σa= 5 and σb= 5 log(Tc) a pp ro xi m at io n er ro r Figure 4.3: Plots of the approximation error versus the logarithm of the breaking time Tc when a and b are both Normal, µa = 42 and µb = 50, and p = 0.2 in the constant load test. The values for σa and σb are shown in the title. Tc’s are the breaking times generated from the U.S. model without approximations and Tc,approx3’s are the approximated breaking times generated from the U.S. model with the approximation steps 1, 2 and 3.The vertical line denotes the logarithm of the median of the breaking times Tc’s. The approximation error is the difference of log(Tc,approx3) and log(Tc). 68 4.5. QUANTILE METHOD three panels. From Figure 4.3, the Taylor series expansion for Tc is accurate for the values which are close to the median of Tc (indicated by the vertical line). As we can infer from (4.18), the Taylor series expansion of a random viable with the expo- nential form is only accurate for the values very close to the fixed point chosen for the Taylor’s expansion. From the upper left panel of Figure 4.3, the approx- imation errors increases slowly when Tc departs from the median. However, the approximation errors increases rapidly when Tc departs from the median in the other three panels. In conclusion, the linearization works acceptably well for Ts’s, but only works well for those Tc’s close to the median of Tc. 4.4.4 Approximate Likelihood Functions The likelihood functions are very complicated for the original solutions for Ts and Tc in the U.S. model as shown in Section 4.3, which makes difficult, the integration and optimization process involved. In contrast, the likelihood functions for the approximate solutions for Ts and Tc are easier to calculate. The results in Section 4.4.1 and 4.4.2 show the approximation steps 1 and 2 to be very accurate for most values of Ts and Tc. So to calculate the approximate likelihood functions, we use the approximate solution for Ts in (4.12) and the approximate solution for Tc in (4.15) after the approximation steps 1 and 2. The simulation results to evaluate the performance of these approximate maximum likelihood estimates are discussed in Chapter 5. 4.5 Quantile Method In this section, we propose a method for parameter estimation based on the mo- ment estimates and the quantile estimates of the distribution of the breaking times. We call this method the quantile method. The quantile method is based on our ap- proximations to Ts and Tc. 69 4.5. QUANTILE METHOD We discuss the quantile method under the constraint that b is positive. In the revised U.S. model, b is defined as positive, which means the Normal assumption on b is not appropriate. However, if the standard deviation σb is small relative to the expected value of b, the probability of b is non-positive is negligible. From (4.16) and (4.18), log(X) and log(Y −T0) can be written as linear com- binations of the random variables a and b, i.e., log(X) ≈ a− (1−1/µb)b+ log(µb)−1, log(Y −T0) ≈ a+ c0a− c0b− c0(µa−µb+1), where c0 = T0 exp(µb−µa). If a and b are independent, then E{log(X)}= µa−µb+ log(µb), var{log(X)}= σ2a +(1−1/µb)2σ2b , E{log(Y −T0)}= µa− c0, var{log(Y −T0)}= (1+ c0)2σ2a + c20σ2b . We denote the means of log(X) and log(Y − T0) as µX and µY respectively, and the variances of log(X) and log(Y −T0) as σ2X and σ2Y respectively. We can estimate µa and µb from the estimates of µX and µY . The equations{ µX = µa−µb+ log(µb), µY = µa−T0 exp(µb−µa), (4.19) yield µX −µY +µb− log(µb)−T0µb exp(−µX) = 0. (4.20) We can numerically solve for µb from (4.20) given the estimates of µX and µY , and then solve for µa from (4.19). 70 4.5. QUANTILE METHOD We can estimate σ2a and σ2b from the estimates of σ 2 X and σ2Y . The equations{ σ2X = σ 2 a +(1−1/µb)2σ2b , σ2Y = (1+ c0) 2σ2a + c 2 0σ 2 b , (4.21) yield{ σ2a = σ 2 X − (1−1/µb)2[(1+ c0)2σ2X −σ2Y ]/{(1+ c0)2(1−1/µb)2− c20}, σ2b = {(1+ c0)2σ2X −σ2Y}/{(1+ c0)2(1−1/µb)2− c20}. (4.22) We can solve for σ2a and σ2b in (4.22) given µa, µb, σ 2 X and σ2Y . We propose the quantile method to estimate µX , µY , σ2X and σ2Y in this section. To estimate µX and µY , we do not need the Normality assumptions on a and b. To estimate σX and σY , we need the Normality assumptions on a and b. The estimation errors of the quantile method consist of the random error and the approximation error. The random error can be reduced by increasing the sam- ple size, and the effect of the approximation error can be reduced by choosing reasonable quantiles. In Section 4.5.1 and Section 4.5.2, we discuss the general quantile method to estimate µX and µY , and σX and σY without considering the approximation error. In Section 4.5.3, we discuss the method to reduce the effect of the approximation error. 4.5.1 Estimates of the Means In this section, we discuss the quantile method to estimate µX and µY from the breaking times Ts and Tc. Let mX , mY , mTs and mTc be the median of the distributions of X , Y , Ts and Tc respectively. 71 4.5. QUANTILE METHOD From the ramp load test, we observe Ts. Ts = X , so we can simply estimate µX by the sample mean of the logarithm of the measured Ts’s. We can also estimate µX by the median of the logarithm of the measured Ts’s if log(Ts) is symmetrically distributed. From the constant load test, we observe the Tc’s, which are equal to X or Y according as X ≤ T0 or X > T0 since by (4.4): Tc = { X , if X ≤ T0, Y, if X > T0, where T0 is set to be the p (p < 0.5) percentile of Ts’s from the ramp load test. However, we cannot estimate µY by the sample mean of the logarithm of the observed Tc−T0 since E{log(Tc−T0)} 6= µY . (4.23) Instead, when log(Y −T0) is symmetrically distributed, we can estimate µY by the logarithm of the median of the Tc−T0’s provided we can show that mTc = mY = exp(µY )+T0. (4.24) Theorem 1. The median of Tc equals to the median of Y when b is positive and p < 0.5. Proof. To show that mTc = mY , we calculate P(Tc ≥ mY ): P(Tc ≥ mY ) = P(Tc ≥ mY ,X ≤ T0)+P(Tc ≥ mY ,X > T0) = P(X ≥ mY ,X ≤ T0)+P(Y ≥ mY ,X > T0) = P(mY ≤ X ≤ T0)+P(Y ≥ mY )P(X > T0|Y ≥ mY ) = P(mY ≤ X ≤ T0)+0.5P(X > T0|Y ≥ mY ). (4.25) First, we notice that T0 is the p (p < 0.5) percentile of X , so the median of X , 72 4.5. QUANTILE METHOD which is the 50-th percentile of X , is larger than T0, i.e., T0 < mX . Second, we notice that if X > T0, then Y > X . (4.26) because according to (4.2), Y = T0+ X b { exp ( b−bT0 X ) −1 } > T0+ X b { 1+b−bT0 X −1 } = T0+X−T0 = X . As a result, the median of X is smaller than the median of Y , i.e., T0 < mX < mY , (4.27) which yields P(mY ≤ X ≤ T0) = 0 in the first part of equation (4.25). Third, we notice that the definition of Y in (4.2) yields X > T0⇔ Y > T0. (4.28) 73 4.5. QUANTILE METHOD The reason is that, b > 0 and according to (4.2), X > T0 ⇔ b > bT0X ⇔ exp ( b−bT0 X ) > 1 ⇔ X b { exp ( b−bT0 X ) −1 } > 0 ⇔ T0+ Xb { exp ( b−bT0 X ) −1 } > T0 ⇔ Y > T0. Therefore, the second part of equation (4.25) can be calculated as 0.5P(X > T0|Y ≥ mY ) = 0.5P(Y > T0|Y ≥ mY ) = 0.5×1 = 0.5, (4.29) since mY > T0 as shown in (4.27). To sum up, we have shown that under the constraint b > 0, P(Tc ≥ mY ) = 0.5. (4.30) As a corollary of Theorem 1, µY can be estimated by the logarithm of the median of Tc minus T0 when log(Y −T0) is symmetrically distributed. The above proof is not true when b is assumed to be Normal since b can be negative if it is Normal distributed. However, if the standard deviation σb is small, the equation (4.30) can be approximated as P(Tc ≥ mY )≈ 0.5. (4.31) Under the assumption that b is Normally distributed, we demonstrate the plau- sibility of (4.31) by simulation plots. The simulation setups are the same as in Sec- tion 4.4.2. We first simulate 1000 X’s to calculate T0, then we simulate another 74 4.5. QUANTILE METHOD 1000 pairs of (a,b)’s and calculate the X’s, Y ’s and Tc’s for each pair. We paint all the simulated points grey, and then paint those points which satisfy Tc > mY black. The plots are shown in Figure 4.4. l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l ll l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l ll l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l ll l l l l l l l l l l ll l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l lll l ll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l ll l l l l ll ll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l ll l l ll l l l l l l l l l l l l l l l l l ll l l l l l l l ll l l l l l l l ll l l l l l l l l l l l l l l l lll l l l l l l l ll l l ll l l l l l l l l l l l l l l l l 40.5 41.0 41.5 42.0 42.5 43.0 49 .0 50 .0 51 .0 σa= 0.4 and σb= 0.4 a b l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l ll l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l ll l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l ll l l l l l l l l l l ll l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l lll l ll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l ll l l l l ll ll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l ll l l ll l l l l l l l l l l l l l l l l l ll l l l l l l l ll l l l l l l l ll l l l l l l l l l l l l l l l lll l l l l l l l ll l l ll l l l l l l l l l l l l l l l l 39 40 41 42 43 44 45 35 40 45 50 55 60 65 σa= 1 and σb= 5 a b l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l ll l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l ll l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l ll l l l l l l l l l l ll l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l lll l ll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l ll l l l l ll ll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l ll l l ll l l l l l l l l l l l l l l l l l ll l l l l l l l ll l l l l l l l ll l l l l l l l l l l l l l l l lll l l l l l l l ll l l ll l l l l l l l l l l l l l l l l 25 30 35 40 45 50 55 47 48 49 50 51 52 53 σa= 5 and σb= 1 a b l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l ll l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l ll l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l ll l l l l l l l l l l ll l l l l l l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l lll l ll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l ll l l l l ll ll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l ll l l ll l l l l l l l l l l l l l l l l l ll l l l l l l l ll l l l l l l l ll l l l l l l l l l l l l l l l lll l l l l l l l ll l l ll l l l l l l l l l l l l l l l l 25 30 35 40 45 50 55 35 40 45 50 55 60 65 σa= 5 and σb= 5 a b Figure 4.4: Plots for illustrating the quantile method when a and b are both Normal, µa = 42 and µb = 50, and p = 0.2 in the constant load test. The values for σa and σb are shown in the title. We paint all the sim- ulated points grey, and then paint those points which satisfy Tc > mY black. These plots support our claim in equation (4.31). From Figure 4.4, we notice that the proportion of the points which satisfy Tc≥ 75 4.5. QUANTILE METHOD mY is around 0.5, although the border line is not of the same shape for different values of the standard deviations. We also confirm this by counting the number of points which satisfy that Tc ≥ mY , which is 500 for every plot. 4.5.2 Estimates of the Variances In this section, we discuss the quantile method to estimate σX and σY from the breaking times Ts and Tc. The Normality assumptions on a and b are needed for this method. If a and b are both assumed to be Normal, then log(X) and log(Y − T0) are both approximately Normal, i.e., log(X)∼ N(µX ,σ2X), (4.32) where µX = µa−µb+ log(µb) and σ2X = σ2a +(1−1/µb)2σ2b , and log(Y −T0)∼ N(µY ,σ2Y ), (4.33) where c0 = T0 exp(µb−µa), µY = µa− c0 and σ2Y = (1+ c0)2σ2a + c20σ2b . For the ramp load test, we can simply estimate σ2X by the sample variance of the logarithm of the Ts’s. We can also estimate σ2X by the quantiles of the logarithm of the Ts’s. Let t1,s and t2,s be the p1,s quantile and the p2,s quantile of the logarithm of the Ts’s. Because log(Ts) is approximately Normal after the linearization under the Normal assumptions on a and b, the difference of t2,s− t1,s is linked to σX , i.e., t2,s− t1,s = csσX where cs is a constant depending on p1,s and p2,s. For example, if we choose p1,s = 0.1587 and p2,s = 0.8413, then t2,s− t1,s ≈ 2σX according to the empirical rule. For the constant load test, we cannot estimate σ2Y by the variance of the loga- rithm of the (Tc−T0)’s since not all the Tc’s are defined as Y ’s. However, we can 76 4.5. QUANTILE METHOD still use the quantile estimates. Let t1,c and t2,c be the p1,c quantile and the p2,c quantile of the logarithmic of the Tc−T0’s. Let q1,Y and q2,Y be the p1,c quantile and the p2,c quantile of the logarithm of Y −T0’s. If we can show that t1,c = q1,Y and t2,c = q2,Y , then t2,c− t1,c = q2,Y −q1,Y = ccσY (4.34) where cc is a constant depending on p1,c and p2,c. We can prove that t1,c = q1,Y and t2,c = q2,Y using the same method as in the proof of the median estimates in Section 4.5.1 under the conditions that b > 0, p1,c > p and p2,c > p. Since the logarithm transformation and the shift transformation do not affect the quantiles, we only need to show that the p1,c (or p2,c) quantile of Tc equals to the p1,c (or p2,c) quantile of Y when b is positive and p1,c > p (or p2,c > p). Theorem 2. The p1,c-th (or p2,c-th) quantile of Tc equals the p1,c-th (or p2,c-th) quantile of Y , when b is positive and p1,c > p (or p2,c > p). The proof of the above Theorem 2 is similar to the proof of Theorem 1. We need only change mY in the proof of Theorem 1 to q1,Y (or q2,Y ), and change 0.5 in the proof of Theorem 1 to p1,c (or p2,c) to get a proof of Theorem 2. Notice that in the above proof, we require p1,c > p and p2,c > p. Therefore, we cannot choose p1,c too small. For example, if p = 0.2, we can choose 0.3 ≤ p1,c ≤ 0.5 and p2,c larger than 0.5. 4.5.3 Approximation Errors As noted in Section 4.4, the linearization is not accurate for all Tc’s, and the farther Tc is from the median, the less accurate the linearization (see Figure 4.3). There- fore, we should choose p1,c and p2,c close to 0.5 to reduce the approximation errors. Choosing reasonable values for p1,s, p2,s, p1,c and p2,c seems quite challeng- ing. In the simulation studies, we chose p1,s = 0.2, p2,s = 0.8, p1,c = 0.45 and p2,c = 0.55, and results are discussed in Chapter 5 with more details. 77 4.6. A METHOD THAT COMBINES THE MAXIMUM LIKELIHOOD AND THE QUANTILE METHOD 4.6 A Method that Combines the Maximum Likelihood and the Quantile Method In this section, we propose a method that combines the approximate maximum likelihood and the quantile method for parameter estimation. As noted in Section 4.5, the quantile method works well for quantiles close to the median of the random variables. Thus the estimates of µa and µb may be accurate but the estimates of σa and σb may not due to large approximation errors. As an alternative, we propose a two-step estimation method that combines the quantile method and the maximum likelihood. First, we estimate µa and µb from the quantile method using the estimates of the median of Ts and Tc, as discussed in Section 4.5.1. Second, we estimate σa and σb by maximizing the likelihood functions using the µ̂a and µ̂b from the first step. This method is more computationally efficient than the likelihood method and it improves the estimates of the standard deviations. The simulation results are discussed in Chapter 5. 78 Chapter 5 Simulation Studies 5.1 Introduction This chapter presents the results of simulation studies to investigate the perfor- mance of the parameter estimation methods proposed for the U.S. model in Chap- ter 4. We simulate the breaking times Ts in the ramp load test and Tc in the con- stant load test from the U.S. model, as discussed in Chapter 3, and then apply the parameter estimation methods in Chapter 4 to estimate the model parameter θ = (µa,µb,σa,σb). Section 5.2 contains the simulation setups and results for the approximate maximum likelihood estimates. Section 5.3 describes the simulation setups and results for the quantile method. Section 5.4 describes the simulation setups and results for the method that combines the approximate maximum likelihood and the quantile method. Section 5.5 summarizes the simulation results from three methods. 79 5.2. GENERAL SETUP FOR THE SIMULATION STUDIES 5.2 General Setup for the Simulation Studies In the simulation studies of this chapter, we set θ = (42,50,0.4,0.4) or θ = (42,50,5,5). We use p = 0.2 in the constant load test. We use ns = nc = 200 in the simulation studies for the three parameter es- timation methods proposed in Chapter 4: the approximate maximum likelihood estimates, the quantile method, and the method that combines the approximate maximum likelihood and the quantile method. The general simulation process is described below. First, for each simulation run, we generate ns = 200 breaking times the Ts’s in the ramp load test and another nc = 200 Tc’s in the constant load test from the U.S. model, as described in Section 3.2. This part is the same for all simulation runs. Second, we estimate the parameter θ using the methods proposed in Chapter 5. This part is different for different parameter estimation methods. Details are discussed in the following sections for each parameter estimation method. We repeat the above simulation process nsim = 100 times. 5.3 Simulation Studies for the Approximate Maximum Likelihood Estimates 5.3.1 Simulation Setups In this section, the parameter setups are the same as that described in Section 5.2. The approximate likelihood functions are calculated after the approximation steps 1 and 2 for the breaking time Ts in the ramp load test and the breaking time Tc in the constant load test, as discussed in Section 4.4.4. To estimate the parameters, we obtain the approximate maximum likelihood estimates numerically from the log-likelihood functions, using the function nlm in R. We set the starting point of θ in the optimization process to be the true value of θ . For the optimization process, the maximum iteration step is set to be 20 and 80 5.3. SIMULATION STUDIES FOR THE APPROXIMATE MAXIMUM LIKELIHOOD ESTIMATES the iteration limit is set to be 200. We record the outputs of the optimization re- sults, which indicates why the optimization process terminates. The output of the optimization results is a number from 1 to 5. Codes 1 and 2 mean the optimization converges, Code 3 means the optimization may not converge, and Codes 4 and 5 mean the optimization does not converge. 5.3.2 Results Figure 5.1 shows the bar plots of the outputs of the optimization results. Figure 5.1 shows that, in quite a few simulation runs, the output is not reliable, since many runs lead to code 3 (may not converge). The codes 4 and 5 do not ap- pear in the 100 simulation runs. The approximate maximum likelihood estimates do not work well in this simulation studies. We will explain possible reasons in next section. Ts only Tc only Ts and Tc 1 : Converge 2 : Converge 3 : May not Converge 0 20 40 60 80 10 0 Figure 5.1: Bar plots of the optimization codes for calculating the approx- imate maximum likelihoodestimates in nsim = 100 simulation runs, when θ = (µa,µb,σa,σb) = (42,50,0.4,0.4), ns = nc = 200 and p = 0.2 in the constant load test. 81 5.3. SIMULATION STUDIES FOR THE APPROXIMATE MAXIMUM LIKELIHOOD ESTIMATES Similar results can be shown when θ = (42,50,5,5). The output of the codes for the approximate maximum likelihood estimates work less well when σa and σb are large. There we see even fewer convergences. 5.3.3 The Shapes of the Likelihood Functions From the results in the previous section, we believe that the reason for the bad performance of the nlm may be that the likelihood functions are flat, or even that the approximate maximum likelihood estimate may not be unique. We perform a simulation study to investigate this idea. More precisely, we generate ns = 200 Ts’s and nc = 200 Tc’s from the U.S. model when θ = (42,50,0.4,0.4) and p = 0.2 in the constant load test. We cal- culate the log-likelihood functions for the Ts’s in the ramp load test, for the Tc’s in the constant load test, and for both the Ts’s and Tc’s, where σa = 0.4 and σb = 0.4 are fixed, while µa varies from 40 to 45 and µb varies from 47 to 52. The corre- sponding plots of the log-likelihood functions are shown in Figure 5.2. Figure 5.2 shows that the log-likelihoods do not have an obvious maximum point when σa and σb are fixed. The log-likelihoods are roughly maximized along the line µb−µa ≈ 8. As a result, we may not be able to estimate µa and µb well, but may be able to estimate µb−µa. Figure 5.2 is not a proof for the statement that θ can not be estimated using the approximate maximum likelihood estimates from the U.S. model. However, it provides a possible explanation for the bad performance of the approximate maximum likelihood estimates in the previous section. In conclusion, the approximate maximum likelihood estimates are not reliable for estimating the parameter θ . 82 5.3. SIMULATION STUDIES FOR THE APPROXIMATE MAXIMUM LIKELIHOOD ESTIMATES E(a) 40 41 42 43 44 45 E(b ) 46 47 48 49 50 51 log−likelihood −400 −200 0 200 400 Ts only E(a) 40 41 42 43 44 45 E(b ) 46 47 48 49 50 51 log−likelihood −4400 −4200 −4000 −3800 −3600 Tc only E(a) 40 41 42 43 44 45 E(b ) 46 47 48 49 50 51 log−likelihood −4500 −4000 −3500 −3000 Ts and Tc Figure 5.2: Perspective plots of the log-likelihoods when σa = 0.4 and σb = 0.4 are fixed, and µa varies from 40 to 45 and µb varies from 47 to 52. 83 5.4. SIMULATION STUDIES FOR THE QUANTILE METHOD 5.4 Simulation Studies for the Quantile Method 5.4.1 Simulation Setups The simulation setups are the same as that discussed in Section 5.2 except that we also consider the case ns = nc = 200,000 for illustrative purposes. For ns = nc = 200, the quantile method estimates µa and µb well, but it pro- duces negative estimates of σa and σb in most simulation runs (see (4.22)). In- creasing the sample size to ns = nc = 200,000 diminishes but does not elimi- nate the problem. The reason is still unclear at this moment. To illustrate the asymptotic performance of the quantile method for estimating σa and σb, we use ns = nc = 200,000 although those numbers are not realistic in practice. For estimation, as in Section 4.5.1, we estimate µa and µb using the estimates of the medians of Ts and Tc, and then as in Section 4.5.2, we estimate σa and σb using the estimates of the quantiles of Ts and Tc. The estimation process is the same as in Section 4.5. We can also estimate µX and σ2X by the sample mean and sample variance of the Ts’s instead of using the median and the quantiles. The results are similar from the two approaches. Here, we only show the simulation results using the median and the quantiles of Ts. The quantiles used in this simulation studies are: p1,s = 0.2, p2,s = 0.8, p1,c = 0.45 and p2,c = 0.55, as discussed in Section 4.5.3. 5.4.2 Results The simulation results using the quantile method are summarized in Figure 5.3 and Figure 5.4. Figure 5.3 contains the box plots for µ̂a and µ̂b for both ns = nc = 200 and ns = nc = 200,000, and the box plots for σ̂a and σ̂b for ns = nc = 200,000 only (when θ = (42,50,0.4,0.4) and p = 0.2 in the constant load test). We do not show the box plots for σ̂a and σ̂b when ns = nc = 200 since the quantile method produces positive estimates for both σ2a and σ2b only in 3 simulation runs out of 84 5.4. SIMULATION STUDIES FOR THE QUANTILE METHOD 100 runs. Figure 5.4 contains the box plots for µ̂a, µ̂b, σ̂a and σ̂b when ns = nc = 200, θ = (42,50,5,5) and p = 0.2 in the constant load test. Figure 5.5 contains the box plots for µ̂a, µ̂b, σ̂a and σ̂b when ns = nc = 200,000, θ = (42,50,5,5) and p = 0.2 in the constant load test. From Figure 5.3, the quantile method estimates µa and µb well when ns = nc = 200, and better when ns = nc = 200,000. The variabilities of the estimates are smaller when ns and nc are larger. This can be explained by the fact that, when ns and nc are larger, the median estimates for Ts and Tc are more accurate, so the estimates for µa and µb are more accurate. Similar results can be shown if we use the mean estimate of Ts instead of the median estimate of Ts. From Figure 5.3, the quantile method estimates σa and σb well when ns = nc = 200,000. The quantile method gives positive estimates for both σ2a and σ2b successfully in 54 simulations runs out of 100 runs, and the box plots show that the estimates are reasonably accurate. From Figure 5.4 and Figure 5.5, the estimates of µa and µb are biased using the quantile method when σa = σb = 5. The standard errors of the estimates are smaller when ns and nc are larger, but the estimates of µa and µb are biased in both figures. This can be explained by the fact that, when σa and σb are large, the approximations 1, 2 and 3 discussed in Section 4.4 are not accurate anymore, as shown in Figure 4.2 and Figure 4.3. As a result, although the estimates of the medians of Ts and Tc are accurate when ns and nc are large, the equations that link µa and µb to the median of log(Ts) and log(Tc−T0) do not hold anymore. As a consequence, the estimates for µa and µb are not accurate. From Figure 5.4, the estimates of σa and σb are still acceptable although the estimates of µa and µb are not very accurate from the quantile method. From Figure 5.5, the estimates of σa and σb are biased from the quantile method. The quantile method gives positive estimates for both σa and σb in 57 simulations runs out of 100 runs. This can also be explained by the fact that, when σa and σb are large, the approximation steps 1, 2 and 3 used for the quantile 85 5.4. SIMULATION STUDIES FOR THE QUANTILE METHOD ll25 35 45 55 ns = nc = 200 µa ll 35 45 55 65 ns = nc = 200 µb 41 .6 42 .0 42 .4 ns = nc = 200,000 µa 49 .6 50 .0 50 .4 ns = nc = 200,000 µb 0. 1 0. 2 0. 3 0. 4 0. 5 ns = nc = 200,000 σa 0. 2 0. 4 0. 6 0. 8 ns = nc = 200,000 σb Figure 5.3: Box lots of the estimates using the quantile method, when θ = (42,50,0.4,0.4) and p = 0.2. The sample size ns and nc used in the simulation runs are shown in plot titles. The grey line indicates the true value of the parameter. 86 5.4. SIMULATION STUDIES FOR THE QUANTILE METHOD l l39 40 41 42 µa l l l 46 47 48 49 50 51 µb 4 5 6 7 σa 1 2 3 4 5 6 7 σb Figure 5.4: Box lots of the estimates using the quantile method, when θ = (µa,µb,σa,σb) = (42,50,5,5) and p= 0.2. The sample size ns = nc = 200. The grey line indicates the true value of the parameter. 87 5.4. SIMULATION STUDIES FOR THE QUANTILE METHOD ll 40 .8 0 40 .8 4 40 .8 8 40 .9 2 µa 48 .7 5 48 .8 0 48 .8 5 48 .9 0 µb l 6. 85 6. 90 6. 95 7. 00 σa 0. 0 0. 5 1. 0 1. 5 σb Figure 5.5: Box lots of the estimates using the quantile method, when θ = (µa,µb,σa,σb) = (42,50,5,5) and p= 0.2. The sample size ns = nc = 200,000. No grey line is shown in the plots because the estimates are too far from the true values. 88 5.5. SIMULATION STUDIES FOR THE COMBINED METHOD method are not accurate. In conclusion, the quantile method estimates µa and µb well when σa and σb are small. The quantile method estimates σa and σb well when σa and σb are small and the sample sizes are extremely large (e.g., 200,000), or when σa and σb are large and the sample sizes are small. Note that a sample of size 200,000 is unreasonable in practice. Here we only use it for illustration. 5.5 Simulation Studies for the Combined Method 5.5.1 Simulation Setups In this simulation study, we use the same setup as that described in Section 5.2. To estimate the parameters, as discussed in Section 4.6, we first estimate µa and µb using the medians of Ts and Tc, and then estimate σa and σb by approximate maximum likelihood assuming the true means are µ̂a and µ̂b. We can also use the sample mean of Ts instead of the median of Ts. The results are similar from the two approaches. Here, we only show the simulation resulting using the median of Ts. To calculate the approximate maximum likelihood estimates of σa and σb, we use the approximate likelihood functions after the approximation steps 1 and 2 for both Ts’s and Tc’s given µ̂a and µ̂b from the quantile estimates, as described in Section 4.6. We choose the starting point (σa,σb) for the optimization process to be the true value, and random values uniformly generated from an area contains the true value for each simulation run. The area is 0 < a < 1 and 0 < b < 1 for (σa,σb) = (0.4,0.4), and 0 < a < 10 and 0 < b < 10 for (σa,σb) = (5,5). 5.5.2 Results The simulation results are summarized in Figure 5.6 to Figure 5.8. Figure 5.6 contains the bar plots of the outputs of the optimization results for 89 5.5. SIMULATION STUDIES FOR THE COMBINED METHOD Fixed Starting Point Random Starting Point 1 : Converge 2 : Converge 3 : May not Converge 4 : Does not Coverge 0 20 40 60 80 10 0 Figure 5.6: Bar plots of the output of the optimization results in calculating the approximate maximum likelihoodestimates for σa and σb in the combined method when θ = (µa,µb,σa,σb) = (42,50,0.4,0.4) and p= 0.2 in the constant load test. Fixed starting point means the starting point in the optimization step of the combined method is set to be the true value (σa,σb) = (0.4,0.4). Random starting point means the starting point in the optimization step of the combined method is set to be different random numbers uniformly generated from the square 0 < σa < 1 and 0 < σb < 1 for different simulation runs. maximizing the approximate likelihood using the combined method when θ = (42,50,0.4,0.4) and p = 0.2 in the constant load test. Comparing Figure 5.6 to Figure 5.1, we notice that the optimization process converges in the combined method for most simulation runs, better than the opti- mization process in the approximate maximum likelihoodmethod in Section 5.1. This means that we are able to estimate σa and σb for most simulation runs given µ̂a and µ̂b from the quantile estimates. Figure 5.7 contains the plots for the estimates of σa and σb from the com- bined method when θ = (42,50,0.4,0.4). In Figure 5.7, the starting point for the optimization process in the second step is chosen to be the true value (σa,σb) = (0.4,0,4), as well as random starting numbers generated uniformly from the square 90 5.5. SIMULATION STUDIES FOR THE COMBINED METHOD 0 < σa < 1 and 0 < σb < 1. The box plots for µ̂a and µ̂b are shown in Figure 5.3. Figure 5.8 contains the plots for the estimates of σa and σb from the combined method when θ = (42,50,5,5). In Figure 5.8, the starting point for the optimiza- tion process in the second step is chosen to be the true value (σa,σb) = (5,5), as well as random starting numbers generated uniformly from the square 0<σa < 10 and 0 < σb < 10. The box plots for µ̂a and µ̂b have been shown in Figure 5.4. From Figure 5.7, the combined method estimates σa and σb well in both the fixed starting point condition and the random starting point condition. The esti- mates are less variable when the starting point is chosen to be the true value. From the two middle panels of Figure 5.7, the combined method still works acceptably well when the starting points are random numbers. However, from the two lower panels of Figure 5.7, the estimates of σa and σb are not the same as those when the starting points are chosen to be the true value for some simulation runs. The start- ing points influence the estimates for some simulation runs, but do not influence the centre of the overall distributions of the estimates much. From Figure 5.8, the combined method estimates σa and σb acceptably well in both the fixed starting point condition and the random starting point condition, al- though the estimate of σb is slightly biased. Combined with Figure 5.4, Figure 5.8 shows that, although the combined method does not estimate µa and µb well in the first step, it still estimates σa and σb acceptably well using the estimates of µa and µb from the first step. From the two lower panels of Figure 5.8, the combined method is not sensitive to the choice of the starting point in the optimization when σa = σb = 5. In conclusion, the combined method estimates the parameters θ well when σa and σb are small. The combined method estimate the parameter σa and σb well even if µ̂a and µ̂b from the first step of the combined method are not very accurate. The combined method is not very sensitive to the choice of the starting point for optimization in the second step. 91 5.5. SIMULATION STUDIES FOR THE COMBINED METHOD l l l l l l l l ll l 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 Fixed Starting Point σa l l l l l 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 Fixed Starting Point σb l l l l l l l 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 Random Starting Point σa l l l l l l l 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 Random Starting Point σb l l l ll l l l l l l ll l l l l l l ll l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l l l l l l 0.0 0.2 0.4 0.6 0.8 1.0 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 σa:fixed starting point σ a :r a n do m s ta rti ng p oi nt l l l l l l l l l l l l l l ll lll l l l l l l l l l ll l l l l l l ll ll l l l l l l l l l l l l l l l l l l 0.0 0.2 0.4 0.6 0.8 1.0 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 σb:fixed starting point σ b: ra n do m s ta rti ng p oi nt Figure 5.7: Plots of the estimates for σa and σb when θ = (42,50,0.4,0.4) and p = 0.2 in the constant load test. Fixed starting point means the starting point in the optimization step of the combined method is set to be the true value (σa,σb) = (0.4,0.4). Random starting point means the starting point in the optimization step of the combined method is set to be different random numbers in the square 0 < σa < 1 and 0 < σb < 1 for different simulation runs. The grey line indicates the true value of the parameter. 92 5.5. SIMULATION STUDIES FOR THE COMBINED METHOD l l l l 3 4 5 6 7 Fixed Starting Point σa l 3 4 5 6 7 Fixed Starting Point σb l l l l l 3 4 5 6 7 Random Starting Point σa l l l 3 4 5 6 7 Random Starting Point σb lll l l l l ll l l l l l l l l 0 2 4 6 8 10 0 2 4 6 8 10 σa:fixed starting point σ a :r a n do m s ta rti ng p oi nt lll l l l l llll ll l ll ll 0 2 4 6 8 10 0 2 4 6 8 10 σb:fixed starting point σ b: ra n do m s ta rti ng p oi nt Figure 5.8: Plots of the estimates for σa and σb when θ = (42,50,5,5) and p= 0.2 in the constant load test. Fixed starting point means the starting point in the optimization step of the combined method is set to be the true value (σa,σb) = (5,5). Random starting point means the starting point in the optimization step of the combined method is set to be different random numbers in the square 0 < σa < 10 and 0 < σb < 10 for different simulation runs. The grey line indicates the true value of the parameter. 93 5.6. SUMMARY OF THE SIMULATION RESULTS 5.6 Summary of the Simulation Results Table 5.1 summarizes main simulation results in this chapter. Methods mean(µ̂a)−µa mean(µ̂b)−µb mean(σ̂a)−σa mean(σ̂b)−σb Success Rates (σa,σb) (s.e.) (s.e.) (s.e.) (s.e.) Approx MLE - - - - 0.51 (0.4, 0.4) Quantile 0.736 0.740 -0.125 0.058 0.03 (0.4, 0.4) (0.652) (0.670) (0.028) (0.028) (1 for µ̂a and µ̂b) Quantile -1.065 -1.101 0.539 -0.913 0.54 (5, 5) (0.078) (0.116) (0.131) (0.192) (1 for µ̂a and µ̂b) Combined 0.736 0.740 0.014 0.016 0.85 (0.4, 0.4) (0.652) (0.670) (0.010) (0.006) (fixed start) 0.736 0.740 0.256 0.450 0.89 (0.652) (0.670) (0.012) (0.006) (random start) Combined -1.065 -1.101 -0.192 0.492 0.95 (5, 5) (0.078) (0.116) (0.059) (0.048) (fixed start) -1.065 -1.101 -0.245 0.476 0.93 (0.078) (0.116) (0.076) (0.053) (random start) Table 5.1: Summary of the simulation results in Chapter 5 when ns = nc = 200. In all simulation runs, (µa,µb) = (42,50) and in the constant load test p = 0.2. The values of (σa,σb) are shown in the table. The success rates denote the proportion of the converged simulation runs for the approximate maximum likelihood estimates and the combined method, and the rates of the simulation runs which produce positive estimates for both σ2a and σ2b for the quantile method. The approximate maximum likelihood estimates are not reliable since the op- timization process in this method only converges in 51 runs out of 100 runs when the likelihood function for both Ts’s and Tc’s are used. We do not report the esti- mates of the parameters because the estimates do not make sense when only half of the results from the simulation runs are reliable. The quantile method estimates µa and µb well but fails to estimate σa and σb in most simulation runs. It produces reasonable estimates for µa and µb in all simulation runs, but only produces positive estimates for σ2a and σ2b in 3 runs out 94 5.6. SUMMARY OF THE SIMULATION RESULTS of 100 runs when σa = σb = 0.4, and only produces positive estimates for σ2a and σ2b in 54 runs out of 100 runs when σa = σb = 5. The combined method works well in estimating µa, µb, σa and σb. In Ta- ble 5.1, the fixed start means the starting points are chosen to be the true values of the parameters, and the random start means the starting points are chosen to be random numbers. Details are contained in Section 5.5. The mean of the estimates in the successful runs are shown in Table 5.1. The medians of the estimates in the successful runs, which are shown in the box plots in the previous sections, are generally closer to the true values than the means of the estimates for most simulation studies. 95 Chapter 6 Experiments and Data Analysis 6.1 Introduction In this section, we will describe and analyze the trial experiment for the duration of load effects conducted at FPInnovations beginning in July of 2010 and also the experiments conducted by Foschi and Barrett in 1980s. Unfortunately, the data from the experiments conducted at FPInnovations is different from the standard ones for the duration of load effects in the literature, since the load does not increase linearly in time t in the ramp load tests. Therefore, we will only study the cumulative distribution functions of the breaking times from the experiments at FPInnovations. Foschi and Yao’s data will be analyzed by a modification of the combined method we proposed in Section 4.6. Section 4.1 contains the discussion and data analysis for the experiments at FPInnovations. Section 4.2 contains the discussion and data analysis for the ex- periments done by Foschi and Barrett. 96 6.2. THE FPINNOVATIONS EXPERIMENT AND DATA ANALYSIS 6.2 The FPInnovations Experiment and Data Analysis 6.2.1 The FPInnovations Experiment The trial experiment for the duration of load project began in July, 2010, at an FPInnovations laboratory in Vancouver. In this experiment, we study the strength of two wood-based products: Douglas-fir lumber (DF) and laminated-veneer lum- ber (LVL) under two moisture content conditions (high and low). In this chapter, we use DF-High, DF-Low, LVL-High and LVL-Low to denote the four resulting groups of lumber. The trial experiments at FPInnovations consist of four stages: the preparation stage, the conditioning stage, the ramp load test stage and the constant load test stage. In the preparation stage, the Douglas-fir lumber test groups are formed by so-called end-matching: a piece of wood is cut across its width, and then one resulting piece is used in the ramp load test while the remaining piece is used in the constant load test group. The laminated-veneer lumber test groups are formed b so-called side-matching: a piece of wood is cut down the length, and then one resulting piece is used in the ramp load test while the remaining piece is used in the constant load test group. The objective of the matching process is to ensure that the groups for the ramp load test and the groups for the constant load test of each of four combinations of product and condition have similar distributions of strength. The size of each final specimen is 38 inches long, 0.875 inches thick, and 1.5 inches wide. The sample size is 20× 2 for each of the four combinations of product and condition (i.e., 20 specimens for the ramp load test and 20 specimens for the constant load test). In the conditioning stage, the specimens that are to be in the low moisture con- tent groups are placed in an environmental chamber at a temperature of 20±3◦C 97 6.2. THE FPINNOVATIONS EXPERIMENT AND DATA ANALYSIS and a relative humidity of 65±5% until no changes in weight are detected. These specimens equilibrate at approximately 8% to 11% moisture content. Similarly, the specimens to be in the high moisture content groups are placed in another en- vironmental chamber at a temperature of 20±3◦C and a relative humidity of 85% until no changes in weight are detected. These specimens equilibrate at approxi- mately 12% to 15% moisture content. In the ramp load test, a specimen is placed between a moving actuator in the centre and immobile supports at both ends (see Figure 1.2). The deflection of the specimen, which is defined as the movement of one point on the specimen relative to another point on the specimen, is controlled to be increasing linearly at a constant rate k, which leads to the failure time Ts around one minute. The deflection rate k is the same for the 20 specimens in each group. But for different combinations of product and condition, the deflection rates k’s are different since the strengths of boards among these combinations are different. However, the resulting breaking times are always around 1 minute. Should we mention that the ramp load test at FPInnovation is done under de- flection control rather than load control. This means that the applied load τ does not increase linearly in time t, as it does in load control ramp tests. Our collabo- rators at FPInnovations use deflection control in the experiments since deflection control is widely used in experiments for testing the strength of wood products. However, descriptions in the literature indicate that load control is widely used in experiments for the duration of load effect, even though load control is not widely used for other experiments. As we discussed in Chapter 2, all accumulation of damage models are based on the load history τ(t), and this is different from the deflection history. An example of the deflection history and the load history of one specimen is shown in Figure 6.1. As we can see from Figure 6.1, when the deflection increases linearly in time t, the load does not increases linearly in time t all the time. As a result, we cannot use the solutions for breaking times in the damage accumulation models as these solutions assumed that τ(t) = kt. In principle, we can still use the 98 6.2. THE FPINNOVATIONS EXPERIMENT AND DATA ANALYSIS damage accumulation models with τ as pictured in Figure 6.1, but the solutions for the breaking times are complicated. 0 10 20 30 40 50 60 70 0 20 40 60 80 Deflection History Time (second) D ef le ct io n (m m) 0 10 20 30 40 50 60 70 0 50 0 10 00 15 00 20 00 Load History Time (second) Lo ad (N ) Figure 6.1: The deflection history and the load history of one specimen in the ramp load test at FPInnovations. In the constant load test, the constant load level is chosen differently as we described in the previous sections. In this experiment, the constant load level is chosen to be 75% of the mean short term strength in the one-minute ramp load tests, and this load is applied to the 20 samples in each of the two Douglas-fir lumber constant load test groups (see Figure 1.3). The constant load tests have not been conducted to the laminated veneer lumber at time of writing this thesis. The constant load tests last for 10 days. 6.2.2 Data Analysis for Experiments at FPInnovations In this section, we analyze the trial data from the experiments at FPInnovations. As noted in the previous section, we cannot use the accumulation of damage mod- els to analyze these data due to the different setups. We compare the empirical cumulative distribution functions of the breaking times Ts’s in the ramp load tests and the breaking times Tc’s in the constant load test from different groups. We fit the breaking times by some distributions. 99 6.2. THE FPINNOVATIONS EXPERIMENT AND DATA ANALYSIS Figure 6.2 contains the empirical cumulative distribution functions of the break- ing times Ts’s (in log-scale) in the ramp load test for each of four combinations of product and condition. Figure 6.3 contains quantile-quantile plots of the breaking times Ts’s (in log-scale) and the fitted log-normal distributions for the four groups. Figure 6.4 shows the empirical cumulative distribution functions of the breaking times Tc’s (in log-scale) in the constant load test for the Douglas-fir lumber groups with high moisture content and with low moisture content respectively. 3.0 3.5 4.0 4.5 5.0 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 log(Ts) Pr op or tio n of th e br ok e n b oa rd s l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l DF−High DF−Low LVL−High LVL−Low Figure 6.2: Empirical cumulative distributions of log(Ts) from the ramp load test results at FPInnovations. In the plot, DF-High, DF-Low, LVL- High and LVL-Low denote Douglas-fir lumber with high moisture content, Douglas-fir lumber with low moisture content, laminated- veneer lumber with hight moisture content and laminated-veneer lum- ber with low moisture content respectively. Figure 6.2 shows that the empirical cumulative distribution functions of log(Ts) for the two Douglas-fir lumber groups with high moisture content and the low moisture content have the same shape but differ by a horizontal shift. The empir- 100 6.3. FOSCHI AND BARRETT’S EXPERIMENTS AND DATA ANALYSIS ical cumulative distribution functions of the two laminated-veneer lumber groups also have similar shapes but differ by a horizontal shift. From Figure 6.2, speci- mens with low moisture content break sooner than those with high moisture con- tent. However, the difference caused by the moisture content for the laminated- veneer lumber is much smaller than that for the Douglas-fir lumber. Figure 6.2 also shows that the empirical cumulative distribution functions of log(Ts) of the Douglas-fir lumber groups and of the laminated-veneer lumber groups differ slightly in shape. The breaking times Ts’s in the two laminated- veneer lumber groups have a smaller variance than those in the Douglas-fir lum- ber groups, which means that the quality of the laminated-veneer lumber is more consistent than the quality of the Douglas-fir lumber, as expected. Figure 6.3 shows that the log-normal distribution fits the breaking times Ts’s reasonably well for each of the four groups. Note that for the Douglas-fir group with low moisture content, two outliers were removed before fitting the data. Figure 6.4 shows that during the 10-day constant load test, more boards of the Douglas-fir lumber group with low moisture content break than those of the Douglas-fir lumber group with high moisture content. Figure 6.4 and Figure 6.2 show that the Douglas-fir lumber is stronger when moisture content is high than when moisture content is low. In conclusion, in the ramp load test, the Douglas-fir lumber and the laminated- veneer lumber are both stronger when moisture content is high. The log-normal distribution fits the breaking times Ts well for each of the four groups. In the 10-day constant load tests, the Douglas-fir lumber with high moisture content is also stronger in the 10-day constant load tests with high moisture content than the Douglas-fir lumber with low moisture content. 6.3 Foschi and Barrett’s Experiments and Data Analysis In this section, we discuss Foschi and Barrett (1982)’s experiments for the du- ration of load effects and analyze the breaking times Tc’s from their constant 101 6.3. FOSCHI AND BARRETT’S EXPERIMENTS AND DATA ANALYSIS l l l ll ll lll ll ll ll l l l l 3.8 4.0 4.2 4.4 4.6 4.8 3. 8 4. 2 4. 6 DF−High log(Ts) Lo ga rit hm o f t he fi tte d di st rib u tio n l l l ll ll ll l ll l l l l 3.7 3.9 4.1 4.3 3. 7 3. 9 4. 1 4. 3 DF−Low log(Ts) Lo ga rit hm o f t he fi tte d di st rib u tio n l l l l ll ll ll lll ll l l l l l 3.9 4.0 4.1 4.2 4.3 3. 90 4. 00 4. 10 4. 20 LVL−High log(Ts) Lo ga rit hm o f t he fi tte d di st rib u tio n l l l lll ll ll lll lll l l l l 3.8 3.9 4.0 4.1 4.2 3. 85 4. 00 4. 15 LVL−Low log(Ts) Lo ga rit hm o f t he fi tte d di st rib u tio n Figure 6.3: The quantile-quantile plots of the breaking times Ts’s (in log- scale) and the fitted log-normal distributions for Ts’s in each group. load tests, where the load level is set to be the 20-th percentile of the short term strength. 6.3.1 Foschi and Barrett’s Experiments In this section, we briefly review Foschi and Barrett’s experiments. In 1980s, Foschi and Barrett conducted an experiment for the duration of load effects on western hemlock lumber of size 2 inches by 6 inches, and grade No. 2 and better. They conducted a ramp load test with a sample size of 150, and then conducted two constant load tests: one with the load level set to be the 20-th percentile and 102 6.3. FOSCHI AND BARRETT’S EXPERIMENTS AND DATA ANALYSIS −4 −2 0 2 4 6 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 log(Tc) Pr op or tio n of th e br ok e n b oa rd s l l l l l l l l l l l l DF−High DF−Low Figure 6.4: Empirical cumulative distribution functions of log(Tc) from the constant load test results in the FPInnovations. In the plot, DF-High and DF-Low denote Douglas-fir lumber with high moisture content and Douglas-fir lumber with low moisture content respectively. There are 20 boards in each group. the other with load level set to be the 5-th percentile of the short-term strength. The initial sample size for each constant load test was 500. However, some boards were discontinued after three months for some reason that was not reported. So the final sample sizes in the one year experiment are 400 boards for the 20-th percentile constant load group and 300 boards for the 5-th percentile constant load group. More details about Foschi and Barrett’s experiments can be found in Foschi and Barrett (1982). In the constant load test of Foschi and Barrett’s experiment when p = 0.2, 207 boards out of 400 boards broke within one year. Foschi and Barrett did not report how many of them broke during the ramp loading part of the constant load tests nor did they specify the value of T0. However, they wrote that 20% of the boards in this test failed during the ramp loading part, as expected. Since the 79- th, 80-th, 81-st and 82-nd order statistics of the breaking times Tc’s are the same 103 6.3. FOSCHI AND BARRETT’S EXPERIMENTS AND DATA ANALYSIS in this dataset, we assume that 82 boards broke during the ramp loading part of the constant load test in our analysis, and assume the 82-nd order statistic of the breaking times Tc’s to be T0. 6.3.2 Data Analysis for Foschi and Barrett’s Experiments In this section, we analyze the breaking times Tc’s from Foschi and Barrett’s ex- periments when the load level is set to be the 20-th percentile of the short term strength. We cannot use the approximate maximum likelihood estimates to estimate parameters for this dataset since the approximate maximum likelihood estimates may not be reliable, as discussed in Section 5.3. On the other hand, we cannot apply the quantile methods and the combined method proposed in Section 4.5 directly to this dataset either, since we do not have the ramp load test results from their experiments. As an alternative, we propose a revised method that combines the approximate maximum likelihood and the quantile method using the breaking times Tc’s only. The only difference between the revised combined method here and the combined method in Section 4.6, is a new way to estimate µX . We recall that the breaking time Tc can be written as Tc = { X , if X ≤ T0, Y, if X > T0, (6.1) where X = exp(a)b exp(b)−1 , (6.2) and Y = T0+ exp(a) exp(b)−1 [ exp { b−T0 exp(b)−1exp(a) } −1 ] . (6.3) In the first step of the combined method discussed in Section 4.6, we estimate 104 6.3. FOSCHI AND BARRETT’S EXPERIMENTS AND DATA ANALYSIS µa and µb from µX and µY , as in (4.19). We estimate µX by the median of log(Ts), and estimate µY by the median of log(Tc−T0). In the second step of the combined method, we estimate σa and σb by the approximate maximum likelihood using µ̂a and µ̂b from the first step. For this dataset, to estimate µY , we can still use the median of log(Tc−T0) as discussed in Section 4.6. To estimate µX , we cannot use the median of log(Ts) since we do not have those Ts’s. However, by definition, those Tc’s, which are less than T0, equal X . Thus, we fit the first 82 order statistics of the logarithm of the breaking times Tc’s to a truncated Normal distribution using the function survreg in R, and estimate µX by the mean of that Normal distribution. After we estimate µa and µb from µ̂X and µ̂Y , we estimate σa and σb by max- imizing the approximate likelihood function for Tc’s as discussed in Section 4.5, using µ̂a and µ̂b from the first step. These steps are all the same as the combined method in Section 4.6. For the optimization process in the second step discussed in Section 4.6, we choose the starting points for σa and σb randomly from the uniform distribution in (0,1). We repeat the optimization process 100 times with different starting points. The results of the optimization show that 89 runs out of 100 runs converge. We estimate σa and σb by the medians of the 89 estimates of σa and σb in these 100 runs. The box plots of the 89 estimates of σa and σb are shown in Figure 6.5. l 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 σa l 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 σb Figure 6.5: Box plots of the 89 estimates of σa and σb. 105 6.3. FOSCHI AND BARRETT’S EXPERIMENTS AND DATA ANALYSIS The resulting estimates are µ̂a = 41.6997, µ̂b = 49.6804, σ̂a = 0.4057 and σ̂b = 0.3105. These estimates of µa and µb are close to Gerhards and Link’s estimates of a and b (â = 43.17 and b̂ = 49.75) when a and b are considered as fixed in their approach. It may be possible to calculate the standard errors of our estimates by the bootstrap method. We do not implement that here. 106 Chapter 7 Conclusions and Future Work 7.1 Conclusions In this thesis, we have reviewed and discussed the accumulation of damage models for the duration of load effects of wood products. We have proposed a method to revise the U.S. model and the Canadian model so that the parameter values do not change when we change units of measurement. We have shown that, for the U.S. model, the distributions of the breaking time Ts’s in the ramp load test and the breaking time Tc’s in the constant load test are not greatly affected by the distributional assumptions for the random effects a and b if the standard deviations of both random effects a and b are relatively small, i.e., if the boards are relatively homogenous. We have shown that the locations of the distributions of the simulated Ts’s and Tc’s are mainly affected by the means of the random effects a and b, and the shapes of the distributions of the simulated Ts’s and Tc’s are mainly affected by the standard deviations of the random effects a and b. Based on simulation studies, we have shown that, breaking times Ts’s gener- ated from the U.S. model seem to approximately follow a log-normal distribution. Breaking times Tc’s generated from the U.S. model or the Canadian model, re- stricted toTc > T0 and Tc − T0 not close to 0, seem to approximately follow a 107 7.2. FUTUREWORK Weibull distribution. For parameter estimation in the U.S. model, we have proposed a maximum likelihood based method, a quantile estimates based method, and a method that combines the maximum likelihood and the quantile estimates for parameter esti- mation in the U.S. model. Through simulation studies, we showed that the max- imum likelihood based method does not work well, the quantile estimates based method estimates the means of the random effects well, but requires an unreason- ably large sample size to estimate the standard deviations of the random effects, and the combined method works very well with a reasonable sample size. We have analyzed the data from the trial experiments conducted at FPInnova- tions. We revised the combined method to estimate parameters in the U.S. model using the breaking times Tc’s from Foschi and Yao’s experiments. Our estimates of the means of the random effects are close to Gerhards and Link’s estimates of the same effects, which they considered as fixed. We also estimated the standard deviations of the random effects. 7.2 Future Work Our estimation methods are based on the same designs of the experiments as Fos- chi and Yao (1986). However, there may be a better design of the experiments which can help to better estimate the parameters. For example, the ramp load tests with different loading rates may be conducted. The damage accumulation models can still be used when the data are collected under the deflection control instead of the load control, however, the solutions of the breaking times are complicated, as discussed in Section 6.1. New models may be developed based on the deflection history of the sample instead of the load history. Our method gives estimates of the parameters in the U.S. model. In the fu- ture, we can extend the method to give the standard errors of the estimates by the bootstrap or other methods. We discussed the possible reasons for the bad performance of the approximate 108 7.2. FUTUREWORK maximum likelihood estimates. There may be potential identifiability issues in the models, i.e., the parameters may not be identifiable in the U.S. model. Our work here mainly focused on the U.S. model. Future work may include the Canadian model, which has more parameters to estimate. The problem of over-parameterization and the problem of non-identifiability may arise in param- eter estimation for the Canadian model. 109 Bibliography [1] Barrett, J. D. and Foschi, R. O. Duration of load and probability of failure in wood. Part I: modelling creep rupture. Canadian Journal of Civil Engeneer- ing, Vol. 5, No. 4: 505-514, 1978. [2] Barrett, J. D. and Foschi, R. O. Duration of load and probability of failure in wood. Part II: constant, ramp and cyclic loadings. Canadian Journal of Civil Engeneering, Vol. 5, No. 4: 515-532, 1978. [3] Cai, Z., Rosowsky, D. V., Hunt, M. O. and Fridley, K. J. Comparison of actual vs. simulated failure distributions of flexural wood specimens subject to 5-day load sequences. Forest Products Journal, 50(1): 7480, 2000. [4] Foschi, R. O. and Barrett, J. D. Load-duration effects in western hemlock lumber. ASCE Journal of the Structural Division, Vol. 108(7): 1494-1510, 1982. [5] Foschi, R. O. and Yao, F. Z. Another look at three duration of load models. In Proceedings of IUFRO Wood Engineering Group Meeting, Florence, Italy, paper: 19-9-1, 1986. [6] Gerhards, C. C. Effects of duration and rate of loading on strength of wood and wood-based materials. Madison: USDA Forest Services Research Report, FPL 283, 1977. [7] Gerhards, C. C. Time-related effects of loading on wood strength: a linear cumulative damage theory. Wood Science, 11(3): 139-144, 1979. 110 BIBLIOGRAPHY [8] Gerhards, C. C. and Link, C. L. A cumulative damage model to predict load duration characteristics of lumber. Wood and Fiber Science, 19(2): 147-164, 1987. [9] Hendrickson, E. M., Ellingwood, B. and Murphy, J. Limit state probabilities for wood structural members. Journal of Structural Engineering, Vol. 113, No. 1, 1987. [10] Markwardt, L.J. and Liska, J.A. Speed of Testing of Wood. Factors in its Control and its Effect on Strength. Proc. ASTM, 48:1139-115, 1948. [11] Wang, J. B. Duration-of-load and creep behavior of thick MPB strand based wood composite. Ph.D. thesis. Faculty of Forestry, the University of British Columbia, 2009. [12] Wood, L. W. Relation of strength of wood to duration of load. Madison: USDA Forest Service Report, No. 1916, 1951. [13] Yao, F. Z. Reliability of structures with load history-dependent strength of an application to wood members. Master thesis. Department of Civil Engi- neering, the University of British Columbia, 1987. 111
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Dynamic duration of load models
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Dynamic duration of load models Zhai, Yongliang 2011
pdf
Page Metadata
Item Metadata
Title | Dynamic duration of load models |
Creator |
Zhai, Yongliang |
Publisher | University of British Columbia |
Date Issued | 2011 |
Description | The duration of load effect is a distinctive and important characteristic of wood strength. It refers to the fact that wood products can usually sustain a high load for a short time but the products may deteriorate and break in the long run. Modelling the duration of load effect and testing wood for specific properties of this effect are important in formulating wood construction standards. Damage accumulation models have been proposed by authors to model the duration of load effects. The models assume that damage is accumulated over time according to the load history, and once the accumulated damage reaches a threshold value, the board will break. Different authors have designed different experiments and proposed different methods for estimating the model parameters. In this work, we consider several damage accumulation models, with a focus on the U.S. model. We investigate the effects of the distributional assumptions for the models, and propose several methods to estimate parameters in the models. Our proposed methods are evaluated via simulation studies. Two real datasets are present for illustration. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-08-29 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0072089 |
URI | http://hdl.handle.net/2429/36958 |
Degree |
Master of Science - MSc |
Program |
Statistics |
Affiliation |
Science, Faculty of Statistics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2011-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2011_fall_zhai_yongliang.pdf [ 5.92MB ]
- Metadata
- JSON: 24-1.0072089.json
- JSON-LD: 24-1.0072089-ld.json
- RDF/XML (Pretty): 24-1.0072089-rdf.xml
- RDF/JSON: 24-1.0072089-rdf.json
- Turtle: 24-1.0072089-turtle.txt
- N-Triples: 24-1.0072089-rdf-ntriples.txt
- Original Record: 24-1.0072089-source.json
- Full Text
- 24-1.0072089-fulltext.txt
- Citation
- 24-1.0072089.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0072089/manifest