Likelihood ratios and Bayesian inference for Poisson channels

The Bayesian Central Limit Theorem,with some intructions on#1(e)of the second problem setWe have seen Bayes formula:posterior probability density=constant prior probability density likelihood function.··The“Bayesian Central Limit Theorem”says that under certain circumstances,the posterior probability distribution is approximately a normal distribution.Actually more than one proposition is called“the Bayesian Central Limit Theorem,”and we will look at two of them.First version:SupposeΘis distributed according to the probability density function fΘ(θ), andX1,X2,X3,......[Θ=θ]∼i.i.d.f X1[Θ=θ](x)||i.e.,X1,X2,X3,...are conditionally independent given thatΘ=θ,and f X1[Θ=θ](x),as a|function of x withθfixed,is the conditional probability density function of X1given thatΘ=θ.Suppose the posterior probability density function fΘ[X1=x1,...,X n=x n]is everywhere|twice-differentiable,for all n.Let c=the posterior expected value=E(ΘX1=x1,...,X n=x n),and|let d=the posterior standard deviation=Var(Θ|X1=x1,...,X n=x n).Then if the distribution of U is the posterior distribution,then the distribution of U−cis dapproximately standard normal;more precisely:that distribution approaches the standard normal distribution as n approaches infinity. More tersely:The posterior distribution approaches a normal distribution as the sample-size n grows.How big a value of n is needed to get a good approximation?For a partial answer to that, consider#1(e)on the second problem set.A90%posterior probability interval found without using the Bayesian Central Limit Theorem is(0.419705351626,0.788091845180).These are correct to12digits after the decimal point if I can trust my software.The endpoints of this interval differ from the approximations given by the above version of the Bayesian Central Limit Theorem by slightly less than0.01.In the case of beta distributions,finding the expected value and the standard deviation is trivial:Just use what is says on page306of Schervish&DeGroot.For some other1distributions,finding the expected value and the standard deviation may require laborious integrals,and hence we might use the next version of the Bayesian Central Limit Theorem: Second version:Start with the same assumptions as in thefirst version.Let c=the posterior mode=the value ofθthat maximizes the posterior density function fΘ[X1=x1,...,X n=x n](θ).|Let k(θ)=log fΘ[X1=x1,...,X n=x n](θ).|Let d=(−k (c)))−1/2.Then if the distribution of U is the posterior distribution,then the distribution of U−cis dapproximately standard normal;more precisely:that distribution approaches the standard normal distribution as n approaches infinity. Important:In#1(e)on the second problem set,try both versions of the Bayesian Central Limit pare the result of each with the result alleged to be correct in the second-to-last paragraph on thefirst page of this handout.Here is a rough sketch of a proof of the second version.Denote fΘX1=x1,...,X n=x n (θ)by f(θ),|and let k(θ)=log f(θ).Thenf(θ)=e log f(θ)k(θ)k(c)+k (c)(θ−c)+k (c)(θ−c)2/2+k (c)(θ−c)3/6+=e=e·········k(c)+k (c)(θ−c)+k (c)(θ−c)2/2=e∼k(c)+0+k (c)(θ−c)2/2=e=[constant]·e k (c)(θ−c)2/2(We must have k (c)=0and k (c)<0since there is a maximum at c.)=[constant]·e−(θ−c)2/(2σ2)provided−1/σ2=k (c),i.e.,σ=(−1/k (c))−1/2,which is the same thing that followed the words“Let d=”above.Why do the3rd-and higher-degree terms evaporate as n→∞above?To answer that you need to know how they depend on n.Those interested in such theoretical questions can work out the details.2。

·1293·二斑素瓢虫线粒体基因组全序列测定和分析朱国渊1，张永科1，2*，孔祥东1，许丽月1，黎小清1，张小娇1，王进强1，段波1*（1云南省热带作物科学研究所，云南景洪666100；2云南农业大学植物保护学院，云南昆明650201）摘要：【目的】在线粒体基因组水平探讨二斑素瓢虫［Illeis bistigmosa （Mulsant ，1850）］完整的线粒体基因组结构特征，分析其在瓢虫科（Coccinellidae ）内的分类地位，为揭示瓢虫科昆虫的系统发育和进化关系提供理论依据。

【方法】通过Illumina 二代测序技术对二斑素瓢虫线粒体基因组进行测序，对基因组序列进行拼装和注释，分析其结构特点和碱基组成；使用最大似然法和贝叶斯法构建瓢虫科16个种线粒体基因组的系统发育进化树，分析其在瓢虫科内的系统发育关系。

【结果】二斑素瓢虫线粒体基因组序列全长17840bp ，包含13个蛋白质编码基因，22个tRNA 基因，2个rRNA 基因和1个A+T 富含区（A+T rich region ）。

基因组全序列的AT 含量为78.44%，表现明显的A+T 偏向性。

13个蛋白质编码基因中除cox1基因以TTG 为起始密码子外，其余12个基因均以ATN 为起始密码子；nad4、nad5、cox2和cox3基因的终止密码子以单独的T 结尾。

除trnS1和trnP 基因外，其余20个tRNAs 的二级结构均为典型的三叶草结构，二级结构中有U-U 或U-G 碱基错配现象。

系统发育分析显示，瓢虫亚科的所有物种聚在同一分支，支持该亚科的单系性；二斑素瓢虫与素菌瓢虫属（Illeis ）另外2个种［柯氏素菌瓢虫（I.koebelei ）和素鞘瓢虫（I.cincta ）］聚为一个分支，组成姐妹关系。

【结论】获得了二斑素瓢虫的线粒体基因组全序列，二斑素瓢虫线粒体基因组符合瓢虫科昆虫线粒体基因组的一般特征；二斑素瓢虫与柯氏素菌瓢虫和素鞘瓢虫的亲缘关系较近，与传统的形态学分类相一致。

