Bayesian changepoint analysis for atomic force microscopy and soft material indentation

Ba yes ian changepoint analy sis for atomic f orce micr oscop y and soft material indentation Daniel Rudo y , Shel ten G . Y u en, Rober t D . Ho we and P atric k J. W olf e † Har v ard University , Cambridge, USA. [Manuscript submitted for publicatio n September 2008. Revised Septembe r 2009.] Summary . Materi al inden tation studies, i n whi ch a pr obe is broug ht into controll ed p hys ical contact wit h an e xperimen tal samp le, hav e l ong been a p rimar y mea ns by which scientists character ize the mecha nical prope r ties of mater ials. More recentl y , the advent of at omic force microscopy , which ope rates o n the same funda mental p rincip le, has in tur n rev olution ized the nan oscale analysis o f soft biomate ria ls such a s cel ls a nd tissues. This p aper ad dresses the inferentia l problems associated with mater ial in dentati on and atomic force microscopy , thro ugh a frame- work for the chan gepoin t an alysis o f pre- a nd post-conta ct data that is ap plicable to experi- ments across a variet y of physical scales. A hie rarchical Bayesian mod el i s p roposed t o ac- count for e xperimen tally obser ved changepoin t smoothness constraints and measurement er- ror v ariab ility , with efficient Monte Carlo methods dev eloped and employed to realize inference via po sterio r sa mpling for pa rameters such as Y oung’ s modul us, a key qua ntifier of mater ial stiffness. These re sults are the first t o provide the materia ls science community w ith r igor- ous inference pro cedures an d uncer tai nty qua ntification , via optimized a nd fully automated high-t hrough put algor ithms, implemented as the publicly av ailable software package Ba yesCP . T o demonstrate the consistent accuracy and wi de appli cability of this appro ach, results are shown f or a variety of data s ets from both macro- and micro-materi als experiments—in cluding silicone, neurons, and red blood cells—conducted by the authors and others. K eywords : Change point de tection; Constraine d switching regressions; Hierarchical Bay esian models; Indenta tion testing ; Markov c hain Mon te Carlo; Materia ls sc ience; Y oung’ s modulu s 1. Introduction This a r ticle develops a hier archical Ba yesian approach for contact p oint determination in material indentation studies and atomic for ce microscopy (AFM). Contemporary applica- tions in mater ials science and bio mechanics range fro m ana lyzing the r esp onse of nov el nanomateria ls to deformation (W ong et al., 1997), to characterizing diseas e through me- chanical pr op erties of cells, tissues, and or gans (Costa, 2004). Exp erimental procedures and analyses, howev er, remain broadly similar acr oss these different material types and physical scales (Gouldstone et a l., 20 0 7), with the scientific aim in all c ases b eing to characterize how a given mater ial sample deforms in resp onse to the applica tion of an external force. As illustrated in Figure 1 overleaf, inden tation exp eriments emplo y a prob e (or can- tilever arm, in the case of AFM) to apply a con trolled force to the ma terial sample. This inden ting prob e displac e s the sample while concurrently mea suring resistive force, with the resultant for ce-p osition da ta used to infer material prop erties suc h as stiffness (in ana logy † A ddr ess for c orr e sp ondenc e: P atric k J. W olfe, Statistics and Information Sciences Lab oratory , Harv ard Universit y , Ox ford St reet, Cam bridge, MA 02138-2901, U SA E-mail: wolfe@stat.har vard.edu; s oftware URL: http://sisl.se as.harvard.edu /Ba yesCP .html 2 D . Rudoy , S. G. Y uen, R. D . Howe and P . J. W olfe 0.1 N Force Sensor Indenter Material Sample 1.3 N x x Fig. 1. Diagram of a macro-scale i ndenta tion experiment in which a spher ical probe, a ttached t o a force sensor , ind ents a materi al s ample and deforms it by a distance δ . I n this h ypothe tical ex ample, a net change of 1 . 2 Newto ns in resistiv e force is consequen tly obser ved. to compressing a s pring in o rder to exp erimentally deter mine its spring constant ). Prior to subsequent data analysis , a key technical pr oblem is to determine pr ecisely the moment at which the pro be comes into contact with the material. Sample pr eparation techniques and siz e s frequently preclude the direc t measur e men t of this contact p oint, and he nc e its inference from indent er for ce-p osition da ta forms the sub ject of this article. A t present , practitioners acr oss fields lack an a greed-up on standard for contact p oint determination; a v ariety of ad ho c data pre-pro cessing metho ds are used, including ev en simple visual inspection (Lin et a l., 2007a). Nevertheless, it is well re cognized that precise contact p oint deter mina tion is necessary to a ccurately infer material prop erties in AFM inden tation exp eriments (Cr ick a nd Yin, 2007 ). F or ex a mple, Dimitriadis et al. (2 002) show that for s ma ll displacements of thin films, e s timation err ors on the o rder of 5 nm for a 2 . 7 µ m sample ca n cause an increase of nearly 200% in the estimated Y oung’s mo dulus—the principal quantifier of material stiffness. When s oft materia ls such as cells are studied a t microscopic scales, for exa mple to determine biomechanical disease mar kers (Costa, 2004 ), the need for robust and rep eata ble AFM analyses beco mes even grea ter (Lin et al., 20 07a). In this article, we present the first fo r mulation of the contact p oint determination tas k as a statistic al changep oint problem, and subsequently employ a switching regr essions mo del to infer Y o ung’s modulus. Section 2 s umma r izes the basic principles o f materia l indentation, showing that the resultant for ce-displacement curves are often well desc r ib e d b y low-order po lynomials. Section 3 intro duces a corr esp onding family of Bay esian mo dels designed to address a wide ra nge of exp erimental conditions, with sp ecialized Markov c hain Monte Carlo samplers for inference developed in Section 4. F o llowing v alidation of the prop osed inference pro cedures in Section 5, they are employ ed in Sectio n 6 to infer material prop er ties of mouse neurons and human r ed blo o d cells from AFM force-p os ition data. The article concludes in Section 7 with a discussio n o f pro mising metho dological and practical extensio ns . 2. Material indentation 2.1. Indentation e xperiments and data Indent ation experiments pro ceed by car efully moving a pro be from an initial noncontact po sition into a material sample, a s shown in Figure 1, while measur ing the resistive for ce at so me pre s crib ed temp oral sampling rate. After a small deforma tion is made, the pro be retracts to its initial p os ition; dur ing this stag e the resis tiv e for c e decreas es with every subsequent mea surement. A t the co nclusion o f ea ch such exp er imen t, tw o force-p osition curves are pro duced, exa mples of which a re sho wn in Figure 2. In this a rticle we consider Bay esian changepoin t analysis f or atomic force microscopy and soft mater ial indenta tion 3 −5 0 5 4 4.5 5 5.5 6 6.5 7 7.5 8 Force (N) −4 −2 0 4.05 4.1 4.15 4.2 4.25 4.3 4.35 4.4 Distance from initial probe position (mm) (a) S ilicone indentation data −3 −2 −1 0 0 200 400 600 800 1000 1200 1400 1600 1800 2000 Force (pN) −3 −2 −1 −20 −10 0 10 20 30 40 50 60 70 80 Distance from initial probe position ( µ m ) (b) Red blo od cell in dentatio n data Fig. 2. Examples of force-position curves from which materi al stiffness proper ties are to be inferred, with force mea surements dur ing indent ation and sub sequent p robe retraction shown in bla ck a nd grey , respectively , and subsets of the in denta tion data shown near the presumed contact points. Experi ments were perfor med using a mechanical ar m f or the soft silicone sa mple, and an a tomic force microscope f or the red b lood cell sample; n ote the di fferences i n physical scale and no ise lev e l. only the forward-indentation cur ves, a s is standard pr actice (Lin et al., 2007 a), though the metho ds we pr esent are extensible to r etraction da ta whenever suitable mo dels are av a ilable. Despite significa n t difference s in physical sca le and no ise level, the curves o f Fig ure 2 (and indeed, most indentation data se ts) feature sev eral commo n characteris tics. In the pre-contact reg ion, force resp ons e appea rs linear in the p osition o f the indent er; drift due to exp erimental co nditions is o ften present, yielding a no n-zero slo pe as in Fig ure 2(a). The po st-contact data app ear well mo deled by low-order p olynomial functions of the co rresp ond- ing displacement, and indeed, the key a s sumption o f such e x per iment s is that, conditioned upo n knowledge of the geo metry of the indenting pr ob e , the relationship be tw een the degr ee of material deforma tion and the measur ed resistive for c e dep ends in a known way on the material stiffness. This is analo gous to the cas e in which an ideal spring is c o mpressed a sp ecified distance δ b y a force of known magnitude; a measure of the spring’s stiffness is given by the spring constant k , whic h may be calculated using Ho o ke’s law: F = − k δ . 2.2. Contact mechanics and the Her tzian model Indent ation data such as thos e shown in Figur e 2 are t ypically acquir e d for line ar-elastic materials. This implies not only that the material instantaneously returns to its orig inal shap e following the cessatio n of an ex ternal forc e , but also that the relatio ns hip b et ween the applied str e ss (force p e r unit a rea) and the resultant strain (deforma tion p er unit length) is linear. This r atio is known as Y oung’s mo dulus —the primary quantifier of mater ial s tiffness int ro duced ab ov e—and is rep or ted in units o f Pascals. In small-defor mation ex p er imen ts to infer Y o ung’s mo dulus E for linear-ela stic materia ls, bo th the indenter geometry and the measured resis tive force come into play , by w ay of the so-called Hertzian mo del (Lin et al., 200 7a). Spec ific a lly , the relationship b etw een the sample deformation depth δ and the mea sured resistive force F takes the form F ∝ E · δ β , (1) 4 D . Rudoy , S. G. Y uen, R. D . Howe and P . J. W olfe where the constant of pro p or tionality and the fixed parameter β depend on the geometry of the indenter tip in a known wa y . Examples include indentation by a shar p pyramid with tip angle 2 φ , or a sphere of r adius R , whereup on (1) takes the follo wing forms: F = 1 . 5 tan φ 2(1 − ν 2 ) E · δ 2 (p yramid) or F = 4 R 1 / 2 3(1 − ν 2 ) E · δ 3 / 2 (sphere), (2) with ν a known, dimensionless quantit y termed Poisson’s ratio. A subsequent fitting of the Hertzian mo del of (1) to exp erimental data th us allows o ne to obtain an estimate for Y oung’s mo dulus E once the po st-contact region has b een iden tified. (Below we retain the standard practitioner notation ( E , F , R, δ, ν, φ ), a s distinct from other v a riables to follow.) 2.3. Pre- and post-contact data regimes Observe that the Hertzian mo del of (1) p osits a relationship be tw een force a nd indentation depth , wherea s the data of Figur e 2 a re seen to b e a function of the p osition o f the indenter. The Hertzian mo del th us describ es the underlying physics of the p ost -contact s tage o f a t ypical indentation exp eriment, while the measur e d data als o compr ise a pr e -contact stage. As w e detail b elow, the union of these tw o regimes is well des crib ed b y a switching r egressio ns scenario, with force meas urements prior to co n tact typically linear in the p osition o f the prob e. The corr esp o nding intercept represe n ts the equilibr ium tare r equired to achieve zer o net force , while the slop e captures constant-velocity drift in for ce measurements that can arise in a v ariety o f exp erimental s ettings. In b oth the pr e- and p os t-contact data regimes, it is standar d to as s ume indep endence o f measurement err ors—a re asonable assumption for a ll practically achiev a ble force sampling rates. Er rors pr ior to co nt act a rise due to force senso r vibration in the exp er imental medium, thermal v aria tions, and other effects—wherea s after contact they also dep end on interactions betw een the prob e and the sample such a s frictio na l forces. Consequently , while erro r v ariances in these reg imes can b e exp ected to differ in pra ctice (analysis of AFM data fro m a r ed blo o d cell, shown in Figure 6 a nd descr ibed la ter, re veals one such example), their relative mag nitudes are not known a prior i. As a fina l consider a tion, uncertain ty in the rep orted position of the indenting probe is typically several or ders of magnitude smaller than the distance b etw een co nsecutive force sa mpling p oints, and can sa fely be disregar ded. As noted by Crick and Yin (2007 ), all of the sources o f v ariability men tioned ab ov e can lead to larg e relative errors in resis tive force measurements near the contact point; Figure 2(b) illustr ates a typical scenario. This makes manual identification of the co n tact po in t difficult in many cases, and motiv ates the mo del-based a pproach that we now descr ibe. 3. Ba yesian changepoint model Having outlined the bas ic principles of material indentation, we now pr o ceed to fo rmulate a Bay es ian mo del for for ce-p osition data that encompasses b oth pre- and p o st-contact regimes. In light of the discuss ion ab ov e, the corr e spo nding task o f contact po int determination is recogniza ble as a changep oint estimatio n problem in the co ntext of s w itching regr essions. In tur n, the Hertzian mo del of (1) implies that an estimate o f Y o ung’s modulus can b e obtained as a line ar funct ion of the lea ding p ost-contact regressio n co efficient . Con trol ov er ex per imen tal conditions implies that each inden tation data set contains precisely one changepoint, thus o b viating a n y need to estimate the pr e s ence or num b er of contact p oints. Bay esian changepoin t an alysis f or atomic force microscopy and soft mater ial indenta tion 5 During an indentation exp eriment, the indenter moves contin uously throug h a seq ue nc e of n positions x = ( x 1 , x 2 , . . . , x n ) ′ and records a for ce measur ement y i at eac h p osition x i , resulting in a se q uence of for ce measurements y = ( y 1 , y 2 , . . . , y n ) ′ . As w e shall la ter treat mo de ls for soft materia l indentation, in whic h the pre- and p ost-co nt act curves are constrained to b e contin uous a t the regres sion ch angep oint, we beg in by introducing a contin uous parameter γ ∈ (1 , n ) denoting the c ontact p oint index , with the co rresp onding c ontact p oint at which the indenter first co nt acts the sample denoted by x γ ∈ ( x 1 , x n ). 3.1. Data lik elihood for i ndentation experiments W e adopt a classica l switching regres sions scenario for our mo del, in which y is as sumed to b e a p olynomia l function of known degree d 1 in p osition x pr ior to contact, and of known degree d 2 in deformation depth δ = x − x γ after contact with the sample is made. This for m ulation enco mpasses the Hertzian mo del of (1) if fractional p ow ers a re allowed; how ever, for clarity o f pr e s ent ation w e consider d 2 to be an integer unless otherw is e noted. Letting p = d 1 + d 2 + 2 denote the num b er of regres sion co efficients in our mo del, and with n the num ber o f da ta p oints, the cor resp onding design matrix is hence of dimension n × p . W e subseq uen tly employ the subscript γ to deno te any qua n tity that depends on γ , and index via subscripts 1 and 2 t he pre- and post- c ontact re gression r egimes. W e denote the r egressio n co efficients by β 1 ∈ R d 1 +1 and β 2 ∈ R d 2 +1 , with design matrices X 1 ,γ , X 2 ,γ defined as follows, for ⌊ γ ⌋ the lar gest integer less than o r equal to γ : X 1 ,γ =      1 x 1 · · · x d 1 1 1 x 2 . . . x d 1 2 . . . . . . . . . . . . 1 x ⌊ γ ⌋ · · · x d 1 ⌊ γ ⌋      , X 2 ,γ =       1 x ⌊ γ ⌋ +1 − x γ · · ·  x ⌊ γ ⌋ +1 − x γ  d 2 1 x ⌊ γ ⌋ +2 − x γ . . .  x ⌊ γ ⌋ +2 − x γ  d 2 . . . . . . . . . . . . 1 x n − x γ . . . ( x n − x γ ) d 2       . (3) The o bserved data y may likewise be partitioned in to pre- and p ost-co n tact vectors y 1 ,γ = ( y 1 , y 2 , . . . , y ⌊ γ ⌋ ) ′ and y 2 ,γ = ( y ⌊ γ ⌋ +1 , y ⌊ γ ⌋ +2 , . . . , y n ) ′ , and, following our discussion in Section 2.3 regarding the noise characteristics typical of inden tation exp eriments, we assume indep endent and Nor mally distr ibuted additive error s, with unknown v a riances σ 2 1 , σ 2 2 . T hus for 1 ≤ i ≤ n , w e have that y i is distributed as follows: y i ∼ ( N ( X 1 ,γ β 1 , σ 2 1 ) if 1 ≤ i ≤ ⌊ γ ⌋ , N ( X 2 ,γ β 2 , σ 2 2 ) if ⌊ γ ⌋ + 1 ≤ i ≤ n . (4) Note that the statistical model of (4) is consistent with the Hertzian mechanics mo del of (1), resulting in a p ost-co n tact force -resp onse curve that is a power of the displac ement δ = x − x γ of the material sample, ra ther than the p osition x of the indenter. Howev er, when d 2 is a n integer, the co efficient s of these tw o p olyno mials ar e related by a simple linea r transformatio n. Consider , for instance , a q uadratic curve in δ given by f ( δ ) = a 0 + a 1 δ + a 2 δ 2 . One may rewrite f ( δ ) a s a q uadratic p olynomial in x as follows: a 0 + a 1 ( x − x γ ) + a 2 ( x − x γ ) 2 = b 0 + b 1 x + b 2 x 2 , (5) where b 2 = a 2 , b 1 = a 1 − a 2 x γ , and b 0 = a 0 − a 1 x γ + a 2 x 2 γ . This tra nsformation enables X 2 ,γ to be r eformulated directly in terms of indenter position x , such that it no longer 6 D . Rudoy , S. G. Y uen, R. D . Howe and P . J. W olfe depe nds c ontinu ously on x γ , in contrast to (3). T ransformatio ns a kin to (5) do not apply , how ever, when the d 2 is a fraction, as in the case o f (2) for a spherica l indenter, or when a contin uity constr aint is e nfo r ced at the c hangep oint; we detail such c a ses b elow. 3.2. A general parame tric Bay esian model for m aterial indentation The likeliho o d of (4), tog ether with the presence of genuin e prior information dictated by the under lying physics of materia l indentation exp eriments, sugg ests a na tural hier archical Bay es ian mo del. In contrast to the semi-conjugate appr oach taken by Ca rlin et al. (19 9 2), we detail b elow a fully conjugate mo del, as this allows for analytica l simplifications that we have observed to be impor tant in practice. Integrating out n uisance para meters improves not only the mixing of the chains under lying the resultant alg orithms and inferential pro cedure s , but also their computational tra ctability when data sizes gr ow large . W e sp ecify pr ior distributions for all mo del par ameters, including the contact p o int index γ ∈ (1 , n ), the pre- and p ost-co n tact reg ression co efficients β 1 and β 2 , and the e r ror v ariances σ 2 1 and σ 2 2 . F or i ∈ { 1 , 2 } , we then as s ume that β i ∼ N ( µ i , σ 2 i Λ − 1 i ), with Λ i a ( d i + 1) × ( d i + 1) dia gonal p ositive definite matrix and µ i ∈ R d i +1 . A standard inv er se- Gamma conjugate prio r I G ( a 0 , b 0 ) is a dopted for b oth v a riances σ 2 1 and σ 2 2 . Finally , recalling that a single re gression changep oint is always a ssumed within the data reco rd, we employ a uniform prior distr ibution on the in terv al (1 , n ) for the contact p oint index γ . In certain exp erimental settings , wher eup on the initial p osition of the indenter is k nown to b e at le a st a cer tain distance from the sample, a n infor mative prior distribution may well b e av aila ble. Because of sensitivity to hyp e rparameter s, we follow standar d practice and adopt h yp er- priors for incr eased mo del robustness. A Gamma prio r is ass umed on b 0 so that b 0 ∼ G ( κ, η ), how ever, we determined fro m simulations that the p os terior estimato rs consider ed were not sensitive to the prior para meter s µ i and Λ i , and s o did not employ a n additional level of hyperprior hier arch y for the re g ression co efficients. F or notatio na l conv enience, we let 1 k denote the 1-vector of dimension k whose entries are all equal to one and 0 the zero matrix of appropriate dimension, and define the v ariables y , β , µ , X γ , Λ , Σ γ , and Σ as follows: y =  y 1 ,γ y 2 ,γ  ∈ R n × 1 , β =  β 1 β 2  ∈ R p × 1 , µ =  µ 1 µ 2  ∈ R p × 1 ; X γ =  X 1 ,γ 0 0 X 2 ,γ  ∈ R n × p , Λ =  Λ 1 0 0 Λ 2  ∈ R p × p ; Σ γ = diag( σ 2 1 1 ⌊ γ ⌋ , σ 2 2 1 n −⌊ γ ⌋ ) ∈ R n × n , Σ = dia g ( σ 2 1 1 d 1 +1 , σ 2 2 1 d 2 +1 ) ∈ R p × p . The p osterior pro ba bilit y distribution o f the mo del para meters ( γ , β , σ 2 1 , σ 2 2 , b 0 ), conditioned on the observ a tions y a nd the fixed model par a meters ψ , ( µ , Λ , a 0 , κ, η ), is then: p ( γ , β , σ 2 1 , σ 2 2 , b 0 | y ; ψ ) ∝ p ( y | β , σ 2 1 , σ 2 2 , γ ) p ( β | σ 2 1 , σ 2 2 ; µ , Λ ) p ( σ 2 1 | b 0 ; a 0 ) p ( σ 2 2 | b 0 ; a 0 ) p ( b 0 ; κ, η ) p ( γ ) ∝ σ − 2( a 0 − 1) 1 e − b 0 /σ 2 1 · σ − 2( a 0 − 1) 2 e − b 0 /σ 2 2 · b κ − 1 0 e − b 0 /η · ( | Σ γ | | Σ | ) − 1 2 exp  − 1 2  ( y − X γ β ) ′ Σ − 1 γ ( y − X γ β ) + ( β − µ ) ′ Σ − 1 Λ ( β − µ )  . (6) T o confirm robustnes s , we a lso studied the effect of repla c ing the diagonal prior cov a r i- ance Λ for the pre- and p ost-contact r egressio n co efficients β 1 and β 2 with an a ppropriately adapted g -prior (Zellner , 1 986) such that β i | ρ i , γ ∼ N ( µ i , σ 2 i ρ 2 i ( X ′ i,γ X i,γ ) − 1 ), with ρ 2 i a scale para meter to which we ascrib ed a diffuse inv erse-Ga mma hyper prior. W e o bserved Bay esian changepoin t an alysis f or atomic force microscopy and soft mater ial indenta tion 7 no meas ur able effect of this change in priors o n the resulting inference—further confirming the inse ns itivit y of the a do pted mo del to the pr ior cov aria nc e of the reg ression co efficients. Moreov er, efficient sa mpling from the conditional distr ibutio n of ρ 2 i is pr ecluded by its de- pendenc e on the co nt act p oint index γ , r educing the overall efficacy of this approach in the Marko v chain Monte Car lo appr oaches to p osterio r sa mpling de s crib ed in Section 4 b elow. 3.3. Smoothness constraints at the changepoi nt It is often the case that force- po sition curves ar e cont inuous a t the contact point x γ . Es- pec ially for soft ma ter ials such as the red blo o d cells we consider in Section 6 , it is to be exp ected that the change in the forc e measurement is smo oth, and a con tinuit y con- straint can serve to regular ize the solution in cases where many differ e n t fits will have high likelihoo d. Imposing smo othness co nstraints dates back to at lea st Hudson (19 66) who co n- sidered this constraint in deriving maximum likelihoo d estima to rs for switching regr e ssions. More r ecently , Stephens (19 9 4) used it in a hierarchical Bayesian setting. In our setting, accor ding to the likelihoo d of (4) , a contin uity constraint o n the pr e- and po st-contact force-p osition curves a t x = x γ implies that β 10 + β 11 x γ + · · · + β 1 d 1 x d 1 γ = β 20 , (7) where β 1 = ( β 10 , β 11 , . . . , β 1 d 1 ) ′ and β 2 = ( β 20 , β 21 , . . . , β 2 d 2 ) ′ denote the vectors of pre- and p ost-contact regres sion co efficients, r esp e ctively . Higher-order smo othness ca n also be imp osed: we say that the for ce-p osition cur ve is s times co nt inuously different iable at x γ if the s th-or der deriv a tives of the pr e- and p o st-contact curves meet at x γ , with (7) corres p onding to the case s = 0 . On the other hand, if X 2 ,γ were a function of the p ositio n x , rather than the displacement x − x γ , then the contin uit y constraint would b eco me: d 1 X i =0 β 1 i x i γ = d 2 X j =0 β 2 j x j γ . (8) Either contin uit y constra in t implies that the likeliho o d function is nonline ar in the contact po in t x γ ; enfor c ing more degree s of smo othness at the changepoint serves to exac e r bate the nonlinearity and ma kes the design of efficien t inference metho ds incr easingly difficult. Stephens (1994) imp osed (8) in a Bayesian switching r egressio ns setting and pr op osed a rejection sampling step within a Gibbs sampler to addr ess the resultant nonlinea rity . Later, in Section 4, we describe a mor e efficien t appr oach that can be applied when either (7) or (8 ) (or higher-o rder analo gues) ar e enforced. 3.4. Changepoint estimation and contact point determination in the literatu re As demonstrated ab ov e, infer ence for material indentation data is well matched to cla ssical statistical frameworks for changepoint estimation. Indep endent of the s pecifics of our con- tact p oint pr oblem, the last ha lf century has seen a v a st b o dy of work in this area . Sequential and fixed-sample- s ize v arieties have b een considered fro m bo th cla ssical and Bayesian view- po in ts, with numerous parametr ic and nonpar ametric approaches prop osed. W e r efer the int erested rea der to several excellent s urveys, including tho se by Hinkley et al. (1980 ), Z a - cks (198 3), W o lfe and Sc hech tman (1984 ), Carlin et a l. (1992), and Lai (19 95). Some o f the ear liest work on maximum likelihoo d estimation of a single changep oint b et ween t wo po lynomial r egimes was done by Qua ndt (1958) and Robison (1964). 8 D . Rudoy , S. G. Y uen, R. D . Howe and P . J. W olfe Historically , Chernoff a nd Za cks (1964) were among the first to consider a pa rametric Bay es ian approa ch to changepo in t estimation. Cha ngep oints aris ing sp ecifically in linear mo dels have b een treated by many autho r s, including Ba con and W a tts (1 971), F er reira (1975), Smith (1 9 75), Choy and Bro emling (1980), Smith (19 80), and Menzefrick e (1981). The intro duction of Marko v chain Mo n te Ca rlo metho ds has led to more sophisticated hi- erarchical Bayesian mo dels for changep oint problems, b eginning with the semi-co njugate approach taken by Car lin et al. (199 2), in which the prior v ar ia nce of the reg ression co- efficients is left unscaled by the noise v ar iance. Adv a nces in trans-dimensiona l simulation metho ds have rekindled interest in m ultiple changep oint pro blems, a s discussed by Stephens (1994), Punsk ay a et al. (20 02) a nd F e a rnhead (2 006), among others. In the co n text o f mater ial indentation, how ever, exis ting a pproaches to contact p o int de- termination do not make us e of changepo int estimation metho dology . In fact, current meth- o ds are erro r-prone and lab or intensive—ev en co ns isting o f vis ual insp ection and manual thresholding (Lin et al., 20 0 7a). How ever, as describ ed ear lier, the many sources of v ar iabil- it y in indentation data imply that one cannot alw ays simply pr o ceed “ b y ey e.” More ov er , in the context of atomic force micr oscopy , mos t exp eriments aiming to character iz e c e ll stiffness, for example, emplo y m ultiple, re p e ate d indentations at different spa tial lo c a tions. These requirements have motiv a ted a more recent desire for effective, high- throughput au- tomate d techniques, a s deta iled in (Lin et al., 20 07a). Int erpreted in a statistical co ntext, the pro cedures thus far adopted by practitioners fall under the gener a l categ ory of likelihoo d fitting. Rotsch et al. (1999) s uggest simply to take tw o p oints in the p os t- c ontact data and solve for E and γ using the appr opriate Hertzian mo del; how ever, the resultant estima te of Y oung’s mo dulus dep ends strongly o n the inden ta tion depth o f the selected p o int s (Costa, 2004). Costa et a l. (20 0 6) pr op ose to minimize the mean-squa r ed er ror of a linear pre-co n tact and quadratic p ost-contact fit to the indentation data, though under the assumption of equal pre - a nd p ost-c ontact v ariance s . None of the existing approaches adopted b y practitioner s, how ever, pr ovides any mea ns of quantif ying uncertaint y in changepoint estimation—an impo rtant considera tion in prac- tice, since measur ement erro rs ca n b e la rge relative to the reaction for ce of the pro be d mate- rial, and consequently may r esult in p o or p oint estimates (Crick and Yin, 2007). Moreov er, such a pproaches fail to capture necessary physical constraints of the mater ial indentation problem, such as the smo othness cons traints describ ed in Section 3.3 ab ove. Such shor t- comings provide strong mo tiv ation for the hier archical model developed a bove, as well as the robust and automated fitting pro cedur es we describ e next. 4. P osterior inference via Marko v chain Monte Carlo The hierarchical Bayesian mo deling framework intro duced ab ove fea tures a large num b er o f unknowns, with constraints o n certain parameters precluding clo sed-form expressio ns fo r the marginal poster iors of in ter est. These consider ations sugg est a simulation-based a pproach to inference; indeed, it is by now s tandard to us e Mar ko v chain Monte Carlo metho ds to draw samples fro m the p osterior in suc h cases. Though widely a v ailable so ftw are pack ag es for Gibbs sampling are adequate for inference in certa in hiera rchical Bay esian settings, the complexity o f the conditional distributions we obtain here (after imp osing constraints and int egra ting out par ameters whenever p oss ible) necessitates explicit a lgorithmic deriv ations on a case- by-case basis. T o this end, w e build upo n the approaches o f Carlin et al. (19 92) and Stephens (1994), and employ Metrop olis-within-Gibbs tec hniques to draw s amples from the p osterior of (6) as well as under the smo othness constra in ts o f Section 3.3 . Bay esian changepoin t an alysis f or atomic force microscopy and soft mater ial indenta tion 9 4.1. Metropoliz ed Gibbs Samplers and V ar iance Reduction The selection of conjugate priors in our mo del allows nuisance pa rameters to b e integrated out, in orde r to reduce the v a riance of the resultant estimator s . F ollowing standa rd manip- ulations, we marginalize over the pr e- and po st-contact regressio n coe fficie n ts β 1 and β 2 , resp ectively . This yields the following ma rginal po sterior probability distribution: p ( γ , σ 2 1 , σ 2 2 , b 0 | y ; ψ ) ∝ σ − 2( a 0 − 1) 1 e − b 0 /σ 2 1 · σ − 2( a 0 − 1) 2 e − b 0 /σ 2 2 · b κ − 1 0 e − b 0 /η · ( | Σ γ | | Σ | | A γ | ) − 1 2 exp  − 1 2  y ′ Σ − 1 γ y + µ ′ Σ − 1 µ − b γ ′ A − 1 γ b γ  , (9) where A γ , X γ ′ Σ − 1 γ X γ + Σ − 1 Λ ∈ R p × p is blo ck-diagonal and b γ , ΛΣ − 1 µ + X γ ′ Σ − 1 γ y ∈ R p × 1 . The ma rginal p osterio r of (9) factors int o a Gamma density in b 0 , and inv er se-Gamma densities in σ 2 1 and σ 2 2 by wa y of the following partitions of A γ and b γ : A γ =  A 1 ,γ 0 0 A 2 ,γ  , A 1 ,γ ∈ R ( d 1 +1) ×⌊ γ ⌋ A 2 ,γ ∈ R ( d 2 +1) × ( n −⌊ γ ⌋ ) ; b γ =  b 1 ,γ b 2 ,γ  , b 1 ,γ ∈ R ⌊ γ ⌋× 1 b 2 ,γ ∈ R ( n −⌊ γ ⌋ ) × 1 . (10) The e xpressions of (9) and (10 ) lead to the following Gibbs sampler : Algorithm 1 Gibbs sampler for changep oint estimation (a) Draw γ ∼ p ( γ | σ 2 1 , σ 2 2 , b 0 , y ; ψ ) a ccording to (9); (b) Draw σ 2 1 ∼ I G  a 0 + 1 2 ⌊ γ ⌋ , b 0 + 1 2 ( y 1 ,γ ′ y 1 ,γ + µ 1 ′ µ 1 − b 1 ,γ ′ A − 1 1 ,γ b 1 ,γ )  ; (c) Draw σ 2 2 ∼ I G  a 0 + 1 2 ( n − ⌊ γ ⌋ ) , b 0 + 1 2 ( y 2 ,γ ′ y 2 ,γ + µ 2 ′ µ 2 − b 2 ,γ ′ A − 1 2 ,γ b 2 ,γ )  ; (d) Draw b 0 ∼ G  κ, η − 1 + σ − 2 1 + σ − 2 2  . T o simulate from the conditio na l distribution of γ , we employ as a Metrop olis -within- Gibbs step a mixture of a lo cal r andom walk mov e with an indep endent Metrop olis step in which the prop osa l density is a p oint wis e ev a luation of (9) on the grid 1 , 2 , . . . , n of inden ter lo cation indices. It is also p ossible to integrate out b oth noise v ariances (or the hyperpara meter b 0 ). In this case, additional Metr op olis steps ar e req uir ed, as the resulting conditional density of b 0 is nonstandar d. Our s imulation studies confirm that these v aria nts exhibit le s s Mont e Carlo v ariatio n than a Gibbs sampler based on the full p osterio r of (6). 4.2. P osterior inference in the presence of smoothness constraints Bearing in mind the underlying ph y sics of soft materials, a sufficien tly high sampling rate of the for c e senso r relative to the s p eed of the indenter may yield data consistent with a smo othness as sumption of a given order s . In this s etting, our inferential pro cedure may be mo dified accor dingly to a ppropriately ta ke this into acco un t. Stephens (1 994) c o nsiders contin uity-constrained switching linea r regr e ssions, and uses re jectio n s ampling to draw fro m the conditional distribution o f the changep oint γ given the remaining mo del par a meters. A more effectiv e pr o cedure, ho wev er , is to transform the data s uch that all but one o f the regres s ion co efficients may b e int egra ted out; this v aria nce r eduction yie lds a Gibbs sampler analogo us to Algor ithm 1 which we detail b elow. In fact, it is p ossible to der ive such an algorithm for an y v alue of s ≥ 0, though for ea se of pres en tation we first descr ibe the case s = 0, whereup on only contin uity is enforced. This approach to v a riance reduction in the presence of smo othness cons tr aints pr o ceeds as follows: define e β 1 = β 1 and, without lo ss of genera lit y , let e β 2 contain the last d 2 elements 10 D . Rud oy , S. G . Y uen , R. D . Howe and P . J . Wolfe of β 2 , so that e β = ( e β ′ 1 , e β ′ 2 ) ′ contains a ll p = d 1 + d 2 + 1 indepe nden t co efficient s. Via the contin uity constr aint of (7), de fine a line ar transfor mation T γ of e β to β as follows: β = T γ e β ; T γ ,   I ( d 1 +1) × ( d 1 +1) 0 ( d 1 +1) × d 2 c 1 × ( d 1 +1) 0 1 × d 2 0 d 2 × ( d 1 +1) I d 2 × d 2   , (11) where I m × m is the m × m identit y matrix, 0 m × n is the m × n matrix o f zeros , and c 1 × ( d 1 +1) , (1 , x γ , x 2 γ , . . . , x d 1 γ ). The choice of which of the d 1 + 2 regr ession co efficients to s e lect a s the depe ndent v ar ia ble is made without loss of genera lity , since a tr ansformation similar to (11 ) can be defined for every such choice, as well as for the contin uit y constraint of (8 ). Since one o f the reg ression co efficients is now a deterministic function of the others and the changep oint, we place a prior directly on e β rather than o n β . W e assume e β 1 ∼ N ( µ 1 , σ 2 1 Λ − 1 1 ) and e β 2 ∼ N ( e µ 2 , σ 2 2 e Λ − 1 2 ), with e Λ 2 a d 2 × d 2 diagonal po sitive definite matrix. Using the transforma tion T γ of (11), we obtain in analog y to (9) the full p oster ior p ( γ , e β , σ 2 1 , σ 2 2 , b 0 | y ; ψ ) ∝ σ − 2( a 0 − 1) 1 e − b 0 /σ 2 1 · σ − 2( a 0 − 1) 2 e − b 0 /σ 2 2 · b κ − 1 0 e − b 0 /η · ( | Σ γ | | e Σ | ) − 1 2 exp h − 1 2 n ( y − f X γ e β ) ′ Σ − 1 γ ( y − f X γ e β ) + ( e β − e µ ) ′ e Σ − 1 e Λ ( e β − e µ ) oi , (12) where f X γ = X γ T γ ∈ R n × p and, in analog y to the quantities ( Σ , µ , Λ ), we define e Σ , diag( σ 2 1 1 d 1 +1 , σ 2 2 1 d 2 ) ∈ R p × p , e µ = ( µ ′ 1 , e µ ′ 2 ) ′ ∈ R p × 1 , and e Λ = diag( Λ 1 , e Λ 2 ) ∈ R p × p . The transforma tion T γ makes it p oss ible to integrate out the r e gression co efficients e β using standar d ma nipulations. Indeed, intro ducing the terms ˜ A γ , f X ′ γ Σ − 1 γ f X γ + e Σ − 1 e Λ ∈ R p × p and ˜ b γ , e Λ e Σ − 1 e µ + f X ′ γ Σ − 1 γ y ∈ R p × 1 as b efore, we obtain the ma r ginal p os terior p ( γ , σ 2 1 , σ 2 2 , b 0 | y ; ψ ) ∝ σ − 2( a 0 − 1) 1 e − b 0 /σ 2 1 · σ − 2( a 0 − 1) 2 e − b 0 /σ 2 2 · b κ − 1 0 e − b 0 /η ·  | Σ γ | | e Σ | | A γ |  − 1 2 exp n − 1 2 ( y ′ Σ − 1 γ y + e µ ′ e Σ − 1 e µ − e b ′ γ e A − 1 γ e b γ ) o . (13) It is straig htforward to g eneralize this notion to any s ∈ {− 1 , 0 , . . . , d 1 + d 2 } ; a prior is put on d 1 + d 2 − s + 1 reg ression co efficients, and a transforma tion T γ analogo us to (11) defined. As no ted previously , the smo othness constraint of (7) introduces dep endence a mong the pr e- and p ost-changep oint regr e ssion co efficients. In contrast to the mar ginal p os terior distribution of (9) derived earlie r fo r the unconstrained case, enforcement of (7) prec ludes int egra ting out the asso ciated nois e v ariances σ 2 1 and σ 2 2 . In the former ca se, the blo ck- diagonal structur e of X γ (and there fo re of A γ ) implies that the induced joint distr ibution on the v a riances is separ able. How ever, in the latter case o f (13 ), f X γ is not block diagonal— owing to the a ction o f T γ —and hence neither is e A γ . Therefore, the v a riances σ 2 1 and σ 2 2 are no lo ng er conditionally indep endent, a nd their joint distribution do es not take the form of known genera lizations of the univ ariate Gamma distribution to the biv aria te cas e (Y ue et al., 2001). Consequently , only the conditional distr ibution of the hyper parameter b 0 is in standard form, a nd sim ulation from (13) pro ceeds with all other v aria bles drawn using Metrop olis-Hasting s (MH) steps, as shown in Algorithm 2 b elow. In co nt rast to the ca se of Algo rithm 1, where a mixture k ernel was employ ed purely fo r r easons of computational efficiency , we emphasize here that such a mov e is in fact re quir e d to sample fro m the full suppo rt (1 , n ) of the changepoint index, otherwise mixing of the underlying c hain is p o or. Bay esian changepoin t an alysis f or ato mic f orce microscopy and soft mater ial indentatio n 11 As befor e , the mixture kernel consisted o f a lo cal random walk mo ve and an indep endent global mov e drawing from a discr ete distribution derived as a p oint wise ev a luation of (13) on the int egers 1 , 2 , . . . , n . The coupling of noise v ariances suggests a joint MH mov e, but this require s sp ecification o f a prop osal cov ariance; in simulations we found sepa rate Normal random w alk mov es for ln( σ 2 1 ) and ln( σ 2 2 ) to be adequate. Algorithm 2 Smo othness-co nstrained Gibbs sampler for changepoint estimation (a) Draw γ ∼ p ( γ | σ 2 1 , σ 2 2 , b 0 , y ; ψ ) accor ding to (13) using a Metrop olis- within-Gibbs step; (b) Draw σ 2 1 ∼ p ( σ 2 1 | γ , σ 2 2 , b 0 , y ; ψ ) likewise; (c) Draw σ 2 2 ∼ p ( σ 2 2 | γ , σ 2 1 , b 0 , y ; ψ ) likewise; (d) Draw b 0 ∼ p ( b 0 | γ , σ 2 1 , σ 2 2 , y ; ψ ) = G  κ, η − 1 + σ − 2 1 + σ − 2 2  . W e re- e mpha size that in our exp erience, int egra ting out the regression co efficient s is essential in o rder to obtain a Gibbs sampler with fav or a ble mixing prop erties. In particular , a sampler drawing for each par ameter of (12) in the presenc e of smo othness constraints is severely restr ic ted in its ability to explo re the state space if one- co ordinate-a t-a-time up dates are e mployed. When s contin uous deriv atives are enforced at the changep oint, the likeliho o d function bec o mes mor e nonlinear in changepoint x γ as s increases. Thus, as the induced constraint set b ecomes more no nlinear, lo cal mov es o n the scale o f the regr e s sion co efficients themselves will lie far from it—leading to small acceptance probabilities. T o overcome these problems, one may des ign a high-dimensional MH move to up date all the smo othness - constrained regre s sion co efficie nts join tly with the changepoint; how ever, a unique mov e m ust b e desig ned for each s to b e considere d. In contrast, integrating o ut the r egressio n co efficients obviates this need by ha ndling such constra int s for a ll s sim ultaneously . 5. Experimental v alidation In o rder to exp e rimentally v alidate Algorithms 1 and 2 prior to their application in the setting of atomic force micro scopy , we designed and p erformed tw o sets of sp ecial macr o - scale inden tation ex per imen ts in which precise contact p oint identification was made po ssible by the use of an imp edance-meas ur ing electro de mounted on the indenter. While this direct measurement pro cedure is pr ecluded in the v ast ma jority of bioma terials ex per imen ts inv olving cells and tissues, as such samples a re submerg ed in an aqueous solution, we were able to measure the true contact p o int for exp eriments inv olving resp ectively cantilev er bending a nd silicone indentation, as descr ibed below. In a ddition, we co nducted a n umber of simulation studies to ch ara c terize uncertaint y in changepoint estimation a s a function of v a r ious mo del parameter s . O n the basis of these simulation studies and subsequen t experimental v alida tion, we found b oth algorithms to be insensitive to the exact choice of h yp erparameter s ψ , a nd hence r etained the fo llowing settings throughout: µ = 0 , Λ = 10 − 5 I p × p , a 0 = 2 , κ = 1, and η = 10 − 2 , with p b eing the total num b er o f r egressio n co efficients. W e set the pre-contact p olynomia l degr ee d 1 = 1 throughout, based on the discuss io n in Section 2 s uppo r ting the as s umption o f a linear pre- contact regime. When smo o thness constraints were us ed, we retained identical pa rameter settings, and decremented the v alue of p appro priately . All p osterior distributions were obtained by running the appropr iate Gibbs samplers for 50 0 00 iterations, a nd discarding the fir st 5 000 samples. Co n vergence was assessed using standa rd metho ds co nfirming that increasing the n umber o f Gibbs iter ations did not a ppreciably change the re s ulting inference. 12 D . Rud oy , S. G . Y uen , R. D . Howe and P . J . Wolfe 5.1. V alidation of changepoint inference via cantile ver bending W e fir st p erformed several trials of an expe riment wher eby a steel cantilev er was b ent through application of a down ward unia xial force. Here, the measured for c e F is expected to c hange linearly with displacement δ = x − x γ according to Ho o ke’s Law, with γ ∈ (1 , n ) representing the contact p oint index. Repres ent ative data shown in Figur e 3 were obtained using a T estBe nch TM Series system with a high-fidelit y linear a ctuator (Bos e Corp or a tion EnduraTEC Systems Group, Minnetonk a, Minneso ta, USA), which mov ed a cylindr ical inden ter int o a ca nt ilevered piece of FSS-0 5/8-1 2 spring steel mea s uring approximately 1 . 27 cm × 2 . 54 cm (ca rb on co nt ent 0 . 9 –1 . 05%; Small Parts, Inc., Miramar , Flo rida, USA) at a spe e d of 10 mm/s. According to the impeda nce measur e men t technique describ ed ab ov e, the contact po in t index γ was determined to co rresp ond to p os itio n index 48 of the indent er, co rresp onding to a co n tact po int x γ ∈ [ − 1 . 3 58 , − 1 . 262] mm. Because of the hardness of steel a nd the sp eed of the indenter, we did not ma ke a s mo othness assumption, and th us used the Gibbs sampler of Algorithm 1 to ev aluate the efficacy of our a pproach, with p = 4 based on the linea r po st-contact r egime implied by Ho oke’s law. A full 1 00% of p os ter ior v alue s for γ a fter a 10% burn- in p or tion took the v alue 48, indicating correct detection of the changep oint. Fig ur e 3 sho ws the re s ults of the curve fitting pro cedure,with similar results obtained for v ary ing indenter spee ds . The minimum mean-squar e error (MMSE) estimate o f Y o ung’s mo dulus was determined to b e 215 . 3 GPa, with an a sso ciated 95% poster ior in terv al of [21 4 . 0 , 21 6 . 6]. By comparison, the range of v alues of Y oung ’s mo dulus for steel with similar carb on conten t is r epo rted in the litera ture to b e 210 ± 12 . 6 GPa (Kala and Kala, 2005). Despite its s imple desig n, this exp e r iment is not far remov ed from practice; a nearly identical pro cedure was employ ed b y W ong et al. (1997) to prob e the mec hanical prop erties of silicone- c arbide nanoro ds, with ea ch na noro d pinned on a substrate and sub jected to a b ending for ce a long the unpinned axis. 5.2. Analysis of silicone indentation data While changes in slop e a t the contact p oint tend to be more pronounced in harder materials such as stee l, the change in measured force is typically smo other in softer ma ter ials such as cells and tissue s , making the contact p oint harder to detect. In earlier work (Y uen et al., 2007) we detaile d an indentation exp eriment using a soft silicone s ample (Aquaflex; Park er Lab ora tories, F airfield, New Jerse y , USA), c hosen both for its similar it y to h uman tissues used in ma terial indentation studies (Chen et al., 1996), a nd b ecaus e it enables direct contact p oint determination via a n imp edance- measuring electro de. T en trials of this exp eriment were conducted, using a sa mple roughly 20 mm in depth, with a maximum indentation of approximately 8 mm; a t ypica l force-displac e ment curve was shown ear lier in Figure 2. A hemispher ical metal indenter of radius R = 87 . 5 mm com- pressed the sample at 10 mm/s, and the resulting forces were measured appr oximately every 10 µ m to yield ≅ 96 0 fo rce-displacement data p oints. Both the unconstraine d ( s = − 1) and contin uity-constrained ( s = 0) mo dels were fitted in this setting, b y wa y of Algo rithms 1 and 2 r e s pec tively . F or a s pherically-tipp ed indenter, the Hertzian mo del of (2 ) indicates that force is prop ortional to ( x − x γ ) 3 / 2 , and hence we employ ed a po st-contact design ma - trix with only the fractional p ow er 3 / 2 . In this regime, w e hav e that p = 4 for Algorithm 1 and p = 3 for Alg orithm 2 . Since the initial distance from the indenter to the sa mple was approximately known, a unifor m prio r on γ ∈ [125 , 2 50] was assumed. F or e ach of the ten data sets collected, the first n = 450 data po int s were taken to represent a conser v ative estimate of op eration within the small-deformatio n regime, and Bay esian changepoin t an alysis f or ato mic f orce microscopy and soft mater ial indentatio n 13 −4 −2 0 0 0.5 1 1.5 Displacement (mm) Force (N) Spring Data −4 −2 0 −0.02 −0.01 0 0.01 0.02 Displacement (mm) Force (N) Residual to MMSE Fit −2.5 −2 −1.5 −1 −0.1 0 0.1 0.2 0.3 0.4 Displacement (mm) Force (N) Spring Data 214 216 218 0 0.05 0.1 0.15 0.2 Young’s Modulus (GPa) Sample Proportion Marginal Posterior of Young’s Modulus Fig. 3. Inference for the can tilev e r experiment of Section 5.1. The force-displacement da ta and poster ior me an re constructi ons obtai ned via Algor ithm 1 are sh own (top left) with a close-up view near the changepoin t (botto m left) and the associated re siduals (top right). T he margin al p osterio r of Y oung’s mo dulus (bottom right ) is sho wn with its mean (black li ne) and 95% i nter val (grey lines). the Gibbs sa mplers of Algorithms 1 and 2 were each run on these data. Res ults for b oth the cases are summar ized in T able 1, a long with the exp erimentally-determined contact po in t, which v ar ied from trial to trial due to vis co elastic effects. Indeed, given that the spacing b etw een data p oints is 0 . 01 mm on a verage, it may b e deduced fr om T able 1 that the average error of 0 . 8– 1 % acro ss tria ls corresp onds to 8–10 sampled data p oints. Marginal contact p oint po sterior distributions for b oth the unco nstrained and contin uity- constrained ca ses a re summarized, pair wise by exp eriment, in the box plo ts of Figure 4 . These a re seen to b e notably more diffuse in the former case (left) tha n the latter (right) ; this is co nsistent with the softness of the silicone sample under study , which makes it difficult to r eject a priori the p oss ibility o f contin uit y at the changep oint. Absent this assumption, the ex per imen tally determined co nt act po in t lies within the 95% p osterior int erv al for ea ch of the ten tria ls. Once a co n tinuit y constraint is imp osed, the marg inal po sterior distributions of the changep oints tighten noticeably; this is consistent with our exp ectation that constra ining the mo del r educes the num b er of high-likelihoo d fits. A more subtle po int ca n also b e deduced from the s light yet consistent left-shift across 95% p oster ior interv a ls under the co n tinuit y-co ns trained regime re la tive to the unco n- strained case. Observe the fourth row of T able 1, whic h rep o rts the v a lues obtained up on extrap olating pre- and p ost-contact MMSE curve fits to their meeting p o int in this latter 14 D . Rud oy , S. G . Y uen , R. D . Howe and P . J . Wolfe T a ble 1. Con tact po int estimates and associat ed er rors ( mm, %) f or the ten sili cone tr ials of Sec- tion 5.2, based on unconstraine d (U) and con tinuity-constrained (C) m odel s, with f or ce-response data sampled ev er y 0.01 mm on av erage. Y oung’s modulu s esti mates b E a nd 9 5% po sterior inter vals for the unconstrained case are also shown, along with av erages ov er all ten tr ials where appropr iate. T rial N o. 1 2 3 4 5 6 7 8 9 10 Avg. T ruth ( x γ ) 5 . 52 5 . 49 5 . 47 5 . 44 5 . 42 5 . 38 5 . 39 5 . 48 5 . 42 5 . 38 – MMSE (U) 5 . 42 5 . 44 5 . 42 5 . 43 5 . 47 5 . 39 5 . 46 5 . 35 5 . 41 5 . 3 9 – MMSE (C) 5 . 40 5 . 43 5 . 41 5 . 39 5 . 46 5 . 37 5 . 41 5 . 34 5 . 34 5 . 3 9 – Extrap. (U) 5 . 39 5 . 42 5 . 40 5 . 37 5 . 42 5 . 36 5 . 40 5 . 32 5 . 34 5 . 3 8 – % Err. (U ) − 1 . 78 − 0 . 79 − 0 . 83 − 0 . 08 0 . 8 8 0 . 19 1 . 40 − 2 . 23 − 0 . 10 0 . 11 0 . 83 % Err. (C) − 2 . 06 − 1 . 02 − 1 . 13 − 0 . 95 0 . 3 1 − 0 . 25 0 . 47 − 2 . 44 − 1 . 41 0 . 03 1 . 00 2.5% 16 . 5 16 . 5 16 . 4 16 . 4 16 . 5 16 . 1 16 . 6 16 . 1 16 . 2 16 . 4 – b E ( kPa ) 17 . 2 17 . 3 17 . 2 17 . 0 17 . 2 16 . 8 17 . 4 16 . 8 17 . 0 16 . 9 17 . 1 97.5% 17 . 9 18 . 4 17 . 9 17 . 6 18 . 0 17 . 6 18 . 1 17 . 7 17 . 8 17 . 9 – 1 2 3 4 5 6 7 8 9 10 5.2 5.25 5.3 5.35 5.4 5.45 5.5 5.55 5.6 5.65 5.7 Contact Point (mm) Trial Number Fig. 4. B ox plots of conta ct point margi nal po sterio r distri butions in ea ch of ten silico ne i ndenta tion tria ls, shown side-by-side for the unconstrain ed (left) and continuity-con strained (right) models, w ith centers of th e grey squares indi cating tr ue con tact point values. Note the consistent decre ase in poster ior v ar iance and slight downward shift of the poster ior under the assumption of continuity . case. The increa se in for ce after contact implies that these v alues will always lie b elow the directly inferred co nt act p oint, as indeed they do. Nev ertheless, enforcing contin uity results in only a small increa se in ov erall contact p oint estimation err or, and pr o duces what practi- tioners might judge to b e more physically feasible curve fits. While indep endent verification of the Y o ung’s mo dulus is not av ailable for the silicone sample used in our study , the quality of our c ontact p oint estimates r elative to a known g round truth leads to high confidence in the inferred v alues of Y o ung’s mo dulus. Overall, the ability to obtain inferential re sults and accompanying uncertaint y quan- tification, as e xemplified by the cantilev er b ending and s ilicone indentation e xpe r iment s, Bay esian changepoin t an alysis f or ato mic f orce microscopy and soft mater ial indentatio n 15 represents a significant improvemen t upo n c ur rent methods ; a more in-depth compar ison to the method of Costa e t a l. (2006) is provided in o ur ear lier work (Y uen et al., 2007 ). In particula r, the latter exper imen t demonstrates that reliable es timates of soft mater ia l prop erties can b e obta ined even in the presence of measurement erro r—a reg ime applicable to many AFM studies of cells a nd tissues , as w e now descr ibe . 6. Inference in the setting of atomic f orce micr oscopy Having v alidated our algorithms on a mac roscopic scale, we tur n to analyz ing cellular bio- materials data collected using atomic fo rce microscopy techniques. In contrast to our earlier examples, no direc t exper imen tal verification is av ailable in this case, though we no te that our resultant estimates of Y o ung’s modulus are considered plausible by ex per imen talists (So crate and Suresh Labs, p ersonal co mmunications, 2 0 08). It is widely b elieved that cell biomechanics can shed light on v ario us disea s es of imp ort, and hence a key r esearch ob jective following the adv ent of AFM technology has been to analyze quantitativ ely the mec ha nical prop erties of v a r ious cell t yp es. Indeed, n umerous pap ers over the past decade hav e linked stiffness and other related mechanical prop erties to cell malfunction and death (Cos ta, 2004 ); for example, the ability of car diac m yocytes in hear t muscle tissue to contract is intimately linked to their cytoskeletal structure and its influence on cellular mechanical r e s po nse (Lieb er et a l., 2004). Here, we are similarly motiv ated to under stand the s tiffness prop erties of neurona l a nd r ed blo o d c ells—currently a topic of intensiv e res earch in the bio mechanics and bio engineering co mm unities. 6.1. Indentation study of embryonic mouse cor tical neu rons W e fir st analyze d ex-vivo live mo use neurons, submerged in cell culture medium and rep eat- edly indented by an AFM (Asylum Resear ch, Santa Ba r bara, California , USA) equipp ed with a spher ically-tipp ed pro b e. The indentation of e ach neuron (So cr ate Lab, Mas- sach use tts Institute of T echnology) yielded appr oximately 700 force measurements, of which the fir st 50 0 were used in subsequent analysis in order to stay within the s ma ll-deformation regime. Such data are of w ide interest to ne ur oscientists and engineer s, as traumatic damage to neurons is hypothesized to b e re la ted to their mec hanical prop erties . As a spherical prob e tip was used for indent ing each cell, we employed a p ost-c ontact design matrix with only the fractional p ow er 3 / 2 , as in the case of our ear lier silicone ex- ample, and the contin uity-constrained sampler of Algorithm 2 . The results of a typical trial are shown in Fig ure 5; the pre- a nd p ost-c o nt act residuals were observed to b e white. The primary inferential quantit y o f interest in such ca ses is Y oung’s mo dulus, the pr inci- pal characterization of cell material stiffness in tro duced in Section 2.2 . Accor ding to the Hertzian mo del of (2) for a s pher ical indenter, the r egressio n co efficient corr esp onding to the ( x − x γ ) 3 / 2 term of the fitted mo del is prop ortio nal to Y oung’s mo dulus, with the constant of prop ortio nality a function of the given radius R = 1 0 µ m and Poisson’s r atio ν = 0 . 5. Thu s, we ca n obtain the p oster ior distribution o f Y oung’s mo dulus by appr opriately scaling the dis tribution of this regress io n co efficient, as shown in the b o ttom left panel o f Fig ure 5. W e r epo rt the MMSE estimate of Y o ung ’s mo dulus as b E = 5 30 . 4 Pa, and the corre s po nd- ing 95 % po sterior interv al a s [518 . 8 , 544 . 1] Pa. This estimate is in reasonable a greement with those previously rep orted for similar neurons (Lu et al., 200 6; Elkin et al., 20 07), with v ariability due primarily to differences in indentation sp eeds (So crate Lab, p erso nal communication, 2008). 16 D . Rud oy , S. G . Y uen , R. D . Howe and P . J . Wolfe 7 8 9 10 5 4.5 4 3.5 Displacement!( µ m) Force!(nN) Neuron!Data 7 8 9 10 0.1 0.05 0 0.05 0.1 0.15 0.2 0.25 Displacement!( µ m) Residual!(nN) Residual!to!MMSE!Fit 510 520 530 540 550 560 0 0.05 0.1 0.15 0.2 0.25 Marginal!Posterior!of!Young’s!Modulus Young’s!Modulus!(Pa) Sample!Proportion 9.39 9.4 9.41 9.42 9.43 0 0.05 0.1 0.15 0.2 0.25 Displacement!( µ m) Sample!Proportion Marginal!Posterior!of!Contact!Point Fig. 5. D ata co llected dur ing a n AFM i ndenta tion of a mouse n euron togeth er with the poster ior mean estimate of the und erlyin g regressions (top left) an d the in duced residua ls (top righ t). Also shown are marginal p osteri or di stributio ns of Y oung’ s m odul us (bo ttom left) and contact point (botto m ri ght), each ov er laid with the poster ior mean (blac k line) and 95% posterior inter val (dashed grey lines). 6.2. Indentation study of red b lood cells The mechanics of r ed blo o d c e lls have also been extensively studied using atomic force microscopy (Radmacher et al., 19 96). In this vein, we next ana lyzed data fr om ex-v ivo live hu man erythro cytes (red blo o d cells, Sures h Lab, Mas sach usetts Institute o f T echnology) submerged in a cell culture medium, indented by an AFM (Asylum Resear ch ) equipp ed with a pyramidally-tippe d prob e. The final 800 of approximately 8 5 00 data p oints were discar ded prior to a nalysis, as they clear ly lay outside the small-deformatio n r egime. Relative to the neuron AFM data c onsidered in Section 6.1, the s ampling rate of re s istive force here is hig h, and c o nsequently it is feasible to enfor ce co nt inuit y ( s = 0) a t the changep oint. Algo rithm 2 was ther efore employ e d, with r esults from a t ypical trial shown in Fig ure 6. In the case a t hand, the r e gression co efficien t β 12 corres p onding to ( x − x γ ) 2 is prop or- tional to the Y oung’s mo dulus as p e r (2), with the consta n t of pro po rtionality dep ending on the inden ter tip angle 2 φ = 70 ◦ and Poisson’s ratio ν = 0 . 5 , and, as we discus s below, no linear p ost-co n tact regres sion ter m was included. The inferred distribution o f Y oung’s mo dulus may b e obtained b y a ppropriately tra nsforming the marg ina l p osterior of β 12 , as detailed in Section 3.1 , and is shown in the b ottom left panel o f Figure 6. The resulta nt MMSE estimate of b E = 25 . 3 Pa and corr esp onding p oster ior interv al of [16 . 0 , 34 . 9] Pa were Bay esian changepoin t an alysis f or ato mic f orce microscopy and soft mater ial indentatio n 17 •5 •4 •3 •2 •1 •60 •40 •20 0 20 40 60 Displacement ( µ m) Force (pN) Red Blood Cell Data −5 −4 −3 −2 −1 −60 −40 −20 0 20 40 60 Residual (pN) Displacement ( µ m) Residual to MMSE Fit 0 4 8 12 16 20 24 28 32 36 40 44 48 0 0.1 0.2 0.3 0.4 Marginal Posterior of Young’s Modulus Young’s Modulus (Pa) Sample Proportion 9.5 10 10.5 11 11.5 12 0 0.05 0.1 0.15 0.2 0.25 Noise Std. Dev. (pN) Sample Proportion Marginal Posterior of Noise Std. Dev. σ 1 σ 2 Fig. 6. Data collected dur ing an AFM inde ntatio n of a human red bloo d ce ll toge ther with th e po sterio r mean estimate of the un derl ying regression s (top lef t) a nd the ind uced residu als (to p righ t). Also shown are margina l po sterio r distr ibution s Y oung’s modulus (b ottom left) and of σ 2 1 and σ 2 2 (bottom righ t), the former ov erla id with the posterio r m ean and 95% poster ior inter v al. confirmed to be consistent with v arious experimental ass umptions (Suresh La b, p ersonal communication, 2008). F urther, the pre - and p ost-co n tact error v aria nces are determined to be unequal, a s shown in the b ottom rig h t panel of Figure 6. Note that in the absenc e of a p ost-contact drift, enforcing contin uous differentiabilit y at the changep oint ( s = 1) constrains the pre-contact linear fit to hav e zero slop e, and is inconsistent with the pre-contact drift clear ly visible in the top left pa nel of Figure 6. On the other hand, if the p ost- c ontact p o lynomial were to include a linea r drift term, then enforc ing contin uous differentiabilit y would imply that the pr e- and p ost-co n tact drifts are identical. Though this is app ealing from a mo deling viewp oint, as it eliminates an additiona l fr e e parameter, practitioners lack e v idence for such an equality . Moreov er , in our exp eriments its inclusion had no appreciable effect on the infere nc e of Y oung ’s mo dulus, and so we did not include a po st-contact drift ter m in o ur final analysis of these data. As in the case of our earlier exper imen ts, we co mpared our approach to the likelihoo d metho d of Co sta et al. (20 06), which yielde d a n estimate of x γ shifted to the rig h t by more than 1 000 da ta p oints r e la tive to tha t shown in Fig ur e 6 . The corr esp onding estimate of Y oung’s mo dulus in turn was found to b e 34 . 8 Pa—an increase o f 37 . 5% relative to the MMSE p oint estimate o f 25 . 3 Pa, and close to the upp er bo undary o f our estimated po sterior int erv al. While in this exp erimental setting one ca nno t conclude that either e s timate is 18 D . Rud oy , S. G . Y uen , R. D . Howe and P . J . Wolfe sup e rior to the o ther, we note that the difference b etw een them can b e attributed in part to our mo del’s incorp oratio n o f differing pre- and p ost-co n tact err or v ariances, and smo othness constraints. Thus, o ne may view our infer e n tial pro cedures as b oth a for malization a nd a n extension of earlier likelihoo d-based approa ches developed by pr actitioners, enabling b oth robust, automated fitting pro cedur es and explicit uncertaint y quantification. 7. Discussion In this article w e have p osed the firs t r igorous formulation of—and solution to—the key inferential proble ms arising in a wide v ariety of material indentation systems and studies. In particular , pr a ctitioners in the materials science communit y to da te hav e la ck ed a ccurate, robust, and a utomated to ols for the estimation of mechanical prop er ties of s o ft materia ls at either macr o - or micro-sca le s (Lin et al., 200 7a; Crick a nd Yin, 200 7). A principal strength of our approa ch is its applicability to the analys is of bio materials da ta obta ined by indenting cells and tissues using atomic force micros copy; con tact p oint determination is even more difficult in this s etting, owing to the gr adual c hange of mea sured re s istive force that is a hallmark of soft mater ials. The Bayesian switching p olynomial reg ression mo del and asso ciated inferential pro cedure s w e have pro po sed provide a means bo th to determine the p oint at which the indenting prob e comes into contact with the sa mple, a nd to estimate the corresp onding materia l proper ties such as Y o ung’s mo dulus. In turn, its careful character iz ation holds op en the even tual pos sibilit y of new biomechanical testing pro cedures for disease (Costa, 20 04). Our parametr ic approach is s trongly motiv ated by—and ex ploits to full adv antage— the Hertzian mo dels gov erning the ph ysical behavior of linear - elastic materials underg oing small deformations . The Bay esia n para digm not o nly enables uncertaint y quantification, crucial in applica tions, but als o allows for the natural incor po ration o f physically-motiv ated smo othness co nstraints at the changep oint. Infer ence is realized through application o f carefully desig ned Markov chain Mo n te Ca rlo metho ds together with classical v a riance re- duction techniques. The r esultant a lgorithms hav e bee n shown here to b e b oth statistically and computationally e fficie n t as w ell as r obust to choice of hyper parameters over a wide range of e xamples, a nd are av ailable o nline for us e by practitioners. Indeed, the direct applicability of our metho ds prec ludes a n y need for data pre- pr o cessing prio r to analysis. Outside of the linear-elastic materials w e consider here, it is of int erest to apply the metho dology to visco ela stic materials (e.g., biop olymers) which return to their pre-contact state slowly over time (Lin et al., 200 7b). In such cases, the a mount of induced deformation depe nds not o nly on the indent er g e o metry , but also o n the rate o f inden tation. Another extension is to incorp or ate multiple spatia lly distributed changepo int s into our mo del, a key construct when atomic force micros c op es are us ed to indent repea tedly a s ample in order to characterize cell stiffness as a function of sur face lo cation (Geisse, 2009). Finally , a sequential estimation sc heme could b e of great use in surgica l r ob otics applications, wher e contact p oint determinatio n plays a key role in enabling tactile sensing —a s ub ject o f curr ent study b y the author s. Ac knowledgement s The a uthors would like to ackno wledge Kris tin Bernick, An thony Gamst, Hedde v an Hoo rn, Petr Jordan, and Thibault P revost for helpful discussions, and would esp ecially like to thank Bay esian changepoin t an alysis f or ato mic f orce microscopy and soft mater ial indentatio n 19 Simona So cra te and Subra Sures h for pr oviding access to data from atomic force micro scop e inden tation exp eriments on neurons and red blo o d c ells, res pectively . The first author is sp onsored b y the National Defense Science and Engineer ing Gradua te F ellowship. The second author is sp onsor e d by United States Na tional Institutes of Health Grant No. NIH R01 HL0736 4 7-01. The author s are gr ateful to anonymous revie w ers for many suggestio ns that have help ed to improv e the clarity of this ar ticle. References Bacon, D. W . and D. G. W atts (1971). Estimating the transition b et wee n tw o intersecting lines. Biometrika, 58, 525–534. Carlin, B. P ., A. E. Gelfand, and A. F. M. S m ith (1992). Hierarc h ical Bay esian analysis of change- p oin t problems. Appl. Statist., 41, 389–405. Chen, E. J., J. Nov akofs ki, W. K. Jenkins, and J. W. D Brien (1996). Youn g’s modu lu s measure- ments of soft tissues with app lication to elasticity imaging. IEEE T r ans. Ultr ason. F err o ele ct. F r e q. Contr ol, 43, 191–194. Chernoff, H. and S . Zacks (1964). Estimating th e current mean of a Normal distribution which is sub jected to changes in time. A nn. Math. Statist., 35, 999–10 18. Cho y , J. H. C. and L. D. Bro emling (1980). Some Bay esian inferences for a changing linear mo del. T e chnometrics, 22, 71–78. Costa, K. D . (2004). Single-cell elastography: Probing for d isease with the atomic force microscop e. Dise ase Markers, 19, 139–154. Costa, K. D., A. J. S im, and F. C.-P . Y in (2006). Non- Hertzian approach to analyzing mechanical prop erties of endoth elial cells prob ed by atomic force microscop y . ASM E J. Biome ch. Eng., 128, 176–184 . Cric k, S. L. and F. C.-P . Y in (2007) Assessing micromec hanical prop erties of cells with atomic force m icroscopy: Importance of the contact p oint. Biome ch. Mo del Me chanobiol., 6, 199–210. Dimitriadis, E. K., F. Hork ay , J. Maresca, B. Kachar, and R. S. Chadwic k (2002). Determination of elastic mo duli of th in lay ers of soft material using the atomic force microscope. Biophys. J., 82, 2798–28 10. Elkin, B. S., E. U. A zeloglu, K. D. Costa, and B. M. Morrison I I I (2007). Mec hanical heterogeneit y of the rat hipp o campus measured by atomic force microscope indentation. J. Neur otr auma, 24, 812–822 . F earnhead, P . (2006). Exact and efficien t Ba yesian in ference for multiple c hangep oint problems. Statist. C omput., 16, 203–213. F erreira, P . E. (1975). A Ba yesian analysis of a switching regression model: A kn o wn num b er of regimes. J. Am. Statist. Ass., 70, 370–374. Geisse, N. A., S. P . Sheehy , and K. K . Park er (2009). Control of my o cyte remo deling in vitro with engineered substrates. I n Vitr o Cel l.Dev.Biol.-Animal, 45 , 343–350. Gouldstone, A., N. Chollacoop, M. Dao, J. Li, A. M. Minor, and Y.-L. S hen (2007). Ind entation across size scales and disciplines: Recent developmen ts in exp erimentation and mo deling. Ac ta Materialia, 55, 4015–4039. Hinkley , D., P . Chapman, and G. Run ger (1980). Change p oint problems. T ec hnical rep ort, Universit y of Minnesota, Minneap olis. Hudson, D. J. (1966). Fitting segmen ted curve s whose join p oints hav e to b e estimated. J. Am. Statist. A ss., 61, 1097–1129. Kala, J. and Z. K ala (2005). I nfluence of y ield strength v ariability o ver cross-section to steel b eam load-carrying capacity . Nonlin. Anal.: M o del. Contr ol, 10, 151–160. Lai, T. L. (1995). Sequential changep oint d etection in quality control and dyn amical systems. J. R. Statist. So c. B , 57, 613–658. Lieber, S., N. Aub ry , J. Pain, G. D iaz, S.- J. Kim, and S. V atner (2004). Aging increases stiffness of cardiac my ocy tes measured by atomic force microscop y nanoindentation. Am. J. Physiol.—He art Cir culat. Physiol., 287, 645–651 . Lin, D. C., E. K. Dimitriadis, and F. Hork a y (2007a). R obust strategies for automated AFM force cu rve analysis—I. Non-ad h esiv e indentation of soft inhomogeneous materials. J. Biome d. Eng., 129, 430–440. Lin, D. C., E. K. Dimitriadis, and F. Hork ay (2007b). R obust strategies for automated AFM force curve analysis—I I. Adhesion-influ enced indentation of soft, elastic materials. J. Biome d. Eng., 129, 904–912. 20 D . Rud oy , S. G . Y uen , R. D . Howe and P . J . Wolfe Lu, Y.-B., K . F ranze, G. Seifert, C. Steinhauser, F. Kirchoff, H. W olburg, J. Guck, P . Janmey , E.-Q. W ei, J. Kas, an d A. Reichen b ac h (2006). Visco elastic prop erties of individual glial cells and neurons in t h e CNS . Pr o c. Natl. A c ad. Sci . USA, 103, 17759–177 64. Menzefric ke, U . (1981). A Ba yesia n analysis of a change in the p recision of a sequ ence of indep end ent Normal random va riables at an unknown time p oint. Appl. Statist., 30, 141–146 . Punsk aya , E., C. Andrieu, A. Doucet, and W . J. Fitzgerald (2002). Bay esian cu rve fitting u sing MCMC with applications to signal segmentation. IEEE T r ans. Signal Pr o c ess., 50, 747–758. Quandt, R. (1958). The estimation of the parameters of a linear regression system ob eying tw o separate regimes. J. Am. Statist. Ass., 53, 873–88 0. Radmacher, M., M. F ritz, C. K acher, J. Clev eland, and P . Hansma (1996). Measuring the visco elas- tic prop erties of human platelets with the atomic force microscope. Biophys. J., 70, 556–567. Robison, D. E. (1964). Estimates for the p oints of interse ction of t wo p olynomial regressions. J. Am. Statist. Ass., 59, 214–224. Rotsch, C., K. Jacobson, an d M. R admac her (1999). Dimensional and mec hanical dynamics of active and stable edges in motile fibroblasts inves tigated by using atomic force microscop y . Pr o c. Natl. A c ad. Sci. USA, 96, 921–926. Smith, A. F. M. (1975). A Bay esian approac h to inference ab out a change-point in a sequence of random v ariables. Biometrika, 62, 407–416. Smith, A. F. M. (1980). Straigh t lines with a change-point: A Ba yesian analysis of some renal transplant data. Appl. Statist., 29, 180–189. Stephens, D. A. (1994). B ay esian retrosp ective multiple-c hangep oint identification. Appl. Statist., 43, 159–178. W olfe, D. A. and E. Schech tman (1984). N onparametric statistical pro cedures for the change p oin t problem. J. Statist. Plan. Inf er., 9, 386–396. W ong, E. W., P . E. Sheehan, and C. M. Liebert (1997). Nanob eam mechanics: Elasticit y , strength, and toughness of nan oro ds and nanotub es. Scienc e, 277, 1971–1 975. Y u e, S., T. B. M. J. Ouarda, and B. Bob eeb (2001). A review of biv ariate Gamma distributions for hydrologi cal app lication. J. Hydr ol . , 246, 1–18. Y u en, S. G., D. G. Rudoy , R. D. How e, and P . J. W olfe (2007). Ba yesian changepoint d etection for switc h ing regressions in material indentation exp eriments. In Pr o c. 14th IEEE Worksh. Statist. Signal Pr o c ess. , p p. 104–108. Zac ks, S. (1983). Survey of classica l and Bay esian approac hes to the change p oint prob lem: Fixed sample and sequ entia l procedu res of testing and estimation. In R e c ent A dvanc es in Statistics: Herman Chernoff F estschrift , pp. 245–269. N ew Y ork: Academic Press. Zellner, A. (1986). On assessing prior distributions and Bay esian regression analysis with g -p rior distributions. In Bayesian I nfer enc e and De cision T e chniques , pp. 233–24 3. New Y ork: Elsevier.