Designing experiments for estimating an appropriate outlet size for a silo type problem
DDesigning experiments for estimating an appropriateoutlet size for a silo type problem
Jesus Lopez-Fidalgo, Caterina May, Jose Antonio Moler
Abstract:
The problem of jam formation during the discharge by gravity of granu-lar material through a two-dimensional silo has a number of practical applications.In many problems the estimation of the minimum outlet size which guarantees thatthe time to the next jamming event is long enough is crucial. Assuming that thetime is modeled by an exponential distribution with two unknown parameters, thisgoal translates to the optimal estimation of a non-linear transformation of the pa-rameters. We obtain c -optimum experimental designs with that purpose, applyingthe graphic Elfving method. Since the optimal designs depend on the nominalvalues of the parameters, a sensitivity study is additionally provided. Finally, asimulation study checks the performance of the approximations made, first withthe Fisher Information matrix, then with the linearization of the function to beestimated. The results are useful for experimenting in a laboratory and translat-ing then the results to a larger scenario. Apart from the application a generalmethodology is developed in the paper for the problem of precise estimation of aone-dimensional parametric transformation in a non-linear model. Keywords:
Elfving graphical procedure; Exponential probability model; FisherInformation Matrix; Granular material; Linearization; Non-linear parameter trans-formation. 1 a r X i v : . [ s t a t . M E ] S e p Introduction
Material in granular form appears in many contexts of applications, as in the phar-maceutical, chemical, food, agricultural and mining industry (see Nedderman,1992). During the discharge by gravity of this material through an outlet, if thesize of the outlet is not large enough, the formation of an arch at some point usu-ally interrupts the flow, causing a jam . An arch is defined as a structure consistingof particles which are mutually stabilized (Janda et al., 2008) until an externalinput of energy breaks their blocking structure and restarts the flow until the nextjam happens.The problem of jam formation during the discharge by gravity of granularmaterial through a two-dimensional silo has been studied in Janda et al. (2008),Amo-Salas et al. (2016b) and Amo-Salas et al. (2016a). In particular, they fo-cus on studying the waiting time that passes between two jamming events, whichdepends on the outlet size, according to some model. This waiting time is alsorelated with the avalanche , that is the amount of material dropped between twojamming events. In Amo-Salas et al. (2016b) and Amo-Salas et al. (2016a) theoptimal experimental designs to estimate the unknown parameters and to discrim-inate between models are obtained.There is a common interest in avoiding a jam at least during a specific periodof time. In fact, the event of breaking the arches may be dangerous, expensiveor just no affordable. Hence, the goal of this paper is the precise estimation ofthe minimum outlet size necessary to guarantee that the expected time betweentwo jamming events will exceed a fixed time of interest. Assuming an exponen-tial model as in Amo-Salas et al. (2016b), this goal determines the problem offinding an optimal design to estimate a non-linear transformation of the unknown2arameters.When the inferential goal is the estimation of a linear combination of the un-known parameters, a c -optimal design minimizes the variance of the maximumlikelihood estimator (classical references on optimal designs are, for instance,Atkinson et al. (2007) and Pukelsheim (2006)). Elfving (1952) provided a graph-ical method to determine c -optimal designs of a linear model on a compact ex-perimental domain, based on the construction of a convex hull . This method isnot easy to use for more than two parameters, but L ´opez-Fidalgo and Rodr´ıguez-D´ıaz (2004) provided an iterative procedure for more than two parameters basedon the graphical Elfving technique. For instance L´opez-Fidalgo and Rodr´ıguez-D´ıaz (2004) used successfully this procedure to compute c–optimal designs formore than two parameters. Since the model considered is non-linear it is possi-ble to determine a local c -optimal design by considering the Fisher InformationMatrix (FIM) for nonlinear models and a first–order linearization of the functionof the parameters to be estimated around some nominal values of the parameters.Moreover, we adopt c -optimality for estimating a non-linear transformation of theparameters by linearizing also this function.The paper is organized as follows. Section 2 presents the problem, sets thebasic notation and explains our method in its generality. Section 3 contains theresults on c -optimal designs. Section 4 provides a sensitivity analysis. Section 5contains a simulation study to check the validity of the approximations applied toobtain the results. Section 6 concludes the paper. All the computations have beendone with Python 3.7. Codes are provided as supplementary material.3 Problem and general method
Consider the problem of falling of particles through a two-dimensional silo pre-sented in Amo-Salas et al. (2016a) and Amo-Salas et al. (2016b) and introducedin Section 1. Denoting by T the time between two jamming events and by φ thesize of the outlet at the bottom of the silo, let E [ T | φ ] = η ( φ ; θθθ ) (1)be the mean time between two jamming events, where θθθ represents the unknownmodel parameter. Following Amo-Salas et al. (2016a) and Amo-Salas et al. (2016b),it is realistic to consider that T has an exponential distribution. In particular, giventhe outlet size φ of the silo, which is a controlled variable, in the next Section wewill assume for the mean function (1) the model in (Amo-Salas et al., 2016a, eq.(3)).There is a common interest in avoiding a jam at least in a period of time.Hence, our goal is the precise estimation of the minimum outlet size necessary toguarantee that the expected time between two jamming events will be greater thana fixed time T of interest: E [ T | φ ] ≥ T . (2)If η ( · , · ) is an invertible function, (2) becomes φ ≥ g ( T , θθθ ) , (3)for some inverse function g . Thus, we are interested in estimating g ( T , θθθ ) whichis a non-linear function of the unknown model parameter. Since T is a fixedconstant from now on we denote it simply by g ( θθθ ) .4o this aim, assume that an experimenter can observe uncorrelated observa-tions from n experiments, t i = η ( φ i ; θθθ ) + ε i , i = , ..., n . (4)Since φ is a controlled variable, the n experimental conditions φ , ..., φ n can bechosen according to a design ξ , that is, a probability distribution on a domain X = [ a , b ] : ξ = (cid:26) φ · · · φ r p · · · p r (cid:27) , with r ≤ n .Observe that the model herein considered is non-linear and the errors ε i havenon-constant variance: Var ( ε i ) = η ( φ i ; θθθ ) The FIM is defined by M ( ξ , θθθ ) = (cid:90) X I ( φ , θθθ ) d ξ ( φ ) , where I ( φ , θθθ ) = − E Y (cid:20) ∂ ∂ θθθ L ( θθθ ; t , φ ) (cid:21) is a two by two matrix and L is the log-likelihood function. Since, for an expo-nential model with mean (1), we have L ( θθθ ; t , φ ) = log (cid:18) η ( φ , θθθ ) exp − t η ( φ , θθθ ) (cid:19) , (5)it follows that the FIM of model (4) at one point φ is I ( φ , θθθ ) = η ( φ , θθθ ) ∇ η ( φ , θθθ ) ∇ η ( φ , θθθ ) T , M ( ξ , θθθ ) = (cid:90) X η ( φ , θθθ ) ∇ η ( φ , θθθ ) ∇ η ( φ , θθθ ) T d ξ ( φ ) , (6)where the transpose is indicated with the superscript T and ∇ stands for the gra-dient.Equation (6) is also the FIM of the following linear gaussian and homoschedas-tic model t i = θθθ T f ( φ i ; θθθ T ) + ε i , (7)with f ( φ ; θθθ ) = η ( φ , θθθ ) ∇ η ( φ , θθθ ) , (8)Our goal is therefore to find an optimal design for precise estimation of g ( θθθ ) ,that is, a design minimizing the variance of the maximum likelihood estimator(MLE) of g ( θθθ ) .When the inferential goal of an experiment is an efficient estimation of a vec-tor of unknown parameters θθθ , an optimal design maximizes a suitable functionalof the FIM, M ( ξ , θθθ ) , because its inverse is asymptotically proportional to the co-variance matrix of the MLE ˆ θθθ , which is asymptotically unbiased. Some classicalreferences on optimal designs are Fedorov (1972), Pzman (1986) and Atkinsonet al. (2007). An optimal design depends on the value of the unknown parametersexcept in the case of linear models.As mentioned above a c -optimal design ξ ∗ c , minimizes the asymptotic varianceof a linear transformation c T θθθ of the unknown model parameters: ξ ∗ c = arg min ξ c T M ( ξ ; θθθ ) − c (9)A very nice way to compute c -optimal designs, especially in two dimensions,is the geometric Elfving procedure (see Elfving, 1952). Such procedure is con-6tructed for estimating a linear transformation c T θθθ of the parameters given a linearhomoschedastic model T = θθθ T f ( φ ) + ε .Remembering that the MLE estimator of g ( θθθ ) is g ( ˆ θθθ ) , let us approximate thenon-linear function g ( ··· ) using Taylor expansion around the true value θθθ t , so thatwe can approximate g ( ˆ θθθ ) with g ( θθθ t ) + ∇ g ( θθθ t )( ˆ θθθ − θθθ t ) . The variance of g ( ˆ θθθ ) canbe then approximated by ∇ g ( θθθ t ) T M ( ξ , θθθ t ) − ∇ g ( θθθ t ) (10)and a c -optimal design for model (4) is a design satisfying (9) with c = c ( θθθ ) givenby c ( θθθ ) = ∇ g ( θθθ ) (11)Notice that two procedures of approximation by linearization have been adopted,and that the c -optimum design satisfying (9) depends on the unknown parametersboth through the vector c and the FIM, M ( ξ , θθθ ) . Hence, a nominal value θθθ guessing the true value θθθ t has to be chosen and the design obtained will be locallyoptimum. Starting from the design space considered in Janda et al. (2008) andthe values obtained in Amo-Salas et al. (2016a), the procedure to obtain c -optimaldesigns is developed in detail in the next section. c -optimal designs Assume that the time T between two jamming events is exponentially distributedwith mean η ( φ ; θθθ ) = C exp ( L φ ) − , φ ∈ X = [ a , b ] , (12)where θθθ T = ( C , L ) , as in (Amo-Salas et al., 2016a, eq. (3)).7ur main goal is the efficient estimation of the minimal diameter φ ∈ X forwhich (2) holds. If the mean of T is given by (12), this means η ( φ ; θθθ ) = C exp ( L φ ) − ≥ T , (13)that is φ ≥ g ( θθθ ) = (cid:114) log ( C ( T + )) L . (14)The extremes of the experimental domain X = [ a , b ] have to satisfy d < a < b < φ C , where d is the diameter of the granular material and φ C is a nominal di-ameter above which jamming is practically impossible. In theory there is not sucha value since there is always a chance of forming an arch, no matter how widethe outlet is. But a practical limit can be assumed and even that value is a param-eter of some models (see Amo-Salas et al., 2016a and Amo-Salas et al., 2016b).Moreover, L > η , between jams is increasing withrespect to φ ; and 0 < C < exp ( L φ ) for any φ since η must be positive. Thisdoes not mean constrained estimation, but just practical limits. The data will bein charge of dealing with them. Following the method presented in Section 2, weare obtaining here c -optimal designs for estimating the bound g ( θθθ ) given in (14). Remark 1
An alternative goal could be the estimation of the minimal diametersuch that, for a given value α , P ( T > T ) ≥ − α . However, for an exponentialmodel, this is equivalent to consider (14) with T = − T / log ( − α ) since − α ≤ P ( T > T ) = exp ( − T / η ( φ ; θθθ )) (15) and then η ( φ ; θθθ ) ≥ − T log ( − α ) . (16)8 or instance, for α = . , T = . T . Thus, the problem is equivalent in termsof estimation and designing and this is the relationship between both thresholds.In this particular case the threshold for the probability is about 20 times the onefor the expectation. The information matrix at a point φ for model (12) is I ( φ , θθθ ) = e φ L C ( e φ L − C ) (cid:32) C − φ − φ C φ (cid:33) . (17)In order to apply the Elfving’s graphical method, we need to obtain the Elfvinglocus, that is, the convex hull of the union of the curve defined by the regressors in(7) and its reflection through the origin; the c -optimal design is then determinedby the crossing point between the line c and the boundary of the Elfving locus(see L ´opez-Fidalgo and Rodr´ıguez-D´ıaz, 2004). The parametric equations thatrepresent the curve are obtained from (8), which becomes, when model (12) isassumed, f ( φ , θθθ ) = G ( φ , θθθ ) ( / C , − φ ) T , where G ( φ , θθθ ) = e φ L e φ L − C . (18)Hence, the parametric equations of f ([ a , b ]) are x ( φ ) = G ( φ , θθθ ) / C , y ( φ ) = − G ( φ , θθθ ) φ , φ ∈ [ a , b ] . (19)According to the experimental case considered in Janda et al. (2008), φ ∈ X = [ . , . ] and the estimates obtained in Amo-Salas et al. (2016a) fromdata will be used as nominal values of the parameters, that is, C = . L = . A A A A obtained in this case. It is worth observing thatthe vertexes of the convex hull in Figure 1 are not tangential points of the curvebut outermost points of the curve. The vector c , defined as the gradient of g ( θθθ ) Figure 1: Convex hull based on the estimates from Amo-Salas et al. (2016a).Green sides represent the possible crossing points of ∇ g ( θθθ ) .evaluated in the nominal values ( C , L ) , is given by c ( θθθ ) T = √ L (cid:32) C (cid:112) log ( C ( T + )) , − (cid:112) log ( C ( T + )) L (cid:33) ; (20)depending on the value of T , c has a different angle and the line directed by ccc crosses the convex hull in A A or A A (equivalently A A or A A ) (see Figure2). The following proposition provides the properties of the Elfving locus for anyvalues of the extremes of the experimental domain and for any possible choice ofthe nominal values of the parameters, ( C , L ) . Proposition 1
Consider the curve (19) and its reflection through the origin. LetA , A , A , A be the outermost points: A = ( − x ( b ) , − y ( b )) , A = ( x ( a ) , y ( a )) , = ( x ( b ) , y ( b )) , A = ( − x ( a ) , − y ( a )) .Then, for any value of ( C , L ) such that < C < exp ( L φ ) and L > , the convexhull of these curves is A A A A . Proof 1
The main point is to prove that the curve (19) is always above the segmentA A and below the A A . The reasoning will be organized in the following steps:1. Observe that x ( φ ) > and y ( φ ) < for any φ ∈ [ a , b ] since L > ande φ L > C > for any φ . Then, the curve (19) and the points A , A are al-ways in the fourth quadrant of the Cartesian plane, while its reflection andthe points A , A are always in the second quadrant.2. We have x (cid:48) ( φ ) < for any φ ∈ [ a , b ] ; it follows that x ( b ) ≤ x ( φ ) ≤ x ( a ) .3. From the first equation of (19) we have G ( φ , θθθ ) = Cx ( φ ) ; moreover, by thedefinition of G ( φ , θθθ ) , φ = L log (cid:18) C xCx − (cid:19) ; then plugging into the second equation of (19), we obtain the cartesianequation of the curve:y ( x ) = − CL x log (cid:18) C xCx − (cid:19) , x ∈ [ x ( b ) , x ( a )] . (21)
4. Notice that y ∈ C ([ x ( b ) , x ( a )]) and thaty (cid:48)(cid:48) ( x ) = − CLx ( Cx − ) < it follows that (21) is concave and therefore (19) is above the segment A A . . In order to prove that (19) is below the segment A A it is enough to provethat the tangent to the curve in A is below A A (which has a negative slopem); this means that the slope of (21) in x = x ( a ) is greater than the slope ofA A .We have y (cid:48) ( x ) = CL (cid:18) Cx − − log C xCx − (cid:19) (22) and then y (cid:48) ( x ) | x = x ( a ) = CL (cid:18) G ( a , θθθ ) − − log C G ( a , θθθ ) G ( a , θθθ ) − (cid:19) = CL (cid:32) e a L − CC − a L (cid:33) = L e a L − CL − a C . (23) At this point there are two cases:(a) If C < e a L / ( + a L ) then y (cid:48) ( x ) | x = x ( a ) > and it is straightforwardthat the slope of the curve is greater than the slope of A A . Note thatin this case we have y (cid:48) ( x ) > for any x ∈ [ x ( b ) , x ( a )] (as in Figure 1).(b) If e a L / ( + a L ) < C < e φ L , then y (cid:48) ( x ) | x = x ( a ) < (as in Figure 3) andwe have to prove thaty (cid:48) ( x ) | x = x ( a ) > m = − C a G ( a , θθθ ) + b G ( b , θθθ ) G ( a , θθθ ) + G ( b , θθθ ) (24) Since (23) can be written asCL (cid:18) G ( a , θθθ ) − − a L (cid:19) , he inequality (24) is equivalent to L G ( a , θθθ ) − − a > − a G ( a , θθθ ) + b G ( b , θθθ ) G ( a , θθθ ) + G ( b , θθθ ) , which gives L G ( a , θθθ ) − + ( b − a ) G ( b , θθθ ) G ( a , θθθ ) + G ( b , θθθ ) > , which is always satisfied since the left term is a sum of two positive quanti-ties. Denote by ( x i , y i ) , i = , ...,
4, the coordinates of the extremes A i of the con-vex hull stated in Proposition 1 and denote by φ i the corresponding values of φ in curve (19) or in its symmetric (from Proposition 1, φ i can be equal to a orequal to b ). Next proposition gives the c -optimal designs obtained by the crossingpoint between ccc = ∇ g ( θθθ ) and the convex hull, according to the Elfving method,depending on the fixed value T . Proposition 2
Depending on the fixed value of T , the convex hull is crossed by c ( θθθ ) through A i A i + , where i = , , , and the c-optimal design is (cid:26) φ i φ i + − p i p i (cid:27) , with p i = (cid:115) ( Kx − y i ) + ( x − x i ) ( x i + − x i ) + ( y i + − y i ) , (25) where K = ( ∂ g / ∂ L ) / ( ∂ g / ∂ C ) and the coordinates of the crossing point P arex = y i − y i + − y i x i + − x i x i K − y i + − y i x i + − x i , y = Kx . (26) In particular, let T i = C exp (cid:18) − y i Lx i C (cid:19) − , then for T ∈ ( max ( , ( − C ) / C ) , T ] the crossing point is in A A ; • for T ∈ ( T , T ] the crossing point is in A A ; • for T > T the crossing point is in A A . Proof
Observe that the lines that contain the segment A i A i + and c ( θθθ ) can berespectively written as y − y i = y i + − y i x i + − x i ( x − x i ) and y = Kx , hence the solution of the crossing point (26) follows straightforwardly.From the Elfving method we have that if the crossing point P is in the side A i A i + , then the c -optimal design is given by (25) with p i = (cid:107) A i P (cid:107) / (cid:107) A i A i + (cid:107) ,where (cid:107) · (cid:107) is the euclidean norm.Finally, taking into account that c ( θθθ ) is given by (20); as log ( C ( T + )) > T > − CC . Moreover, since ∂ g / ∂ L < ∂ g / ∂ C > c ( θθθ ) always moves into the fourthquadrant. As only the vertices A and A can be in the fourth quadrant, then P = A i , i = ,
3, are the only two situations where the optimal design reduces toone point. In such a case, y i / x i = K , and then T i = e − y i Lx i CC − . xample 1 Let X = [ . , . ] , C = . , L = . and T = ,then, according to Proposition 2, the convex hull is crossed by c ( θθθ ) in A A andthe optimal design is ξ ∗ c = (cid:26) .
53 5 . . . (cid:27) It is almost equally weighted as the D-optimal design. Figure 2 represents theconvex hull and ∇ g ( θθθ ) in this example. Figure 2: Location of the main points addressed in Proposition 2 for the nominalvalues from Amo-Salas et al. (2016a) and T = Example 2
In the proof of Proposition 1 two situations are distinguished depend-ing on whether point A is a maximum or not. One has been illustrated in Exam-ple 1 and the second one is illustrated in this example. For L = . and φ ∈ [ . , . ] , a value of C in the interval ( . , . ) must be chosen, sayC = . . Here T = .Figure 3 represents the convex hull and ∇ g ( θθθ ) for T = in this example T = C = . L and experimental domain as in Amo-Salas et al. (2016a). where A is not a maximum. Now, Proposition 2 holds, and then ξ ∗ c = (cid:26) .
53 5 . . . (cid:27) It is interesting to stress that this design put more weight in the right extreme,and therefore longer experimentation times are required, although the limit T ismuch smaller than in the previous example. Assume that T is a given value; the following steps describe the procedure toperform a sensitivity study for the choice of the nominal values of the parameters. Step 1 : Consider φ ∈ [ . , . ] and the nominal values ( C , L ) . The c-optimal design is obtained in Proposition 2, ξ ( ) c = (cid:40) φ ( ) i φ ( ) i + − p ( ) i p ( ) i (cid:41) tep 2 : We consider a grid where the parameters C and L take potential actualvalues ( C ∗ , L ∗ ) in a neighborhood of the nominal values ( C , L ) . Thus we obtainthe c-optimal design ξ ∗ c = (cid:26) φ ∗ i φ ∗ i + − p ∗ i p ∗ i (cid:27) , when the true values of the parameters is a pair ( C ∗ , L ∗ ) in the grid. Step 3 : For each ( C ∗ , L ∗ ) in the grid, the following values are obtained: • M = ( − p ∗ i ) I ( φ ∗ i , C ∗ , L ∗ ) + p ∗ i I ( φ ∗ i + , C ∗ , L ∗ ) Var ( g ) = ∇ ( g ( C ∗ , L ∗ )) T M − ∇ ( g ( C ∗ , L ∗ )) • Consider the nominal values ( C , L ) where p ( ) i and φ ( ) i where ob-tained in Step 1 and obtain M = ( − p ( ) i ) I ( φ i , C ∗ , L ∗ ) + p ( ) i I ( φ i + , , C ∗ , L ∗ ) Var ( g ) = ∇ ( g ( C ∗ , L ∗ ) T M − ∇ g ( C ∗ , L ∗ ) • Compute the relative efficiency given by
Var ( g ) / Var ( g ) . Example 3
From Proposition 2, three different situations can be distinguisheddepending on the nominal values chosen for ( C , L ) . In particular, when the nomi-nal values ( C , L ) are as in the Example 1, T can be in the intervals (0.49, 2.57],(2.57, 203603.03] or (203603.03, ∞ ). From here, we consider the following threecases: T = , and , . The first one is too small to have a practicalinterest, while the last one needs a diameter longer than those ones in the designspace. Thus, they are extreme cases, but interesting to be considered in this study.Consider a grid of points ( C ∗ , L ∗ ) appropriate to detect sensitive changes inthe efficiencies. In Figures 4, 5 and 6 the efficiencies for the three cases considered re shown. In all the three cases C ∗ varies in the interval ( C − . , C + . ) whileL ∗ varies in the interval ( L − . , L + . ) in cases 1 and 2 and in the interval ( L − . , L + . ) in case 3.Observe that as we change C ∗ and L ∗ in the grid, also the three intervalsstated in Proposition 2 change. Since the value T is fixed, the crossing point canbe in a different segment A i A i + for ( C , L ) and the point of the grid ( C ∗ , L ∗ ) .For instance, in Figures 4 and 5, the largest decrement of the efficiency happensfor large values of C ∗ combined with small values of L ∗ , and it can be checkedthat they provide values of T smaller than 2 in the case 1 and T values smalleror slightly larger than 200 in the case 2. Finally, in Figure 6, a smaller intervalis chosen to vary L ∗ because dramatic changes of the efficiency are observed forfurther values of L ∗ ; indeed, T is also very sensitive to small changes in theparameters. Now, the efficiency decreases when L ∗ grows and C ∗ decreases (topleft on Figure 6) and when L ∗ decreases and C ∗ grows (bottom right on the table).These two situations correspond, respectively, with values of T much smaller ormuch larger than 300,000. In other words, the cross points of the gradient withthe convex hull are far away from the cross point of ( C , L ) or in other segment.We could say that, in the three cases, when both, L ∗ and C ∗ , grow or decrease,the efficiency is more stable; but changes of C ∗ and L ∗ in opposite directions makethe efficiency to reduce quicker. T = T = T = , Observe that expressions (7) and (10) in Section 2 show the two linear approxi-mation procedures that have been adopted to solve the problem. The main goalof this section is to compare the a priori approximated variances and covariancesof the estimates with the empirical variances and covariances of the estimate ob-tained by simulation. This is to have an insight of the accuracy of linearizing whenlooking for c-optimal designs in non-linear models.The simulations will be performed in the following steps:
Step 1 : Consider the nominal values ( C , L ) and obtain the optimal design(25) for a fixed value T . Following the notation in Proposition 2, n i observationsare randomly allocated at φ i and n i + = n − n i at φ i + . Step 2 : We obtain the MLEs of C , L and g ( θθθ ) . As the responses follow anexponential distribution with mean (12), we have n i responses from an exponentialdistribution with parameter λ i that are denoted by t ( i ) k , k = , · · · , n i and n i + withparameter λ i + , which are denoted by t ( i + ) k , k = , · · · , n i + where20 j = Ce L φ j − C , j = i , i + . (27)The likelihood function depends on the sample obtained, t , and the parametervalues C and L , L n = L n ( t , θ ) = λ n i i e − λ i n i ∑ k = t ( i ) k λ n i + i + e − λ i + n i + ∑ k = t ( i + ) k . By solving the equations ∂ L n / ∂ C = ∂ L n / ∂ L = λ j = n jn j ∑ k = t jk = T j ; j = i , i + C and L :ˆ C = (cid:32) ( + T i + ) φ i ( + T i ) φ i + (cid:33) φ i + − φ i , ˆ L = log (cid:18) ( + T i + )( + T i ) (cid:19) φ i + − φ i (28)The MLE of g ( θθθ ) in (14) is given by g ( ˆ θθθ ) , where ˆ θθθ T = ( ˆ C , ˆ L ) . Step 3 : Step 2 is repeated m times obtaining three m -vectors, ˆC , ˆL , ˆg , whichcontain, respectively, the MLEs of C , L and g ( θθθ ) computed at each step. Step 4 : To study the accuracy of the approximation, the covariance matrixof ˆ θθθ is approximated by the empirical covariance matrix of ( ˆ CCC , ˆ LLL ) . Since theMLE is asymptotically efficient, the covariance matrix of ˆ θθθ should be similarto the Frechet-Cramer-Rao bound for n sufficiently large. In the multiparametercase, this bound is equal to ℑ = ∂ Ψ / ∂ θθθ T × I ( φ , θθθ ) ∗ ∂ Ψ T / ∂ θ , where I ( φ , θθθ ) isdefined in (17) and Ψ ( θθθ ) = E ( ˆ θθθ ) .Observe that ( ∂ Ψ / ∂ θθθ ) i j = ∂ Ψ i / ∂ ( θ j ) = Cov ( ˆ θ j , ∂ log ( L n ) / ∂ θ j )) , where L n is the likelihood function. In order to approximate this matrix, in step 2 we21ill also obtain, in each run, the bidimensional vector: ∂ log ( L n ) / ( ∂ θθθ ) = (cid:18) ∂ log ( L n ) / ∂ C ∂ log ( L n ) / ∂ L (cid:19) = i + ∑ j = i n j ∑ k = C + e L φ j − C + e L φ j ( e L φ j − C ) t ( j ) k − φ j e L φ j e L φ j − C (cid:34) − C e L φ j − C t ( j ) k (cid:35) (29)then, we approximate ( ∂ Ψ / ( ∂ θθθ )) i j with the corresponding sample covariance. Example 4
Consider the setup of Example 1.
Step 1:
Let C = . and L = . as in Janda et al. (2008), andconsider several values of T (see Table 1). Step 2:
We allocate randomly n = experimental points following theoptimal design obtained from Proposition 2. The MLE values of C, L and g ( θθθ ) are obtained jointly with the pair of values of the vector (29) that we denote,respectively, f ( ) n and f ( ) n . Step 3:
Step 2 is repeated m = times and the 1000-dimensional vectors ˆ C, ˆ L, g ( ˆ θθθ ) , f n and f n are stored. Step 4:
Table 1 shows a high similitude between the target value g ( θθθ ) andits MLE ˆ g. Also, between the variance obtained with the simulated Cov ( ˆ CCC , ˆ LLL ) denoted in the table as ˆ Var ( ˆ g ) and the variance obtained with ℑ , which is denotedin the table as Var ( ˆ g ) . The numbers must be multiplied by − .Observe that for T = . nor the estimator, neither the variance are similar.As T is in the interval (0.4887, ∞ ), values close to the boundary carry out a slowerconvergence of the estimators. In Table 2 we study the approach for T = . ofg = . and ˆ g and Var [ ˆ g ] and ˆ Var [ ˆ g ] for increasing values of the sample sizen. < T T < T < T T > T T × × × × g g Var ( g ) ∗
87 5.8 1.4 0.9 0.8 0.6 0.6 0.6 1.0 1.4
Var ( ˆ g ) ∗
562 6.1 1.4 0.9 0.8 0.6 0.6 0.6 1.0 1.4 ∗ The variances must be multiplied by 10 − Table 1: Simulation performance for Example 1 for several values of T n bias = ˆ g − g Var [ ˆ g ] ∗ Var [ ˆ g ] ∗ ∗ The variances must be multiplied by 10 − Table 2: Accuracy of the approximations for different values of n , T = . C = . L = . he decreasing rate is smaller for the bias than it is for the variance. In this paper we consider the problem of estimating the parameters of a non-linearmodel for the time between two jams in the emptying of a silo. This may be ap-plied to a number of phenomena such as delivering some material on a mine ona vertical tunnel. In most of the cases a jam might be rather dramatic involvingsome expense procedure to break the jam. In the case of the mine some explosivehas to be use including risks and delays. Then a very important aim is to deter-mine the diameter of the outlet, say φ , in order to guarantee a period of time longenough. This could be considered as a specific expected time, say T , or else aspecific probability of reaching a specific time without jams. This entails the es-timation of a lower bound expressed as a non-linear function that depends on theunknown parameters and T . For both situations, expected time and probability,give the same function of the parameters to be estimated tuning adequately thethree specific constants mentioned above. In order to obtain an analytical solutionof the problem, first we use the Fisher Information approximation for the covari-ance matrix of the estimates of the paramateres. Then the non-linear lower bound,which is the target for estimation, is linearized being the its gradient the c-vectorfor c–optimality. A model with two parameters is chosen, and, so, the graphicElfving procedure to find the c-optimal design is used.Propositions 1 and 2 establish, respectively, the main characteristics of theconvex hull depending on the parameter values and then an explicit expression forthe c-optimal design. Moreover, the latter indicates that the c-vector may intersectthe convex hull in three sides of the convex hull depending on three intervals where24 can lie. The vertices produce c-optimal designs with only one–point designs,otherwise two points are needed.The vertices of the convex hull are critical points in the sensitivity analysissince they indicate a change of the type of design. For this study a uniform gridwith values for the parameters around the nominal values was considered in orderto detect big changes in the efficiency. A dramatic loss of efficiency happenswhen the parameter values considered in the grid produce a change of edge forthe the crossing point of the c-vector. A smaller decreasing is observed when thecrossing point moves away on the same edge of the convex hull. Both facts implya very important change of the weights of the c-optimal design in Proposition 2.Besides this, for very large values of T , the sensitivity of the design with respectto the selection of the nominal values is large, in fact, a small change of one of theparameters gives place to a dramatic decreasing of the efficiency, this is why thesensitivity study requires a reduced scale on this parameter.A simulation study is carried out to check the accuracy of the double proce-dure to linearize the problem. So that, given the original non-linear model theML estimators are obtained in a simulation procedure with a large number n ofobservations allocated in the c-optimal design, given a T value and usual nominalvalues taken from the literature. Results show very close results, in general, theapproximation procedure produces slightly higher variances of the lower boundfor the silo outlet size than the simulated one. When T is close to its lower bound,the convergence is slower and n must be enlarged.25 uppementary material All the computations have been done with Python 3.7. Codes are provided in twofiles.
Acknowledgements
The first author was sponsored by Ministerio de Ciencia y Tecnolog´ıa MTM2016-80539-C2-1-R and the third one by Ministerio de Ciencia y Tecnolog´ıa MTM2017-83812-P and MTM2016-77015-R. The second author thanks the Departamento deEstad´ıstica, Inform´atica y Matem´aticas of the Universidad P´ublica de Navarra forenduring his scientific visit to the department.
References
M. Amo-Salas, E. Delgado-M´arquez, L. Filov´a, and J L´opez-Fidalgo. Optimaldesigns for model discrimination and fitting for the flow of particles.
Statist.Papers , 57(4):875–891, 2016a.M. Amo-Salas, E. Delgado-M´arquez, and J. L´opez-Fidalgo. Optimal experimen-tal designs in the flow rate of particles.
Technometrics , 58(2):269–276, 2016b.ISSN 0040-1706.A. C. Atkinson, A. N. Donev, and R. D. Tobias.
Optimum experimental designs,with SAS , volume 34 of
Oxford Statistical Science Series . Oxford UniversityPress, Oxford, 2007. ISBN 978-0-19-929660-6.G. Elfving. Optimum allocation in linear regression theory.
The Annals of Math-ematical Statistics , 84(4):44002–1–44002–6, 1952.26. Fedorov.
Theory of Optimal Experiments . Acadimic Press, New York, 1972.A. Janda, I. Zuriguel, A. GArcimartin, L. A. Pugnaloni, and D. Maza. Jammingand critical outlet size in the discharge of a two-dimensional silo.
Eurphysicsletters , 23:255–262, 2008.J. L´opez-Fidalgo and J. M. Rodr´ıguez-D´ıaz. Elfving’s method for m -dimensionalmodels. Metrika , 59(3):235–244, 2004. ISSN 0026-1335.R.M. Nedderman.
Statics and Kinematics of Granular Materials . CambridgeUniversity Press, Cambridge, 1992.F. Pukelsheim.
Optimal design of experiments , volume 50 of
Classics in Ap-plied Mathematics . Society for Industrial and Applied Mathematics (SIAM),Philadelphia, PA, 2006. ISBN 0-89871-604-7. doi: 10.1137/1.9780898719109.Reprint of the 1993 original.A. Pzman.