Package‘sEparaTe’August18,2023Title Maximum Likelihood Estimation and Likelihood Ratio TestFunctions for Separable Variance-Covariance StructuresVersion0.3.2Maintainer Timothy Schwinghamer<***************************.CA>Description Maximum likelihood estimation of the parametersof matrix and3rd-order tensor normal distributions with unstructuredfactor variance covariance matrices,two procedures,and for unbiasedmodified likelihood ratio testing of simple and double separabilityfor variance-covariance structures,two procedures.References:Dutilleul P.(1999)<doi:10.1080/00949659908811970>,Manceur AM,Dutilleul P.(2013)<doi:10.1016/j.cam.2012.09.017>,and Manceur AM,DutilleulP.(2013)<doi:10.1016/j.spl.2012.10.020>.Depends R(>=4.3.0)License MIT+file LICENSEEncoding UTF-8LazyData trueRoxygenNote7.2.0NeedsCompilation noAuthor Ameur Manceur[aut],Timothy Schwinghamer[aut,cre],Pierre Dutilleul[aut,cph]Repository CRANDate/Publication2023-08-1807:50:02UTCR topics documented:data2d (2)data3d (2)lrt2d_svc (3)lrt3d_svc (5)mle2d_svc (7)12data3d mle3d_svc (8)sEparaTe (10)Index11 data2d Two dimensional data setDescriptionAn i.i.d.random sample of size7from a2x3matrix normal distribution,for a small numerical example of the use of the functions mle2d_svc and lrt2d_svc from the sEparaTe packageUsagedata2dFormatA frame(excluding the headings)with42lines of data and4variables:K an integer ranging from1to7,the size of an i.i.d.random sample from a2x3matrix normal distributionId1an integer ranging from1to2,the number of rows of the matrix normal distributionId2an integer ranging from1to3,the number of columns of the matrix normal distribution value2d the sample data for the observed variabledata3d Three dimensional data setDescriptionAn i.i.d.random sample of size13from a2x3x2tensor normal distribution,for a small numerical example of the use of the functions mle3d_svc and lrt3d_svc from the sEparaTe packageUsagedata3dFormatA frame(excluding the headings)with156lines of data and5variables:K an integer ranging from1to13,the size of an i.i.d.random sample from a2x3x2tensor matrix normal distributionId3an integer ranging from1to2,the number of rows of the3rd-order tensor normal distribution Id4an integer ranging from1to3,the number of columns of the3rd-order tensor normal distribu-tionId5an integer ranging from1to2,the number of edges of the3rd-order tensor normal distribution value3d the sample data for the observed variablelrt2d_svc Unbiased modified likelihood ratio test for simple separability of avariance-covariance matrix.DescriptionA likelihood ratio test(LRT)for simple separability of a variance-covariance matrix,modified tobe unbiased infinite samples.The modification is a penalty-based homothetic transformation of the LRT statistic.The penalty value is optimized for a given mean model,which is left unstruc-tured here.In the required function,the Id1and Id2variables correspond to the row and column subscripts,respectively;“value2d”refers to the observed variable.Usagelrt2d_svc(value2d,Id1,Id2,subject,data_2d,eps,maxiter,startmat,sign.level,n.simul)Argumentsvalue2d from the formula value2d~Id1+Id2Id1from the formula value2d~Id1+Id2Id2from the formula value2d~Id1+Id2subject the replicate,also called the subject or individual,thefirst column in the matrix (2d)datafiledata_2d the name of the matrix dataeps the threshold in the stopping criterion for the iterative mle algorithm(estimation) maxiter the maximum number of iterations for the mle algorithm(estimation)startmat the value of the second factor variance-covariance matrix used for initialization,i.e.,to start the mle algorithm(estimation)and obtain the initial estimate of thefirst factor variance-covariance matrixsign.level the significance level,or rejection rate in the testing of the null hypothesis of simple separability for a variance-covariance structure,when the unbiased mod-ified LRT is used,i.e.,the critical value in the chi-square test is derived bysimulations from the sampling distribution of the LRT statistic n.simul the number of simulations used to build the sampling distribution of the LRT statistic under the null hypothesis,using the same characteristics as the i.i.d.random sample from a matrix normal distributionOutput“Convergence”,TRUE or FALSE“chi.df”,the theoretical number of degrees of freedom of the asymptotic chi-square distribution that would apply to the unmodified LRT statistic for simple separability of a variance-covariance structure“Lambda”,the observed value of the unmodified LRT statistic“critical.value”,the critical value at the specified significance level for the chi-square distribution with“chi.df”degrees of freedom“mbda”will indicate whether or not the null hypothesis of separability was rejected, based on the theoretical LRT statistic“Simulation.critical.value”,the critical value at the specified significance level that is derived from the sampling distribution of the unbiased modified LRT statistic“mbda.simulation”,the decision(acceptance/rejection)regarding the null hypothesis of simple separability,made using the theoretical(biased unmodified)LRT“Penalty”,the optimized penalty value used in the homothetic transformation between the biased unmodified and unbiased modified LRT statistics“U1hat”,the estimated variance-covariance matrix for the rows“Standardized_U1hat”,the standardized estimated variance-covariance matrix for the rows;the standardization is performed by dividing each entry of U1hat by entry(1,1)of U1hat“U2hat”,the estimated variance-covariance matrix for the columns“Standardized_U2hat”,the standardized estimated variance-covariance matrix for the columns;the standardization is performed by multiplying each entry of U2hat by entry(1,1)of U1hat“Shat”,the sample variance-covariance matrix computed from the vectorized data matrices ReferencesManceur AM,Dutilleul P.2013.Unbiased modified likelihood ratio tests for simple and double separability of a variance-covariance structure.Statistics and Probability Letters83:631-636.Examplesoutput<-lrt2d_svc(data2d$value2d,data2d$Id1,data2d$Id2,data2d$K,data_2d=data2d,n.simul=100)outputlrt3d_svc An unbiased modified likelihood ratio test for double separability of avariance-covariance structure.DescriptionA likelihood ratio test(LRT)for double separability of a variance-covariance structure,modified tobe unbiased infinite samples.The modification is a penalty-based homothetic transformation of the LRT statistic.The penalty value is optimized for a given mean model,which is left unstructured here.In the required function,the Id3,Id4and Id5variables correspond to the row,column and edge subscripts,respectively;“value3d”refers to the observed variable.Usagelrt3d_svc(value3d,Id3,Id4,Id5,subject,data_3d,eps,maxiter,startmatU2,startmatU3,sign.level,n.simul)Argumentsvalue3d from the formula value3d~Id3+Id4+Id5Id3from the formula value3d~Id3+Id4+Id5Id4from the formula value3d~Id3+Id4+Id5Id5from the formula value3d~Id3+Id4+Id5subject the replicate,also called individualdata_3d the name of the tensor dataeps the threshold in the stopping criterion for the iterative mle algorithm(estimation) maxiter the maximum number of iterations for the mle algorithm(estimation)startmatU2the value of the second factor variance-covariance matrix used for initialization startmatU3the value of the third factor variance-covariance matrix used for initialization,i.e.,startmatU3together with startmatU2are used to start the mle algorithm(estimation)and obtain the initial estimate of thefirst factor variance-covariancematrix U1sign.level the significance level,or rejection rate in the testing of the null hypothesis of simple separability for a variance-covariance structure,when the unbiased mod-ified LRT is used,i.e.,the critical value in the chi-square test is derived bysimulations from the sampling distribution of the LRT statistic n.simul the number of simulations used to build the sampling distribution of the LRT statistic under the null hypothesis,using the same characteristics as the i.i.d.random sample from a tensor normal distributionOutput“Convergence”,TRUE or FALSE“chi.df”,the theoretical number of degrees of freedom of the asymptotic chi-square distribution that would apply to the unmodified LRT statistic for double separability of a variance-covariance structure“Lambda”,the observed value of the unmodified LRT statistic“critical.value”,the critical value at the specified significance level for the chi-square distribution with“chi.df”degrees of freedom“mbda”,the decision(acceptance/rejection)regarding the null hypothesis of double sep-arability,made using the theoretical(biased unmodified)LRT“Simulation.critical.value”,the critical value at the specified significance level that is derived from the sampling distribution of the unbiased modified LRT statistic“mbda.simulation”,the decision(acceptance/rejection)regarding the null hypothesis of double separability,made using the unbiased modified LRT“Penalty”,the optimized penalty value used in the homothetic transformation between the biased unmodified and unbiased modified LRT statistics“U1hat”,the estimated variance-covariance matrix for the rows“Standardized_U1hat”,the standardized estimated variance-covariance matrix for the rows;the standardization is performed by dividing each entry of U1hat by entry(1,1)of U1hat“U2hat”,the estimated variance-covariance matrix for the columns“Standardized_U2hat”,the standardized estimated variance-covariance matrix for the columns;the standardization is performed by multiplying each entry of U2hat by entry(1,1)of U1hat“U3hat”,the estimated variance-covariance matrix for the edges“Shat”,the sample variance-covariance matrix computed from the vectorized data tensorsReferencesManceur AM,Dutilleul P.2013.Unbiased modified likelihood ratio tests for simple and double separability of a variance-covariance structure.Statistics and Probability Letters83:631-636.Examplesoutput<-lrt3d_svc(data3d$value3d,data3d$Id3,data3d$Id4,data3d$Id5,data3d$K,data_3d=data3d,n.simul=100)outputmle2d_svc Maximum likelihood estimation of the parameters of a matrix normaldistributionDescriptionMaximum likelihood estimation for the parameters of a matrix normal distribution X,which is char-acterized by a simply separable variance-covariance structure.In the general case,which is the case considered here,two unstructured factor variance-covariance matrices determine the covariability of random matrix entries,depending on the row(one factor matrix)and the column(the other factor matrix)where two X-entries are.In the required function,the Id1and Id2variables correspond to the row and column subscripts,respectively;“value2d”indicates the observed variable.Usagemle2d_svc(value2d,Id1,Id2,subject,data_2d,eps,maxiter,startmat)Argumentsvalue2d from the formula value2d~Id1+Id2Id1from the formula value2d~Id1+Id2Id2from the formula value2d~Id1+Id2subject the replicate,also called individualdata_2d the name of the matrix dataeps the threshold in the stopping criterion for the iterative mle algorithmmaxiter the maximum number of iterations for the iterative mle algorithmstartmat the value of the second factor variance-covariance matrix used for initializa-tion,i.e.,to start the algorithm and obtain the initial estimate of thefirst factorvariance-covariance matrixOutput“Convergence”,TRUE or FALSE“Iter”,will indicate the number of iterations needed for the mle algorithm to converge“Xmeanhat”,the estimated mean matrix(i.e.,the sample mean)“First”,the row subscript,or the second column in the datafile“U1hat”,the estimated variance-covariance matrix for the rows“Standardized.U1hat”,the standardized estimated variance-covariance matrix for the rows;the stan-dardization is performed by dividing each entry of U1hat by entry(1,1)of U1hat“Second”,the column subscript,or the third column in the datafile“U2hat”,the estimated variance-covariance matrix for the columns“Standardized.U2hat”,the standardized estimated variance-covariance matrix for the columns;the standardization is performed by multiplying each entry of U2hat by entry(1,1)of U1hat“Shat”,is the sample variance-covariance matrix computed from of the vectorized data matrices ReferencesDutilleul P.1990.Apport en analyse spectrale d’un periodogramme modifie et modelisation des series chronologiques avec repetitions en vue de leur comparaison en frequence.D.Sc.Dissertation, Universite catholique de Louvain,Departement de mathematique.Dutilleul P.1999.The mle algorithm for the matrix normal distribution.Journal of Statistical Computation and Simulation64:105-123.Examplesoutput<-mle2d_svc(data2d$value2d,data2d$Id1,data2d$Id2,data2d$K,data_2d=data2d) outputmle3d_svc Maximum likelihood estimation of the parameters of a3rd-order ten-sor normal distributionDescriptionMaximum likelihood estimation for the parameters of a3rd-order tensor normal distribution X, which is characterized by a doubly separable variance-covariance structure.In the general case, which is the case considered here,three unstructured factor variance-covariance matrices determine the covariability of random tensor entries,depending on the row(one factor matrix),the column (another factor matrix)and the edge(remaining factor matrix)where two X-entries are.In the required function,the Id3,Id4and Id5variables correspond to the row,column and edge subscripts, respectively;“value3d”indicates the observed variable.Usagemle3d_svc(value3d,Id3,Id4,Id5,subject,data_3d,eps,maxiter,startmatU2,startmatU3)Argumentsvalue3d from the formula value3d~Id3+Id4+Id5Id3from the formula value3d~Id3+Id4+Id5Id4from the formula value3d~Id3+Id4+Id5Id5from the formula value3d~Id3+Id4+Id5subject the replicate,also called individualdata_3d the name of the tensor dataeps the threshold in the stopping criterion for the iterative mle algorithmmaxiter the maximum number of iterations for the iterative mle algorithmstartmatU2the value of the second factor variance covariance matrix used for initialization startmatU3the value of the third factor variance covariance matrix used for initialization,i.e.,startmatU3together with startmatU2are used to start the algorithm andobtain the initial estimate of thefirst factor variance covariance matrix U1Output“Convergence”,TRUE or FALSE“Iter”,the number of iterations needed for the mle algorithm to converge“Xmeanhat”,the estimated mean tensor(i.e.,the sample mean)“First”,the row subscript,or the second column in the datafile“U1hat”,the estimated variance-covariance matrix for the rows“Standardized.U1hat”,the standardized estimated variance-covariance matrix for the rows;the stan-dardization is performed by dividing each entry of U1hat by entry(1,1)of U1hat“Second”,the column subscript,or the third column in the datafile“U2hat”,the estimated variance-covariance matrix for the columns“Standardized.U2hat”,the standardized estimated variance-covariance matrix for the columns;the standardization is performed by multiplying each entry of U2hat by entry(1,1)of U1hat“Third”,the edge subscript,or the fourth column in the datafile“U3hat”,the estimated variance-covariance matrix for the edges“Shat”,the sample variance-covariance matrix computed from the vectorized data tensorsReferenceManceur AM,Dutilleul P.2013.Maximum likelihood estimation for the tensor normal distribution: Algorithm,minimum sample size,and empirical bias and dispersion.Journal of Computational and Applied Mathematics239:37-49.10sEparaTeExamplesoutput<-mle3d_svc(data3d$value3d,data3d$Id3,data3d$Id4,data3d$Id5,data3d$K,data_3d=data3d) outputsEparaTe MLE and LRT functions for separable variance-covariance structuresDescriptionA package for maximum likelihood estimation(MLE)of the parameters of matrix and3rd-ordertensor normal distributions with unstructured factor variance-covariance matrices(two procedures),and for unbiased modified likelihood ratio testing(LRT)of simple and double separability forvariance-covariance structures(two procedures).Functionsmle2d_svc,for maximum likelihood estimation of the parameters of a matrix normal distributionmle3d_svc,for maximum likelihood estimation of the parameters of a3rd-order tensor normaldistributionlrt2d_svc,for the unbiased modified likelihood ratio test of simple separability for a variance-covariance structurelrt3d_svc,for the unbiased modified likelihood ratio test of double separability for a variance-covariance structureDatadata2d,a two-dimensional data setdata3d,a three-dimensional data setReferencesDutilleul P.1999.The mle algorithm for the matrix normal distribution.Journal of StatisticalComputation and Simulation64:105-123.Manceur AM,Dutilleul P.2013.Maximum likelihood estimation for the tensor normal distribution:Algorithm,minimum sample size,and empirical bias and dispersion.Journal of Computational andApplied Mathematics239:37-49.Manceur AM,Dutilleul P.2013.Unbiased modified likelihood ratio tests for simple and doubleseparability of a variance covariance structure.Statistics and Probability Letters83:631-636.Index∗datasetsdata2d,2data3d,2data2d,2data3d,2lrt2d_svc,3lrt3d_svc,5mle2d_svc,7mle3d_svc,8sEparaTe,1011。

Package‘BayesMixSurv’October12,2022Type PackageTitle Bayesian Mixture Survival Models using AdditiveMixture-of-Weibull Hazards,with Lasso Shrinkage andStratificationVersion0.9.1Date2016-09-08Author Alireza S.Mahani,Mansour T.A.SharabianiMaintainer Alireza S.Mahani<**************************>Description Bayesian Mixture Survival Models using Additive Mixture-of-Weibull Hazards,with Lasso Shrinkage andStratification.As a Bayesian dynamic survival model,it relaxes the proportional-hazard sso shrinkage controlsoverfitting,given the increase in the number of free parameters in the model due to pres-ence of two Weibull componentsin the hazard function.License GPL(>=2)Depends survivalNeedsCompilation noRepository CRANDate/Publication2016-09-0810:24:27R topics documented:bayesmixsurv (2)bayesmixsurv.crossval (5)plot.bayesmixsurv (7)predict.bayesmixsurv (8)summary.bayesmixsurv (10)Index121bayesmixsurv Dynamic Bayesian survival model-with stratification and Lassoshrinkage-for right-censored data using two-component additivemixture-of-Weibull hazards.DescriptionBayesian survival model for right-censored data,using a sum of two hazard functions,each hav-ing a power dependence on time,corresponding to a Weibull distribution on event density.(Note that event density function for the mixture model does NOT remain a Weibull distribution.)Each component has a different shape and scale parameter,with scale parameters each being the ex-ponential of a linear function of covariates specified in formula1and formula2.Stratification is implemented using a common set of intercepts between the two sso shrinkage-using Laplace prior on coefficients(Park and Casella2008)-allows for variable selection in the presence of low observation-to-variable ratio.The mixture model allows for time-dependent(and context-dependent)hazard ratios.Confidence intervals for coefficient estimation and prediction are generated using full Bayesian paradigm,i.e.by keeping all samples rather than summarizing them into mean and sd.Posterior distribution is estimated via MCMC sampling,using univariate slice sampler with stepout and shrinkage(Neal2003).Usagebayesmixsurv(formula1,data,formula2=formula1,stratCol=NULL,weights,subset,na.action=na.fail,control=bayesmixsurv.control(),print.level=2)bayesmixsurv.control(single=FALSE,alpha2.fixed=NULL,alpha.boundary=1.0,lambda1=1.0 ,lambda2=lambda1,iter=1000,burnin=round(iter/2),sd.thresh=1e-4,scalex=TRUE,nskip=round(iter/10))##S3method for class bayesmixsurvprint(x,...)Argumentsformula1Survival formula expressing the time/status variables as well as covariates usedin thefirst component.data Data frame containing the covariates and response variable,as well as the strat-ification column.formula2Survival formula expressing the covariates used in the second component.Noleft-hand side is necessary since the response variable information is extractedfrom formula1.Defaults to formula1.stratCol Name of column in data used for stratification.Must be a factor or coerced intoone.Default is no stratification(stratCol=NULL).weights Optional vector of case weights.*Not supported yet*subset Subset of the observations to be used in thefit.*Not supported yet*na.action Missing-datafilter function.*Not supported yet(only na.fail behavior works)*control See bayesmixsurv.control for a description of the parameters inside the control list.print.level Controlling verbosity level.single If TRUE,a single-component model,equivalent to Bayesian Weibull survival re-gression,with Lasso shrinkage,is implemented.Default is FALSE,i.e.a two-component mixture-of-Weibull model.alpha2.fixed If provided,it specifies the shape parameter of the second component.Defaultis NULL,which allows the MCMC sampling to estimate both shape parameters.alpha.boundary When single=FALSE and alpha2.fixed=NULL,this parameter specifies an up-per bound for the shape parameter of thefirst component,and a lower bound forthe shape parameter of the second component.These boundary conditions areenforced in the univariate slice sampler function calls.lambda1Lasso Shrinkage parameter used in the Laplace prior on covariates used in thefirst component.lambda2Lasso Shrinkage parameter used in the Laplace prior on covariates used in thesecond component.Defaults to lambda1.iter Number of posterior MCMC samples to generate.burnin Number of initial MCMC samples to discard before calculating summary statis-tics.sd.thresh Threshold for standard deviation of a covariate(after possible centering/scaling).If below the threshold,the corresponding coefficient is removed from sampling,i.e.its value is clamped to zero.scalex If TRUE,each covariate vector is centered and scaled before model estimation.The scaling parameters are saved in return object,and used in subsequent callsto predict ers are strongly advised against turning this feature off,since the quality of Gibbs sampling MCMC is greatly enhanced by covariatecentering and scaling.nskip Controlling how often to print progress report during MCMC run.For example,if nskip=10,progress will be reported after10,20,30,...samples.x Object of class’bayesmixsurv’,usually the result of a call to bayesmixsurv....Arguments to be passed to/from other methods.ValueThe function bayesmixsurv.control return a list with the same elements as its input parameters.The function bayesmixsurv returns object of class bayesmixsurv,with the following components: call The matched callformula1Same as input.formula2Same as input.weights Same as input.*Not supported yet*subset Same as input.*Not supported yet*na.action Same as input.*Not supported yet*(current behavior is na.fail)control Same as input.X1Model matrix used for component1,after potential centering and scaling.X2Model matrix used for component2,after potential centering and scaling.y Survival response variable(time and status)used in the model.contrasts1The contrasts used for component1(where relevant).contrasts2The contrasts used for component2(where relevant).xlevels1A record of the levels of the factors used infitting for component1(where relevant).xlevels2A record of the levels of the factors used infitting for component2(where relevant).terms1The terms object used for component1.terms2The terms object used for component2.colnamesX1Names of columns for X1,also names of scale coefficients for component1. colnamesX2Names of columns for X1,also names of scale coefficients for component2. apply.scale.X1Index of columns of X1where scaling has been applied.apply.scale.X2Index of columns of X2where scaling has been applied.centerVec.X1Vector of centering parameters for columns of X1indicated by apply.scale.X1. centerVec.X2Vector of centering parameters for columns of X2indicated by apply.scale.X2. scaleVec.X1Vector of scaling parameters for columns of X1indicated by apply.scale.X1. scaleVec.X2Vector of scaling parameters for columns of X2indicated by apply.scale.X2. Xg Model matrix associated with stratification(if any).stratContrasts The contrasts used for stratification model matrix,if any.stratXlevels A record of the levels of the factors used in stratification(if any)). stratTerms The terms object used for stratification.colnamesXg Names of columns for Xg.idx1Vector of indexes into X1for which sampling occured.All columns of X1whose standard deviation falls below sd.thresh are excluded from sampling and theircorresponding coefficients are clamped to0.idx2Vector of indexes into X2for which sampling occured.All columns of X2whose standard deviation falls below sd.thresh are excluded from sampling and theircorresponding coefficients are clamped to0.median List of median values,with elements including alpha1,alpha2(shape param-eter of components1and2),beta1,beta2(coefficients of scale parameter forcomponents1and2),gamma(stratification intercept adjustments,shared by2comoponents),and sigma.gamma(standard deviation of zero-mean Gaussiandistribution that is the prior for gamma’s).max Currently,a list with one element,loglike,containing the maximum sampled log-likelihood of the model.smp List of coefficient samples,with elements alpha1,alpha2(shape parametersfor components1and2),beta1,beta2(scale parameter coefficients for com-ponents1and2),loglike(model log-likelihood),gamma(stratification interceptadjustments,shared by2comoponents),and sigma.gamma(standard deviationof zero-mean Gaussian distribution that is the prior for gamma’s).Each param-eter has iter samples.For vector parameters,first dimension is the number ofsamples(iter),while the second dimension is the length of the vector.Author(s)Alireza S.Mahani,Mansour T.A.SharabianiReferencesNeal R.M.(2003).Slice Sampling.Annals of Statistics,31,705-767.Park T.and Casella G.(2008)The Bayesian Lasso.Journal of the American Statistical Association,103,681-686.Examples#NOTE:to ensure convergence,typically more than100samples are needed#fit the most general model,with two Weibull components and unspecified shape parametersret<-bayesmixsurv(Surv(time,status)~as.factor(trt)+age+as.factor(celltype)+prior,veteran,control=bayesmixsurv.control(iter=100))#fix one of the two shape parametersret2<-bayesmixsurv(Surv(time,status)~as.factor(trt)+age+as.factor(celltype)+prior,veteran ,control=bayesmixsurv.control(iter=100,alpha2.fixed=1.0))bayesmixsurv.crossval Convenience functions for cross-validation-based selection of shrink-age parameter in the bayesmixsurv model.Descriptionbayesmixsurv.crossval calculates cross-validation-based,out-of-sample log-likelihood of a bsgwmodel for a data set,given the supplied folds.bayesmixsurv.crossval.wrapper applies bayesmixsurv.crossval to a set of combinations of shrinkage parameters(lambda1,lambda2)and produces the resultingvector of log-likelihood values as well as the specific combination of shrinkage parameters asso-ciated with the maximum log-likelihood.bayesmixsurv.generate.folds generates random par-titions,while bayesmixsurv.generate.folds.eventbalanced generates random partitions withevents evenly distributed across partitions.The latter feature is useful for cross-valiation of smalldata sets with low event rates,since it prevents over-accumulation of events in one or two partitions,and lack of events altogether in other partitions.Usagebayesmixsurv.generate.folds(ntot,nfold=5)bayesmixsurv.generate.folds.eventbalanced(formula,data,nfold=5)bayesmixsurv.crossval(data,folds,all=FALSE,print.level=1,control=bayesmixsurv.control(),...)bayesmixsurv.crossval.wrapper(data,folds,all=FALSE,print.level=1,control=bayesmixsurv.control(),lambda.min=0.01,lambda.max=100,nlambda=10,lambda1.vec=exp(seq(from=log(lambda.min),to=log(lambda.max),length.out=nlambda)) ,lambda2.vec=NULL,lambda12=if(is.null(lambda2.vec))cbind(lambda1=lambda1.vec,lambda2=lambda1.vec)else as.matrix(expand.grid(lambda1=lambda1.vec,lambda2=lambda2.vec)),plot=TRUE,...) Argumentsntot Number of observations to create partitions for.It must typically be set tonrow(data).nfold Number of folds or partitions to generate.formula Formula specifying the covariates to be used in component1,and the time/statusresponse variable in the survival model.data Data frame containing the covariates and response,used in training and predic-tion.folds An integer vector of length nrow(data),defining fold/partition membershipof each observation.For example,in5-fold cross-validation for a data set of200observations,folds must be a200-long vector with elements from theset{1,2,3,4,5}.Convenience functions bayesmixsurv.generate.folds andbayesmixsurv.generate.folds.eventbalanced can be used to generate thefolds vector for a given survival data frame.all If TRUE,estimation objects from each cross-validation task is collected and re-turned for diagnostics purposes.print.level Verbosity of progress report.control List of control parameters,usually the output of bayesmixsurv.control.lambda.min Minimum value used to generate lambda.vec.lambda.max Maximum value used to generate lambda.vec.nlambda Length of lambda.vec vector.lambda1.vec Vector of shrinkage parameters to be tested for component-1coefficients.lambda2.vec Vector of shrinkage parameters to be tested for component-2coefficients.lambda12A data frame that enumerates all combinations of lambda1and lambda2to betested.By default,it is constructed from forming all permutations of lambda1.vecand lambda2.vec.If lambda2.vec=NULL,it will only try equal values of the twoparameters in each combination.plot If TRUE,and if the lambda1and lambda2entries in lambda12are identical,aplot of loglike as a function of either vector is produced....Further arguments passed to bayesmixsurv.plot.bayesmixsurv7ValueFunctions bayesmixsurv.generate.folds and bayesmixsurv.generate.folds.eventbalanced produce integer vectors of length ntot or nrow(data)respectively.The output of these functionscan be directly passed to bayesmixsurv.crossval or bayesmixsurv.crossval.wrapper.Func-tion bayesmixsurv.crossval returns the log-likelihood of data under the assumed bsgw model,calculated using a cross-validation scheme with the supplied fold parameter.If all=TRUE,the esti-mation objects for each of the nfold estimation jobs will be returned as the"estobjs"attribute of thereturned value.Function bayesmixsurv.crossval.wrapper returns a list with elements lambda1and lambda2,the optimal shrinkage parameters for components1and2,respectively.Additionally,the following attributes are attached:loglike.vec Vector of log-likelihood values,one for each tested combination of lambda1andlambda2.loglike.opt The maximum log-likelihood value from the loglike.vec.lambda12Data frame with columns lambda1and lambda2.Each row of this data framecontains one combination of shrinkage parameters that are tested in the wrapperfunction.estobjs If all=TRUE,a list of length nrow(lambda12)is returned,with each elementbeing itself a list of nfold estimation objects associated with each call to thebayesmixsurv function.This object can be examined by the user for diagnosticpurposes,e.g.by applying plot against each object.Author(s)Alireza S.Mahani,Mansour T.A.SharabianiExamples#NOTE:to ensure convergence,typically more than30samples are neededfolds<-bayesmixsurv.generate.folds.eventbalanced(Surv(futime,fustat)~1,ovarian,5)cv<-bayesmixsurv.crossval(ovarian,folds,formula1=Surv(futime,fustat)~ecog.ps+rx,control=bayesmixsurv.control(iter=30,nskip=10),print.level=3)cv2<-bayesmixsurv.crossval.wrapper(ovarian,folds,formula1=Surv(futime,fustat)~ecog.ps+rx ,control=bayesmixsurv.control(iter=30,nskip=10),lambda1.vec=exp(seq(from=log(0.1),to=log(1),length.out=3)))plot.bayesmixsurv Plot diagnostics for a bayesmixsurv objectDescriptionFour sets of MCMC diagnostic plots are currently generated:1)log-likelihood trace plots,2)coef-ficient trace plots,3)coefficient autocorrelation plots,4)coefficient histograms.Usage##S3method for class bayesmixsurvplot(x,pval=0.05,burnin=round(x$control$iter/2),nrow=2,ncol=3,...)Argumentsx A bayesmixsurv object,typically the output of bayesmixsurv function.pval The P-value at which lower/upper bounds on coefficients are calculated andoverlaid on trace plots and historgrams.burnin Number of samples discarded from the beginning of an MCMC chain,afterwhich parameter quantiles are calculated.nrow Number of rows of subplots within eachfigure,applied to plot sets2-4.ncol Number of columns of subplots within eachfigure,applied to plot sets2-4....Further arguments to be passed to/from other methods.Author(s)Alireza S.Mahani,Mansour T.A.SharabianiExamplesest<-bayesmixsurv(Surv(futime,fustat)~ecog.ps+rx,ovarian,control=bayesmixsurv.control(iter=800,nskip=100))plot(est)predict.bayesmixsurv Predict method for bayesmixsurv modelfitsDescriptionCalculates log-likelihood and hazard/cumulative hazard/survival functions over a user-supplied vec-tor time values,based on bayesmixsurv model object.Usage##S3method for class bayesmixsurvpredict(object,newdata=NULL,tvec=NULL,burnin=object$control$burnin,...)##S3method for class predict.bayesmixsurvsummary(object,idx=1:dim(object$smp$h)[3],burnin=object$burnin,pval=0.05,popmean=identical(idx,1:dim(object$smp$h)[3]),make.plot=TRUE,...)Argumentsobject For predict.bayesmixsurv,an object of class"bayesmixsurv",usually the re-sult of a call to bayesmixsurv;for summary.predict.bayesmixsurv,an objectof class"predict.bayesmixsurv",usually the result of a call to predict.bayesmixsurv.newdata An optional data frame in which to look for variables with which to predict.Ifomiited,thefitted values(training set)are used.tvec An optional vector of time values,along which time-dependent entities(haz-ard,cumulative hazard,survival)will be predicted.If omitted,only the time-independent entities(currently only log-likelihood)will be calculated.If a singleinteger is provided for tvec,it is interpreted as number of time points,equallyspaced from0to object$tmax:tvec<-seq(from=0.0,to=object$tmax,length.out=tvec).burnin Number of samples to discard from the beginning of each MCMC chain beforecalculating median value(s)for time-independent entities.idx Index of observations(rows of newdata or training data)for which to generatesummary statistics.Default is the entire data.pval Desired p-value,based on which lower/upper bounds will be calculated.Defaultis0.05.popmean Whether population averages must be calculated or not.By default,populationaverages are only calculated when the entire data is included in prediction.make.plot Whether population mean and other plots must be created or not....Further arguments to be passed to/from other methods.DetailsThe time-dependent predicted objects(except loglike)are three-dimensional arrays of size(nsmpx nt x nobs),where nsmp=number of MCMC samples,nt=number of time values in tvec,andnobs=number of rows in newdata.Therefore,even for modest data sizes,these objects can occupylarge chunks of memory.For example,for nsmp=1000,nt=100,nobs=1000,the three objects h,H,S have a total size of2.2GB.Since applying quantile to these arrays is time-consuming(asneeded for calculation of median and lower/upper bounds),we have left such summaries out ofthe scope of predict ers can instead apply summary to the prediction object to obtainsummary statistics.During cross-validation-based selection of shrinkage parameter lambda,thereis no need to supply tvec since we only need the log-likelihood value.This significantly speeds upthe parameter-tuning process.The function summary.predict.bayesmixsurv allows the user tocalculates summary statistics for a subset(or all of)data,if desired.This approach is in line with theoverall philosophy of delaying the data summarization until necessary,to avoid unnecessary loss inaccuracy due to premature blending of information contained in individual samples.ValueThe function predict.bayesmixsurv returns as object of class"predict.bayesmixsurv"with thefollowingfields:tvec Actual vector of time values(if any)used for prediction.burnin Same as input.median List of median values for predicted entities.Currently,only loglike is pro-duced.See’Details’for explanation.smp List of MCMC samples for predicted entities.Elements include h1,h2,h(haz-ard functions for components1,2and their sum),H1,H2,H(cumulative hazardfunctions for components1,2and their sum),S(survival function),and loglike(model log-likelihood).All functions are evaluated over time values specified intvec.10summary.bayesmixsurvkm.fit Kaplan-Meyerfit of the data used for prediction(if data contains responsefields).The function summary.predict.bayesmixsurv returns a list with the followingfields:lower A list of lower-bound values for h,H,S,hr(hazard ratio of idx[2]to idx[1]observation),and S.diff(survival probability of idx[2]minus idx[1]).Thelast two are only included if length(idx)==2.median List of median values for same entities described in lower.upper List of upper-bound values for same entities described in lower.popmean Lower-bound/median/upper-bound values for population average of survival prob-ability.km.fit Kaplan-Meyerfit associated with the prediction object(if available).Author(s)Alireza S.Mahani,Mansour T.A.SharabianiExamplesest<-bayesmixsurv(Surv(futime,fustat)~ecog.ps+rx+age,ovarian,control=bayesmixsurv.control(iter=400,nskip=100))pred<-predict(est,tvec=50)predsumm<-summary(pred,idx=1:10)summary.bayesmixsurv Summarizing BayesMixSurv modelfitsDescriptionsummary method for class"bayesmixsurv".Usage##S3method for class bayesmixsurvsummary(object,pval=0.05,burnin=object$control$burnin,...)##S3method for class summary.bayesmixsurvprint(x,...)Argumentsobject An object of class’bayesmixsurv’,usually the result of a call to bayesmixsurv.x An object of class"summary.bayesmixsurv",usually the result of a call to summary.bayesmixsurv.pval Desired p-value,based on which lower/upper bounds will be calculated.Defaultis0.05.burnin Number of samples to discard from the beginning of each MCMC chain beforecalculating median and lower/upper bounds....Further arguments to be passed to/from other methods.summary.bayesmixsurv11ValueAn object of class summary.bayesmixsurv,with the following elements:call The matched call.pval Same as input.burnin Same as input.single Copied from object$control$single.See bayesmixsurv.control for explana-tion.coefficients A list including matrices alpha,beta1,beta2,and gamma(if stratification is used).Each matrix has columns named’Estimate’,’Lower Bound’,’UpperBound’,and’P-val’.alpha has two rows,one for each components,while eachof beta1and beta2has one row per covariate.gamma has one row per stratum(except for the reference stratum).Author(s)Alireza S.Mahani,Mansour T.A.SharabianiSee AlsoSee summary for a description of the generic method.The modelfitting function is bayesmixsurv.Examplesest<-bayesmixsurv(Surv(futime,fustat)~ecog.ps+rx,ovarian,control=bayesmixsurv.control(iter=800,nskip=100))summary(est,pval=0.1)Indexbayesmixsurv,2,8,10,11 bayesmixsurv.control,6,11 bayesmixsurv.crossval,5 bayesmixsurv.generate.folds(bayesmixsurv.crossval),5plot.bayesmixsurv,7predict.bayesmixsurv,8print.bayesmixsurv(bayesmixsurv),2print.summary.bayesmixsurv(summary.bayesmixsurv),10 summary,11summary.bayesmixsurv,10summary.predict.bayesmixsurv(predict.bayesmixsurv),812。

基于似然比的英文笔记检验方法The assessment of English writing proficiency has long been a critical aspect of language education and evaluation. Traditional methods, such as holistic scoring or rubric-based assessments, have their limitations in providing comprehensive and objective feedback to learners. In recent years, the development of computational linguistics and natural language processing has paved the way for more advanced techniques in writing evaluation. One such approach is the likelihood ratio-based English writing evaluation method, which offers a robust and data-driven way to assess the quality of written texts.The likelihood ratio-based method relies on the fundamental principles of statistical inference and hypothesis testing. The underlying premise is that well-written English texts exhibit certain linguistic patterns and characteristics that can be quantified and compared to a reference corpus of high-quality writing. By calculating the likelihood ratio between the target text and the reference corpus, we can determine the degree to which the target text conforms to the expected norms of proficient English writing.The likelihood ratio-based method involves several key steps. First, the target text is preprocessed to extract a set of linguistic features, such as lexical diversity, sentence complexity, and grammatical accuracy. These features are carefully selected to capture the nuances of effective written communication in English. Next, a reference corpus of high-quality English writing, typically comprising texts from professional publications or academic sources, is established. This corpus serves as the benchmark against which the target text will be evaluated.The likelihood ratio is then calculated by comparing the distribution of the linguistic features in the target text to the distribution of the same features in the reference corpus. The likelihood ratio is a statistical measure that quantifies the relative likelihood of the target text belonging to the reference corpus versus belonging to a hypothetical corpus of poorly written English. A likelihood ratio greater than 1 indicates that the target text is more likely to be from the reference corpus, suggesting a higher level of writing proficiency, while a likelihood ratio less than 1 suggests that the target text is more likely to be from a corpus of poor writing.One of the key advantages of the likelihood ratio-based method is its ability to provide detailed and actionable feedback to learners. By analyzing the specific linguistic features that contribute to thelikelihood ratio, instructors can identify the strengths and weaknesses of a student's writing and provide targeted feedback for improvement. For example, if the likelihood ratio reveals that the target text has a lower-than-expected lexical diversity, the instructor can suggest strategies for expanding the writer's vocabulary and using more varied word choices.Moreover, the likelihood ratio-based method is highly scalable and can be applied to a large volume of written texts, making it particularly useful in educational settings with large class sizes or high-stakes writing assessments. The automated nature of the analysis also ensures consistency and objectivity, reducing the potential for subjective biases that can arise in manual scoring.However, it is important to note that the likelihood ratio-based method should not be seen as a replacement for more holistic assessments of writing quality. While it provides valuable quantitative insights, it does not fully capture the nuances of effective communication, such as the coherence of ideas, the persuasiveness of arguments, or the creativity of expression. Therefore, the likelihood ratio-based method is best used in conjunction with other assessment techniques, such as rubric-based scoring or peer review, to provide a comprehensive evaluation of a writer's abilities.In conclusion, the likelihood ratio-based English writing evaluation method offers a promising approach to assessing the quality of written texts. By leveraging the power of statistical inference and computational linguistics, this method can provide detailed and objective feedback to learners, enabling them to identify and address their writing weaknesses more effectively. As the field of language assessment continues to evolve, the likelihood ratio-based method is likely to play an increasingly important role in supporting the development of English writing proficiency.。

Expert elicitation of recharge model probabilities for the Death Valley regional flow systemMing Yea,*,Karl F.Pohlmann b ,Jenny B.ChapmanbaSchool of Computational Science and Department of Geologic Sciences,Florida State University,Tallahassee,FL 32306,USA bDesert Research Institute,Nevada System of Higher Education,755East Flamingo Road,Las Vegas,NV 89119,USA Received 12January 2008;received in revised form 29February 2008;accepted 3March 2008KEYWORDSModel uncertainty;Prior model probability;Model averaging;Expert elicitation;Recharge estimates;Death Valley regional flow systemSummaryThis study uses expert elicitation to evaluate and select five alternativerecharge models developed for the Death Valley regional flow system (DVRFS),covering southeast Nevada and the Death Valley area of California,USA.The five models were developed based on three independent techniques:an empirical approach,an approach based on unsaturated-zone studies and an approach based on saturated-zone studies.It is uncertain which recharge model (or models)should be used as input for groundwater models simulating flow and contaminant transport within the DVRFS.An expert elicitation was used to evaluate and select the recharge models and to determine prior model prob-abilities used for assessing model uncertainty.The probabilities were aggregated using simple averaging and iterative methods,with the latter method also considering between-expert variability.The most favorable model,on average,is the most compli-cated model that comprehensively incorporates processes controlling net infiltration and potential recharge.The simplest model,and the most widely used,received the sec-ond highest prior probability.The aggregated prior probabilities are close to the neutral choice that treats the five models as equally likely.Thus,there is no support for selecting a single model and discarding others,based on prior information and expert judgment.This reflects the inherent uncertainty in the recharge models.If a set of prior probability from a single expert is of more interest,we suggest selecting the set of the minimum Shannon’s entropy.The minimum entropy implies the smallest amount of uncertainty and the largest amount of information used to evaluate the models.However,when enough data are available,we prefer to use a cross-validation method to select the best set of prior model probabilities that gives the best predictive performance.ª2008Elsevier B.V.All rights reserved.0022-1694/$-see front matter ª2008Elsevier B.V.All rights reserved.doi:10.1016/j.jhydrol.2008.03.001*Corresponding author.Tel.:+18506444587.E-mail address:mingye@ (M.Ye).Journal of Hydrology (2008)354,102–115a v a i l ab l e a t w w w.sc i e n c ed i re c t.c o mjou rnal homep age:www.elsevier.c om/locate/jhydro lIntroductionUncertainty analysis of hydrologic models is an essential element for decision-making in water resource manage-ment.This paper is focused on conceptual model uncer-tainty,which arises when multiple conceptualizations of a hydrologic system(or its processes)are all acceptable given available knowledge and data.A model averaging concept has been developed to assess the conceptual model uncer-tainty by averaging predictions of multiple models using appropriate weights associated with each model.The weights can be calculated using likelihood functions(Beven, 2006and its references therein)in the chi-square sense,the information criterion of AIC(Akaike,1974)or AICc(Hurvich and Tsai,1989)in the Kullback–Leibler sense(Burnham and Anderson,2002,2004;Poeter and Anderson,2005),or the information criterion of BIC(Schwarz,1978)or KIC(Kash-yap,1982)in the Bayesian sense(Draper,1995;Hoeting et al.,1999;Neuman,2003;Ye et al.,2004,2005,2008; Vrugt et al.,2006;Vrugt and Robinson,2007).This paper ad-dresses conceptual model uncertainty and model averaging in the Bayesian context.In Bayesian model averaging(BMA)(Draper,1995;Hoet-ing et al.,1999)or its maximum likelihood version(MLBMA) (Neuman,2003),if D is a quantity that one wants to predict, then its posterior distribution given conditioning data D (including measurements of model parameters and observa-tions of state variables)is the average of the distributions p(D|M k,D)under each model M k weighted by the posterior model probability p(M k|D),i.e.,pðD j DÞ¼X Kk¼1pðD j M k;DÞpðM k j DÞð1ÞThe posterior model probability,p(M k|D),is estimated via the Bayes’theorempðM k j DÞ¼pðD j M kÞpðM kÞP Kl¼1pðD j M lÞpðM lÞð2Þwhere p(D|M k)is the model likelihood function and can be approximated by p(D|M k)=exp(ÀKIC k/2)or p(D|M k)= exp(ÀBIC k/2)(Ye et al.,2004),and p(M k)is prior probability of model M k.Summation of the prior probabilities of all the alternative models is one,X Kk¼1pðM kÞ¼1ð3Þimplying that all possible models of potential relevance to the problem at hand are under study,and that all models differ from each other sufficiently to be considered mutu-ally exclusive(the joint probability of two or more models being zero).The question of how to assign prior probabili-ties p(M k)to models M k remains largely open.A common practice is to adopt a‘‘reasonable‘neutral’choice’’(Hoet-ing et al.,1999),according to which all models are initially considered to be equally likely,there being insufficient prior reason to prefer one over another.However,the neu-tral choice of prior model probabilities ignores expert knowledge of the system to be modeled,thereby implying maximum ignorance on the part of the analyst.Generally speaking,the prior model probability is an analyst’s(or a group of analysts’)subjective degree of reasonable belief(Jeffreys,1957)or confidence(Zio and Apostolakis,1996)in a model.The belief or confidence is ideally based on expert ing expert judgments is prevalent in uncertainty and risk analysis(Cooke,1991; Ayyub,2001;Bedford et al.,2006),especially when experi-mental and statistical evidence is insufficient(Refsgaard et al.,2006).For a complicated hydrologic system,expert judgment or experience is the basis of conceptual model development,and may be more informative than limited observations.This is particularly true for subsurface hydrol-ogy,where hydraulic parameters are measured from sparse samples(boreholes)and mathematical models may disagree with geologic rules(Wingle and Poeter,1993;Lele and Das, 2000).Garthwaite et al.(2005)argue that a better use of expert judgment could add more information than slight improvement of data analysis techniques.Hence,we view integrating expert judgment in BMA(by specifying subjective prior probabilities)to be a strength rather than a weakness.Madigan et al.(1995)and Zio and Apostolakis(1996)demonstrated that using informative prior model probabilities(in contrast to equal ones)on the basis of expert judgment can improve model simulation and uncertainty assessment.Ye et al.(2005)developed a constrained maximum entropy method,which estimates informative prior model probabilities through the maximiza-tion of the Shannon’s entropy(Shannon,1948)subject to constraints reflecting a single analyst’s(or group of ana-lysts’)prior perception about how plausible each alterna-tive model(or a group of models)is relative to others, and selection of the most likely among such maxima corre-sponding to alternative perceptions of various analysts(or groups of analysts).By running cross-validation,Ye et al. (2005)demonstrated that,in comparison to using equal prior model probabilities,using informative probabilities improves model predictive performance.The subjective prior model probabilities can be directly obtained through expert elicitation.The expert elicitation has been applied to many studies,for example,future cli-mate change(Arnell et al.,2005;Miklas et al.,1995),perfor-mance assessment of proposed nuclear waste repositories (Hora and Jensen,2005;McKenna et al.,2003;Draper et al.,1999;Hora and von Winterfeldt,1997;Zio and Apos-tolakis,1996;Morgan and Keith,1995;DeWispelare et al., 1995;Bonano and Apostolakis,1991;Bonano et al.,1990), estimation of parameter distributions(Parent and Bernier, 2003;Geomatrix Consultants,1998;O’Hagan,1998),devel-opment of Bayesian network(Pike,2004;Stiber et al., 1999,2004;Ghabayen et al.,2006),and interpretation of seismic images(Bond et al.,2007).Formal expert elicitation processes have been proposed by Hora and Iman(1989)and Keeney and von Winterfeldt(1991),among others.Although expert elicitation is criticized in various aspects,such as selection of experts and accurate expression of experts’knowledge and belief in probability forms(O’Hagan and Oak-ley,2004),the quality of educing expert judgments can be controlled by a formal procedure of expert elicitation and documentation(Garthwaite et al.,2005).Nevertheless,ex-pert judgments should be used with caution,not to replace ‘‘hard’’science(Apostolakis,1990).When assessing concep-tual model uncertainty,it is essential to adjust the prior probability to obtain the posterior model probability by con-ditioning of on-site measurements and observations.Expert elicitation of recharge model probabilities for the Death Valley regionalflow system103Different from general uses of expert elicitation for model parameterization and development,this paper uses the ex-pert elicitation to estimate prior model probabilities of alter-native models.With few examples of such an application of expert elicitation in model uncertainty assessment (Zio and Apostolakis,1996;Draper et al.,1999;Curtis and Wood,2004),this study is expected to provide theoretical and prac-tical guidelines for future applications of expert elicitation.This paper is focused on development of prior model proba-bilities using expert elicitation;discussion of using on-site data to further evaluate the alternative models is beyond our scope here.The expert elicitation is used in this paper to estimate prior probabilities of five recharge models developed for the Death Valley regional flow system (DVRFS),covering southwestern Nevada and the Death Valley area of eastern California,USA (Fig.1a).Due to existing and potential radionuclide contamination at the US Department of En-ergy’s Nevada Test Site (NTS)and the proposed Yucca Moun-tain high-level nuclear waste repository in the DVRFS,it is critical to predict contaminant transport in the region.Hydrologic and geologic conditions in the DVRFS are compli-cated,rendering multiple conceptualizations of the system based on limited data and information.Because conceptual model uncertainty can be significant,ignoring it (focusingonly on parametric uncertainty)may result in biased predic-tions and underestimation of uncertainty.While expert elic-itation was used for evaluating uncertainty of recharge and geological models (Pohlmann et al.,2007),this paper fo-cuses on the recharge models applied throughout the DVRFS.In the past few decades,several recharge models have been independently developed for Nevada by different researchers based on different scientific theories.These in-clude the Maxey–Eakin model (Maxey and Eakin,1949),the discrete-state compartment model (Kirk and Campana,1990;Carroll et al.,2007),the elevation-dependent chlo-ride mass balance model (Russell and Minor,2002;Russell,2004;Minor et al.,2007)and the distributed parameter wa-tershed model (Hevesi et al.,2003).It is unclear to scien-tists working in the DVRFS which recharge model should be used for groundwater flow and contaminant transport modeling.As recharge is the major driving force of ground-water flow,and thus contaminant transport,in the arid environment of the DVRFS,it is important to understand re-charge model uncertainty.Our ultimate goal is to incorpo-rate the recharge model uncertainty in our uncertainty analysis of DVRFS groundwater models.It is worth pointing out that recharge model uncertainty is prevalent and not limited to the DVRFS.Recharge is a fun-damental component of groundwater systems,andwithFigure 1(a)Boundaries of the Death Valley regional flow system,the Nevada Test Site,the proposed Yucca Mountain nuclear waste repository,and recharge rate estimates (m/d)of models (b)MME (modified Maxey–Eakin model),(c)NIM1(net infiltration model with runon–runoff component),(d)NIM2(net infiltration model without runon–runoff component),(e)CMB1(chloride mass balance model with alluvial mask),and (f)CMB2(chloride mass balance model with alluvial and elevation masks).104M.Ye et al.multiple recharge estimation methods(or models)avail-able,it is nontrivial to select the recharge estimation meth-od appropriate for a given environment(see review articles of Scanlon et al.,2002;Scanlon,2004).Scanlon et al.(2002) suggested using multiple methods to enhance reliability of recharge estimates.This is in line with the new concept of model averaging discussed above.The second section of this paper introduces the recharge models considered in the expert elicitation.Recharge esti-mates of the models are briefly compared in terms of their values,spatial distributions and statistical characteristics. In particular,we explain the reasons for treating recharge uncertainty as conceptual model uncertainty,rather than as parametric uncertainty.The process of expert elicitation is listed in the third section,followed by discussion of elic-itation results in the fourth section.Our conclusions are summarized in thefifth section.Description of thefive alternative recharge modelsThefive recharge models considered for the DVRFS are de-scribed briefly below;details of the models can be found in their original publications.Additional comparison of the models can be found in Rehfeldt(2004)and Pohlmann et al.(2007).Description of the geologic,hydrologic and hydrogeologic conditions of the DVRFS is beyond the scope of this paper,and the reader is referred to D’Agnese et al.(1997)and Belcher(2004)for further information on these topics.Modification of the Maxey–Eakin method(MME) Maxey and Eakin(1949)presented an empirical method (known as the Maxey–Eakin method)for estimating ground-water recharge as a function of precipitation.Since its inception,the Maxey–Eakin method has become the pre-dominant technique used for estimating annual groundwa-ter recharge in Nevada.The method estimates recharge viaR¼X Ni¼1C i P ið4Þwhere R is the estimated recharge,C i are the percentage adjustment coefficients,P i are the annual precipitation val-ues within zones of precipitation and N is the number of pre-cipitation zones.Maxey and Eakin(1949)utilized the precipitation map for Nevada developed by Hardman (1936)that includes hand-drawn contours based on weather station records and topography.The precipitation is distrib-uted amongfive isohyets(N=5)of5,8,12,15and20in. Assuming a steady-state basinflow condition in which dis-charge from a basin is approximately the same as recharge into the basin,the coefficients,C i,were developed through a trial-and-error method to attain a general agreement be-tween the volumes of estimated recharge and measured dis-charge for13basins in eastern and central Nevada.The coefficients,listed in Table1,increase in magnitude as the amount of precipitation increases while evapotranspira-tion and surface water runoff presumably decline.Note that the precipitation zone receiving less than8in./yr rainfall does not contribute to groundwater recharge.Given the incomplete coverage of the DVRFS domain by the Hardman precipitation map,Epstein(2004)modified the Maxey–Eakin model,hereinafter referred to as the modified Maxey–Eakin model(MME).The method uses the PRISM map(Precipitation Estimation on Independent Slopes Model)(Daly et al.,1994)so that the recharge is estimated in a consistent way over both the Nevada and California por-tions of the DVRFS.Considering uncertainty in the PRISM estimates of precipitation,the MME evaluates uncertainty of the recharge coefficients,C i,using an automated calibra-tion method based on91basins.Table1lists the mean coef-ficients of four precipitation zones(thus N=4in MME)used to estimate recharge of the DVRFS.Different from the Max-ey–Eakin method,the coefficient for the lowermost precip-itation zone is allowed to be nonzero.Although the MME model is more complicated than the original ME model,it is still the simplest model in the model set.The recharge map of the DVRFS estimated using the MME(with the mean coefficients)is shown in Fig.1b.Two net infiltration models(NIM)Hevesi et al.(2003)developed a distributed-parameter wa-tershed model,INFILv3,for estimating temporal and spatial distribution of net infiltration and potential recharge in the Death Valley region,including the DVRFS.The estimates of net infiltration quantify downward drainage of water across the lower boundary of the root zone,and are used as an indication of potential recharge under current climate con-ditions.Based on the daily average water balance at the root zone,the model comprehensively represents processes controlling net infiltration and potential recharge.The daily water balance includes the major components of the water balance for arid to semiarid environments,including precip-itation;infiltration of rain;snowmelt and surface water into soil or bedrock;runoff(excess rainfall and snowmelt);Table1Recharge coefficients for the Maxey–Eakin method and the modified Maxey–Eakin method(Epstein,2004)Maxey–Eakin method Modified Maxey–Eakin methodPrecipitation zone(in./yr)Coefficient Precipitation zone(in./yr)Coefficient 0to less than80.000to less than100.0198to less than120.0310to less than200.04912to less than150.0720to less than300.19515to less than200.15Greater than300.629 Greater than200.25Expert elicitation of recharge model probabilities for the Death Valley regionalflow system105surface water runon(overlandflow and streamflow);bare-soil evaporation;transpiration from the root zone;redistri-bution or changes in water content in the root zone;and net infiltration across the lower boundary of the root zone.Var-ious techniques were developed to estimate these quanti-ties and their spatial and temporal variability,which renders this method comprehensive but complicated.The model parameters(e.g.,bedrock and soil saturated hydrau-lic conductivity and root density)were adjusted through model calibration by comparing simulated and observed streamflow as well as basin-wide average net infiltration and previous estimates of basin-wide recharge.Two alternative net infiltration models with and without runon–runoff component(Hevesi et al.,2003)are consid-ered in this paper to represent the two opposite conceptu-alizations.Fig.1c and d depicts the averaged annual net infiltration estimates of the two models.Groundwater re-charge can be estimated from the net infiltration estimates by multiplying the net infiltration with coefficients related to rock hydraulic conductivity at the water table,since the net infiltration distribution only accounted for surficial characteristics of the system.For more details about the determination of the coefficients,the reader is referred to Belcher(2004).For convenience in this discussion,the two net infiltration models are also referred to as recharge models.Two elevation-dependent chloride mass balance models(CMB)The chloride mass balance(CMB)method estimates re-charge in basins(or any hydrologic systems)based on a bal-ance between chloride mass within hydrologic input and output components.The method assumes that chloride in groundwater within the basins originates from chloride in precipitation in mountain uplands and dry-fallout and is transported to adjacent valleys by steady-state groundwa-terflow(Dettinger,1989).At its most fundamental level, the method requires only estimates of annual precipitation in the recharge areas,total chloride input(chloride concen-trations in precipitation and recharge water)and total chlo-ride output(chloride concentrations in adjacent basin groundwater).The rate of recharge,R,can be calculated as(Maurer et al.,1996)R¼C p PC rÀC sw S wC rð5Þwhere C p is the combined wet-fall and dry-fall atmospheric chloride concentration normalized to precipitation,P is the mean annual precipitation rate,C r is the chloride concentra-tion in recharge water and C sw is the chloride concentration in surface water runoff S w.For individual basins,recharge rate can be estimated from this information if the following assumptions are met(Dettinger,1989):(1)there are no other major sources or sinks for chloride in the system;(2)surface runoff is small in comparison to groundwaterflow;and(3) the recharge areas are correctly delineated.Russell and Min-or(2002)extended the chloride mass balance approach to account for the elevation of precipitation,the limited quan-tities of recharge that are thought to occur on low-elevation alluvial surfaces,and uncertainty inherent in the data.This elevation-dependent chloride mass balance approach was applied by Russell and Minor(2002)to a7900-km2region of the Nevada Test Site(NTS)and vicinity within the DVRFS.Although this recharge/elevation relationship simulates recharge at all elevations,several studies suggest that sig-nificant groundwater recharge does not occur through low-elevation alluvial sediments in southern Nevada.Russell and Minor(2002)thus developed two models to address this uncertain conceptualization of low-elevation recharge.The first model assumes that all land surface areas covered by alluvial sediments receive negligible recharge based on the results of previous studies and soil-water chloride pro-files of40boreholes completed in unsaturated alluvium within the NTS(Russell and Minor,2002).This model is called the CMB model with alluvial mask.The second model assumes that the elevation of the lowest perennial spring that discharges from a perched groundwater system in the study area represents the lowest elevation at which signifi-cant recharge occurs.This spring is Cane Spring,which is lo-cated at an elevation of1237m above mean seas level. Coincidentally,this is approximately the same elevation (1200m)that Harrill(1976)and Dettinger(1989)consider to be the minimum at which precipitation makes a signifi-cant contribution to recharge in desert basins of central and southern ing the concept of a recharge cut-off elevation,Russell and Minor(2002)define a zone of zero recharge that encompasses all elevations below1237m plus elevations above1237m that are covered by alluvium.This model is called CMB with both elevation and alluvial masks. To assess uncertainty in the model parameters and mea-surements(e.g.,precipitation and chloride concentration in spring water),Russell and Minor(2002)developed a Monte Carlo method to estimate multiple realizations of the recharge estimates.The two models were further ex-tended in Russell(2004)and this study to include more ba-sins in Nevada and cover the DVRFS.Fig.1e and f depicts mean recharge estimates of the two CMB models. Summary and discussionThefive recharge models are summarized as follows: MME(Fig.1b):modified Maxey–Eakin model using the mean coefficients.NIM1(Fig.1c):net infiltration model with runon–runoff component.NIM2(Fig.1d):net infiltration model without runon–run-off component.CMB1(Fig.1e):chloride mass balance model with alluvial mask(mean estimates only).CMB2(Fig.1f):chloride mass balance model with alluvial and elevation masks(mean estimates only).Fig.1illustrates similarities and differences of the re-charge rate estimates(m/d)of thefive models,and Table 2lists the total recharge estimates(m3/d)for the entire DVRFS by each method.The MME gives the highest recharge estimate,and the CMB models give higher estimates than the NIM models.Due to the runon–runoff component con-sidered in NIM1,the recharge estimate of NIM1is higher than that of NIM2,while spatial patterns of the recharge estimate are similar in the two models.Because of the extra106M.Ye et al.elevation mask considered in CMB2,the recharge estimate of CMB2is lower than that of CMB1;for the same reason, spatial patterns of the recharge estimate are different in the two models(less recharge is estimated in southern Ne-vada in CMB2).The recharge estimate of the MME has the smoothest spatial distribution,due to the four precipitation zones.The different recharge estimates are viewed as a re-sult of conceptual model uncertainty,rather than paramet-ric uncertainty,since they are caused by simplification and inadequacy/ambiguity in describing the recharge process and not by uncertainty in recharge measurements them-selves(Wagener and Gupta,2005).Given thefive recharge models,which model(or models) should be used for groundwater modeling?Is it reasonable and justifiable to select a single model and to discard others based on expert judgment?How should uncertainty of the recharge models be assessed?The expert elicitation is used to answer these questions,and ultimate results of this ex-pert elicitation are the prior model probabilities essential to the BMA for assessing the conceptual model uncertainty. Process of the expert elicitationWhile several processes of expert elicitation have been sug-gested in the literature(e.g.,Hora and Iman,1989;Bonano et al.,1990),the process proposed by Keeney and von Win-terfeldt(1991)was followed,since it is closely pertinent to eliciting probability from experts and has been applied to model probability elicitation(Zio and Apostolakis,1996). The formal process consists of the seven steps listed below. Implementation of the process for the recharge models is also described.Step1:Identification and selection of elicitation issues The elicitation issues are the questions posed to the ex-perts that require their answers.The following three issues are considered for assessing the recharge model uncertainty: (1)Is the model set complete,given the objective of theanalysis?BMA requires that alternative models arecomprehensively exhaustive(all alternative modelsare included in the model set).Since this requirementcannot be satisfied in an absolute sense,we elicitfrom the experts whether there are other alternativemodels that are comparable in importance to thefivemodels and should be considered.(2)What are the plausibility ranks of these models,giventhe objective of the analysis?Whereas ranking ofmodel plausibility is qualitative and the ranks cannotdirectly give the prior model probability,the modelranking helps experts evaluate relative plausibilityof the models before they estimate prior modelprobability.(3)What is the probability value that best represents theconfidence you would place on each recharge model,given the objective of the analysis?Model probabili-ties are the ultimate goal of the expert elicitation,and will be used directly in the BMA to calculate theposterior model probability through Eq.(2).Step2:Identification and selection of expertsExpert elicitation requires three types of experts:gener-alists,specialists and normative experts.In this study,the generalists should be knowledgeable about various aspects of the recharge models and the broader study goals(in this case,assessing groundwaterflow and contaminant transport in the DVRFS).They typically have substantive knowledge in one discipline(e.g.,geology or hydrology)and a general understanding of the technical aspects of the problem. While the generalists are not necessarily at the forefront of any specialty within their main discipline,the specialists should be at the forefront of one specialty relevant to the recharge models.The specialists often do not have the gen-eralists’knowledge about how their expertise contributes to the broader study with respect to recharge model uncer-tainty analysis.Normative experts typically have training in probability theory,psychology and decision analysis.They assist generalists and specialists in articulating their profes-sional judgments and thoughts so that they can be used in a meaningful way in the conceptual model uncertainty assess-ment.A high-quality elicitation requires the teamwork of all three types of experts.Selecting experts is a time-consuming process,and may take more than a year for a full-scale elicitation(e.g.,having international nomination of experts and forming an expert panel of international scientists,as in Hora and Jensen, 2005).With practical limitations,we selected national and state experts,who were believed well-qualified owing to their familiarity with the hydrogeologic conditions of the DVRFS and their research at the forefront of recharge esti-mation in semi-arid environments of the southwestern US. Five specialists,two generalists and one normative expert were identified.The normative expert had an advisory role and was not involved in evaluating the recharge model uncertainty.Step3:Discussion and refinement of elicited issuesThis step allows discussion and refinement,if necessary, of the issues and quantities that will be elicited.While Kee-ney and von Winterfeldt(1991)suggest completing this step by a1-day meeting of all experts,such a meeting was con-sidered unnecessary for this project.Instead,one month before the elicitation,the experts received the three clearly stated elicitation issues,as well as original publica-tions of thefive recharge models and references about con-ceptual model uncertainty,BMA,prior model probability and expert judgment.The experts studied these materials, and some discussed details of the models with us and requested more reading materials.Step4:Training for the elicitationLed by the normative expert,the training was conduct in two meetings in thefirst half day of elicitation.In thefirst training meeting,the normative expert introduced theTable2Recharge estimates(m3/d)of thefive rechargemodels in the DVRFSRecharge model DVRFS(m3/d)MME596,190.8NIM1341,930.6NIM2282,223.1CMB1385,213.7CMB2365,647.2Expert elicitation of recharge model probabilities for the Death Valley regionalflow system107。

Bayesian Inference1. Bayes’ Theorem and Heuristic Applications Ex. 1. (Learning) W is either tossing a coin or a die. H tries to guess which, but W only reports 1 for Heads, 2 for Tails, if she is tossing a coin; or 1,..., 6 if she is tossing a die. Initially H gives equal chances for both behaviors: P(coin) = P(die) = 1/2. Suppose W reports: 1, 2, 5, 1,... (1) What does H think when the first 1 is reported; (2) after 1,2; (3) after 1,2,5; (4) after 1,2,5,1?(1) P(first is 1) = 1/2×1/2 + 1/2×1/6 = 1/3;P(first is 1 and W uses coin) = 1/2×1/2 = 1/4; soP(coin*first is 1) = 1/4÷1/3 = 3/4; up from 1/2!(2) P(coin*1, 2) = 3/4×1/2÷(3/4×1/2+1/4×1/6) = 9/10.(3) But after 5 he knows it must be a die after all! Math:P(coin*1, 2, 5) = 9/10×0÷(9/10×0+1/10×1/6) = 0.(4) Later results won’t matter, as shown by:P(coin*1, 2, 5, 1) = 0×1/2÷(0×1/2+1×1/6) = 0. This uses Bayes’ Theorem. Suppose the state of thecworld is either A or A. We observe data D. Recall,P(D*A) = P(A1D)/P(A). Therefore (Venn diagram!),c cP(A*D) = P(D*A)P(A)/(P(D*A)P(A) + P(D*A)P(A)). This is valid whenever D and A are events whose probabilities are defined.12Ex. 2. (Diagnostic testing) A = person is ill,D = medical test says she is ill (result “positive”),P(D *A) / se = sensitivity , P(D *A )/ sp = specificity ;c c P(A) / p = prevalence of illness.How likely it is that a person is ill if her test is positive?P(A *D) = se×p/(se×p + (1 - sp)×(1 - p)),so if sp < 1, then for rare illnesses (or small p) the probability can be very small. E.g., if se = 1, sp = .9,and p = .01, then P(A *D) = .01/(.01 + .1×.99) ..1 only!Ex. 3. (Learning about prevalence) Study a population with unknown prevalence (of an illness), so initially p -i U[0,1]. Examine individuals one by one, X = 1, if ith is i 1n ill, otherwise X = 0, i = 1,..., n. Define S = X +...+ X -Bin(n, p). What should we think of p now? The prevalence takes some value in [p, p+dp] with probability = length of the interval = dp, so the probability of seeing S cases and p in this range is Divide this with its integral and cancel the factorials to see that the density is proportional to p (1 - p). This S n - S must be the density of the beta distribution Be(", $),'(" + $)/ '(") '($)×p (1 - p), with " = S + 1, $ = n "-1$-1- S + 1. (Graph for n = 10, 100, and S = n/10, n/2!)3Here we used the Bayes’ theorem in the form,where g(p) = 1, for p ~ U[0,1], is the so-called prior distribution of the prevalence.Ex. 4. (Backcasting) Consider a formal language that has words of random length that consist of vowels (v)and consonants (c): vcc, cvcvcv,... Suppose (1) that the probability of a vowel is 0 < B < 1 in every position,and (2) that a vowel is followed by a consonant with c probability 0 < p < 1, and a consonant if followed by a v vowel with probability 0 < p < 1. Thus the words are realisations of a Markov chain with state space {v, c}and transition matrixAssumption (1) implies that B = (B , 1 - B ) is theT invariant distribution of the chain satisfying B = B M .T T v v c Therefore,B = p /(p + p ). Suppose the letter number t> 1 is a vowel. What is the probability that the previous letter was also a vowel? By Bayes’ theorem,c c v P(t-1 is v *t is v) = (1 - p )B /((1 - p )B + p (1 - B)).This the posterior probability of a vowel at t-1, given that at t there is a vowel; B is the prior probability.4。

姜雪萍,陈艳君,朱富成,王芳,孙传伯,赵群*(皖西学院生物与制药工程学院,安徽六安237012)摘要：通过D N A 条形码I T S 序列,探讨鉴别中药百合及其混伪品的有效方法。

从G en B an k 下载了8种百合属(Lilium )、2种假百合属(Notholirion )和1种老鸦瓣属(Amana )物种的I T S 部分核D N A 序列。

经过比对后,I T S 序列长度为610b p ,G +C 平均含量为60.1%。

用贝叶斯法(B ayesian in feren ce ,B I )和简约法(M axim u m p arsim on y ,M P )分别构建系统发育树。

结果表明,百合药材的3种正品基源植物的单倍型在B I 和M P 2种系统发育树中都各自聚类为单系,并不和其他混伪品的物种混合。

因此,I T S 序列可以作为D N A 条形码对中药百合进行分子鉴定,能够将中药百合的正品基源植物和其他混伪品的植物种类区分开来,从而鉴定中药百合的真伪。

关键词：百合;I T S ;D N A 条形码;分子鉴定中图分类号：R 282.5文献标识码：A 文章编号：2095-0896(2019)09-001-03JIANG Xue-ping et al.(West Anhui University ,Lu ’an ,Anhui 237012)Abstract An effective method to identify the traditional Chinese medicine (TCM )Lily and its adulterants was studied by DNA barcoding ITS sequence in this study.Sequences of nuclear DNA from ITS of 8species from Lilium ,2species from Notholirion ,1species from Amana were downloaded from the GenBank.These sequences were compared and edited ,the length of the sequence fragment was 610bp ,and the percentage of G +C was 60.1%.Two methods were used to condstruct the phylogenetic trees ,including Bayesian inference (BI )and maximum parsimony (MP ).The results showed that the haplotypes of the three genuine Lily species formed their own monophyletic group respectively in both the BI and MP phylogenetic trees ,and distinguished from its adulterants ,so ITS sequence could be used as the DNA barcoding for the identification of Lily and its adulterants.Key words Lily ;ITS ;DNA Barcoding ;Molecular identification NA Barcode Identification of Lily Medicinal Materials百合药材的DNA 条形码鉴定中药百合为百合科(L iliaceae )百合属(Lilium )植物卷丹(Lilium lancifolium )、百合(L.brownii )或细叶百合(L.pumilum )的干燥肉质鳞叶,为我国传统的名贵中药,具有药、食两用功能,有“蔬菜人参”的称号[1-4]。