--- title: "Maximum Approximate Bernstein/Beta Likelihood Estimation in R package 'mable'" author: "Zhong Guan" date: "24/09/2024" output: pdf_document: number_sections: yes toc: yes toc_depth: 3 #keep_tex: true vignette: > %\VignetteIndexEntry{Maximum Approximate Bernstein/Beta Likelihood Estimation in R package 'mable'} %\VignetteEncoding{UTF-8} %\VignetteEngine{knitr::knitr} header-includes: \usepackage{bm} link-citations: yes bibliography: - bernstein.bib biblio-style: chicago --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` \setcounter{section}{0} # Introduction Any continuous density function $f$ on a known closed interval $[a,b]$ can be approximated by Bernstein polynomial $f_m(x; p)=(b-a)^{-1}\sum_{i=0}^m p_i\beta_{mi}[(x-a)/(b-a)]$, where $p_i\ge 0$, $\sum_{i=0}^m p_i=1$ and $\beta_{mi}(u)=(m+1){m\choose i}u^i(1-u)^{m-i}$, $i=0,1,\ldots,m$, is the beta density with shapes $(i+1, m-i+1)$. This provides a way to approximately model the unknwon underlying density with a mixture beta model with an appropriate \emph{model degree} $m$ and solve a nonparametric or semiparametric statistical problem ``parametrically'' using the maximum likelihood method. For instance, based on one-sample data, $x_1,\ldots,x_n$, one can estimate a nonparametric density $f$ parametrically by maximizing the approximate likelihood $\ell(p)=\sum_{j=0}^n\log f_m(x_j; p)$ with respect to $p$ [@Guan-jns-2015]. Since the Bernstein polynomial model of degree $m$ is nested in the model of degree $m+1$, the maximum likelihood is increasing in $m$. The increase is negligible when the model becomes overfitting. Therefore an optimal degree can be chosen as the change-point of the log likelihood ratios over a set of consecutive candidate model degrees. This approach works surprisingly well for even more complicated models and data. With an estimate $\hat p=(\hat p_0,\ldots,\hat p_m)$ of $p$ one can estimate the cumulative distribution function $F$ by $\hat F(x)=F_m(x;\hat p)=\sum_{i=0}^m \hat p_i B_{mi}[(x-a)/(b-a)]$, where $B_{mi}(u)=\int_0^u\beta_{mj}(t)dt$, $i=0,\ldots,m$, is the beta distribution function with shapes $(i+1,m-i+1)$. This manual will illustrate the use of the R package \texttt{mable} for obtaining not only smooth estimates of density, cumulative distribution, and survival functions but also estimates of parameters such as regression coefficients. ```{r echo = FALSE} library(mable) ``` # One-sample Problems ## Raw Data Let $x_1,\ldots,x_n$ be a sample from a population with cumulative distribution function $F$ and density function $f$ on $[a,b]$. If $[a,b]$ is unknown we choose $[a,b]\supset [x_{(1)},x_{(n)}]$, where $x_{(1)}$ and $x_{(n)}$ are the minimum and maximum statistics, respectively. A lower bound for unimodal density is $m_b=\mu(1-\mu)/\sigma^2-3$, where $\mu$ and $\sigma^2$ are, respectively, the mean and variance of the distribution after transformation $Y=(X-a)/(b-a)$. ### Example: Vaal River Annual Flow Data For the annual flow data of Vaal River at Standerton as given by Table 1.1 of @Linhart-and-Zucchini-model-selection-book give the flow in millions of cubic metres, ```{r Vaal-Flow-data} data(Vaal.Flow) head(Vaal.Flow, 3) ``` we want to estimate the density and the distribution functions of annual flow ```{r results = "hide", warning = FALSE, message=FALSE, eval=FALSE} vaal<-mable(Vaal.Flow$Flow, M = c(2,100), interval = c(0, 3000), IC = "all", controls = mable.ctrl(sig.level = 1e-8, maxit = 2000, eps = 1.0e-9)) ``` % the following hided data contain object vaal obtained by dput(vaal) ```{r echo=FALSE} vaal<-structure(list(p = c(6.66395112071704e-142, 0.104927778267183, 0.759431290544939, 2.3090355435799e-08, 2.50844096058296e-16, 7.19449983525583e-21, 1.71511749770294e-20, 1.04929393999386e-07, 0.120235385140723, 1.09761247086527e-22, 3.68974628342609e-96, 2.26150112042313e-228, 0, 0, 0, 5.03890934145742e-234, 1.22687397237239e-39, 0.015405418027405, 4.73545737093756e-197, 0), m = 19L, mloglik = 56.8193876978043, lk = c(41.7728484408995, 46.2771336776274, 47.5575720764407, 47.6039418528322, 48.5523633781241, 50.0345791055467, 51.5963161932995, 52.9697499606394, 54.1010902527862, 54.8308740677684, 55.2629525882091, 55.3303154820947, 55.5239461544182, 55.7498362671488, 56.0447387198335, 56.3500212973192, 56.5685300185606, 56.8193876978043, 56.9055159628123, 57.0383849150017, 57.1452107710574, 57.2967564537662, 57.3728544332059, 57.5168152409349, 57.5743141045582, 57.6944740233863, 57.7495227312013, 57.8363891782184, 57.9032785845485, 58.0240496191399, 58.1061070158426, 58.1552998278603, 58.2422277990116, 58.3429453821828, 58.4797317517046, 58.5424968600091, 58.6506339000825, 58.7235093706809, 58.8483282630303, 58.9139053283805, 58.9608054306796, 59.0001871250332, 59.0609102368121, 59.1276739281344, 59.1742625540794, 59.2091853525299, 59.2297688141381, 59.2731312947835, 59.3234170202736, 59.3593795886725, 59.3884266355113, 59.4426647316214, 59.4916516076017, 59.522517849804, 59.5364175591817, 59.5794196363324, 59.6176201933844, 59.6758686843946, 59.7037694629491, 59.7386829043747, 59.7837475696609, 59.8444248872729, 59.9057645787641, 59.946513968038, 59.9958410065159, 60.0268585687398, 60.0748399114557, 60.1117605895372, 60.1436353540656, 60.1740335482315, 60.2432060831327, 60.3023533782479, 60.3613696685541, 60.4010676861201, 60.4502166815203, 60.4783364018285, 60.5205697911987, 60.5570413921821, 60.5863606010283, 60.6229278768838, 60.6856622848798, 60.7508809200547, 60.7981961100826, 60.8346639232035, 60.8776955928298, 60.9181895275888, 60.9702875905045, 61.0126462333645, 61.0594245354054, 61.1087161234944, 61.1665319394629, 61.2379063045393, 61.2911041469829, 61.3367859527923, 61.3733789317552, 61.4027683168247, 61.4328921127991, 61.4629513058874, 61.499450177195 ), lr = c(21.0400267707794, 26.0092133279022, 23.5318492590107, 27.1434666969398, 35.0730762259421, 45.0102014514883, 55.0427504607432, 64.3160676551957, 70.1759232548958, 72.8411351282936, 70.8390029810927, 70.6058943957985, 70.913116777793, 72.3210882063501, 74.0712067197948, 74.6759893592703, 75.9233214868143, 74.6125325680211, 74.11480170489, 73.2428708875183, 73.1731617162639, 71.8783447150091, 71.7897434521724, 70.2500262869158, 69.8389254768924, 68.3221259221015, 67.3977774586726, 66.1517508921878, 65.9001754860353, 64.9889379933534, 63.5082861192643, 62.7409911064193, 62.2599432429687, 62.5024555203726, 61.3737106030307, 61.1511118160824, 60.2716006916461, 60.4655773598734, 59.4958453862877, 58.1640101620904, 56.6982670143826, 55.688825543864, 54.8243992024348, 53.562035565918, 52.0769299250977, 50.319120648889, 49.0440768198218, 47.9255334501533, 46.5338116404446, 45.0213744746735, 44.0275676832284, 42.9438737512301, 41.5143117622436, 39.7718738672515, 38.6141051565547, 37.3778072126404, 36.5402780374695, 35.1330173615615, 33.8755938452882, 32.8231186752621, 32.0773625475224, 31.3581065782058, 30.2572295681155, 29.3315255554462, 28.0691814610304, 27.1394226728944, 26.012535363232, 24.8046264740693, 23.5842670932705, 23.0832108801469, 22.4117764825536, 21.7504220040654, 20.7340621741221, 19.9056013250303, 18.6922002503777, 17.7528450976514, 16.7174752298143, 15.5655406148263, 14.5559949910125, 14.0167055588444, 13.5382702672744, 12.7432972366204, 11.7568937867725, 10.899274118421, 10.0039350267711, 9.32550672065937, 8.47840255713019, 7.71875131060652, 7.01266859001961, 6.4780592938132, 6.25920963599765, 5.70827336281335, 5.01235671765166, 4.10225359885751, 3.00896962305403, 1.94465924213493, 0.892233157440845, 0), M = c(2L, 100L), interval = c(0, 3000), pval = c(1, 1, 1, 0.135747987070765, 0.270328157892676, 0.349180754949455, 0.417371289620143, 0.460063253179317, 0.479448013538027, 0.469838893452404, 0.441298754319351, 0.349364191505256, 0.249319863136916, 0.175688495338965, 0.139862041470063, 0.114616780419735, 0.0855621111525051, 0.0679560880517338, 0.0452015731977289, 0.0324854801085018, 0.0230538559228871, 0.0175784438393906, 0.012427418152239, 0.00961177306691063, 0.00669187842802754, 0.00515277872079156, 0.00367314466606472, 0.00276885425454754, 0.00205277822191852, 0.00165612736486742, 0.00127881266634433, 0.000952575127481037, 0.000753455575955519, 0.000611934040600226, 0.000525146786109709, 0.000411029686438358, 0.000343852178177007, 0.000276276502896855, 0.000238893196351286, 0.000192329997576457, 0.000151963830109736, 0.000119503958089573, 9.70962248756368e-05, 7.98883244934601e-05, 6.43153503118166e-05, 5.12122767023504e-05, 4.01749532250584e-05, 3.26037760055575e-05, 2.68062881254583e-05, 2.17041778300953e-05, 1.74722970925911e-05, 1.45979096953797e-05, 1.21515894859758e-05, 9.90397528388698e-06, 7.91263008670384e-06, 6.59506291256218e-06, 5.47716145693489e-06, 4.68789467344966e-06, 3.8606809591446e-06, 3.21841458639227e-06, 2.72758066366396e-06, 2.36745136505956e-06, 2.052434471711e-06, 1.70484313655184e-06, 1.44424936665555e-06, 1.18425169826075e-06, 1.00636820077327e-06, 8.3940975292851e-07, 6.95241539672153e-07, 5.75707648220458e-07, 5.1494484343273e-07, 4.52737899725442e-07, 3.98772506349232e-07, 3.3940110522046e-07, 2.94708836445778e-07, 2.46536952075438e-07, 2.12232688534542e-07, 1.81119941378149e-07, 1.52830039401586e-07, 1.30983814239372e-07, 1.18092237877399e-07, 1.07116784331396e-07, 9.41610780458291e-08, 8.1287126474372e-08, 7.1146715407977e-08, 6.20893869651695e-08, 5.54313032141707e-08, 4.86971299951122e-08, 4.31916912235764e-08, 3.85396806690252e-08, 3.4965970452383e-08, 3.2534852589805e-08, 2.9347718255579e-08, 2.61600523465688e-08, 2.29820830144334e-08, 1.99653764632046e-08, 1.7391248330334e-08, 1.51680690230194e-08, 1.33969156879132e-08), ic = list(BIC = c(37.5984611710038, 40.0155527727839, 43.383184806545, 41.3423609479887, 40.2035888383328, 43.7729982007032, 41.1603480185604, 46.7081690557959, 47.8395093479428, 48.569293162925, 46.9141780484178, 44.8943473073557, 47.1751716146269, 47.4010617273576, 47.6959641800422, 43.8268594876323, 44.0453682088737, 44.2962258881174, 44.3823541531254, 40.3408358354192, 42.5348553264226, 42.6864010091314, 44.849692623519, 44.993653431248, 42.9639586599235, 43.0841185787516, 43.1391672865666, 43.2260337335837, 41.205729504966, 41.3265005395573, 41.40855793626, 41.4577507482778, 41.5446787194291, 37.4710090327046, 39.6949890371742, 41.8449477804266, 35.6915039156565, 37.8515730212027, 37.9763919135522, 38.0419689789023, 38.0888690812015, 36.0410571406072, 38.188973887334, 36.1685439437084, 36.2151325696534, 36.2500553681039, 34.1834451947643, 32.1396140404619, 34.2770934008998, 28.0514750644552, 34.3421030161375, 26.0475665724564, 28.1837470833844, 24.0402260556911, 24.0541257650688, 22.0099342072717, 24.1353283992715, 22.1063832553339, 20.0470903989406, 22.169197475314, 20.1270685056524, 24.36213309316, 22.3362791497034, 26.5514158088729, 22.4263555774552, 24.5445667746269, 22.505354482395, 26.7166624303722, 22.5741499250049, 22.6045481191708, 24.7609142890198, 24.820061584135, 24.8790778744412, 27.005969526955, 24.9679248874074, 24.9960446077156, 22.951084362138, 22.9875559631213, 25.1040688069154, 23.0534424478231, 23.116176855819, 27.3557827608896, 25.3159043159697, 23.2651784941427, 25.3954037987169, 25.4358977334759, 25.4879957963916, 27.6175480741994, 27.6643263762403, 27.7136179643293, 25.68424014535, 25.7556145104264, 27.8960059878178, 25.8544941586794, 23.8038935026945, 23.833282887764, 21.7762130487905, 19.719078606931, 21.8427711131865), AIC = c(39.7728484408995, 43.2771336776273, 45.5575720764407, 44.6039418528322, 44.5523633781241, 47.0345791055467, 46.5963161932995, 49.9697499606394, 51.1010902527862, 51.8308740677684, 51.2629525882091, 50.3303154820948, 51.5239461544182, 51.7498362671488, 52.0447387198335, 50.3500212973192, 50.5685300185606, 50.8193876978043, 50.9055159628123, 49.0383849150017, 50.1452107710574, 50.2967564537662, 51.3728544332059, 51.5168152409349, 50.5743141045582, 50.6944740233863, 50.7495227312013, 50.8363891782184, 49.9032785845485, 50.0240496191399, 50.1061070158426, 50.1552998278603, 50.2422277990116, 48.3429453821828, 49.4797317517046, 50.5424968600091, 47.6506339000825, 48.7235093706809, 48.8483282630303, 48.9139053283805, 48.9608054306796, 48.0001871250332, 49.0609102368121, 48.1276739281344, 48.1742625540794, 48.2091853525299, 47.2297688141381, 46.2731312947835, 47.3234170202736, 44.3593795886725, 47.3884266355113, 43.4426647316214, 44.4916516076017, 42.522517849804, 42.5364175591817, 41.5794196363324, 42.6176201933844, 41.6758686843946, 40.7037694629491, 41.7386829043747, 40.7837475696609, 42.8444248872729, 41.9057645787641, 43.946513968038, 41.9958410065159, 43.0268585687398, 42.0748399114557, 44.1117605895372, 42.1436353540656, 42.1740335482315, 43.2432060831327, 43.3023533782479, 43.3613696685541, 44.4010676861201, 43.4502166815203, 43.4783364018285, 42.5205697911987, 42.5570413921821, 43.5863606010283, 42.6229278768838, 42.6856622848798, 44.7508809200547, 43.7981961100826, 42.8346639232035, 43.8776955928298, 43.9181895275888, 43.9702875905046, 45.0126462333645, 45.0594245354054, 45.1087161234944, 44.1665319394629, 44.2379063045393, 45.2911041469829, 44.3367859527923, 43.3733789317552, 43.4027683168247, 42.4328921127991, 41.4629513058874, 42.499450177195), QHC = c(38.9149132692382, 41.9902309201355, 44.6996369047794, 43.3170390953403, 42.8364930348016, 45.7476763480548, 44.4514782641464, 48.6828472031475, 49.8141874952943, 50.5439713102766, 49.5470822448866, 48.1854775529416, 49.8080758110957, 50.0339659238263, 50.328868376511, 47.7762157823354, 47.9947245035769, 48.2455821828206, 48.3317104478285, 45.6066442283567, 47.142437670243, 47.2939833529518, 48.7990489182221, 48.9430097259512, 47.5715410037439, 47.691700922572, 47.746749630387, 47.833616077404, 46.4715378979035, 46.5923089324949, 46.6743663291976, 46.7235591412153, 46.8104871123667, 44.0532695238766, 45.619023479229, 47.1107561733641, 42.9319904559457, 44.4338335123747, 44.5586524047241, 44.6242294700743, 44.6711295723734, 43.2815436808963, 44.7712343785059, 43.4090304839976, 43.4556191099425, 43.490541908393, 42.0821577841707, 40.6965526789854, 42.1758059903062, 37.9248658012131, 42.2408156055438, 36.5791833583315, 38.0571378201423, 35.2300688906834, 35.2439686000612, 33.8580030913812, 35.3251712342638, 33.9544521394434, 32.5533853321673, 34.0172663594235, 32.6333634388791, 35.5519759281523, 34.1843480338129, 37.0830325947481, 34.2744244615647, 35.7344096096192, 34.3534233665045, 37.2482792162473, 34.4222188091144, 34.4526170032803, 35.9507571240121, 36.0099044191273, 36.0689207094335, 37.5375863128302, 36.1577677223998, 36.1858874427079, 34.7991532462475, 34.8356248472309, 36.2939116419077, 34.9015113319326, 34.9642457399286, 37.8873995467647, 36.505747150962, 35.1132473782522, 36.5852466337092, 36.6257405684682, 36.677838631384, 38.1491648600745, 38.1959431621154, 38.2452347502044, 36.8740829803423, 36.9454573454187, 38.4276227736929, 37.0443369936717, 35.651962386804, 35.6813517718735, 34.2825079820172, 32.8835995892749, 34.3490660464132)), chpts = c(2L, 2L, 2L, 4L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 12L, 12L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L), convergence = 2L, delta = 2.65321666859109e-07, xNames = "Vaal.Flow$Flow", data.type = "raw"), class = "mable") ``` Here we truncate the density by \texttt{interval=c(0, 3000)} and choose an optimal degree $m$ among the candidate degrees M[1]:M[2] using the method of change-point. The maximum number of iterations is \texttt{maxit} and the convergence criterion is \texttt{eps} for each $m$ of M[1]:M[2]. The search of an optimal degree stops when the $p$-value \texttt{pval} of change-point test falls below the specified significance level \texttt{sig.level} or the largest degree, M[2], has been reached. If the latter occurs a warning message shows up. In this case we should check the last value of \texttt{pval}. In the above example, we got warning message and the last \texttt{pval}, `r vaal$pval[vaal$M[2]-vaal$M[1]+1]`, which is small enough. The selected optimal degree is $m=$ `r vaal$m[1]`. One can also look at the Bayesian information criteria, \texttt{BIC}, and other information criteria, Akaike information criterion \texttt{AIC} and Hannan–Quinn information criterion \texttt{QHC}, at each candidate degree. These information criteria are not reliable due to the difficulty of determining the model dimension. The \texttt{plot} method for \texttt{mable} class object can visualize some of the results. ```{r vaal-river-data-plot, fig.align='center', fig.cap="Vaal River Annual Flow Data", fig.width=7, fig.height=7} op <- par(mfrow = c(1,2)) layout(rbind(c(1, 2), c(3, 3))) plot(vaal, which = "likelihood", cex = .5) plot(vaal, which = "change-point", lgd.x = "topright") hist(Vaal.Flow$Flow, prob = TRUE, xlim = c(0,3000), ylim =c(0,.0022), breaks =100*(0:30), main = "Histogram and Densities of the Annual Flow of Vaal River", border = "dark grey",lwd = 1, xlab = "Flow", ylab = "Density", col ="light grey") lines(density(x<-Vaal.Flow$Flow, bw = "nrd0", adjust = 1), lty = 2, col = 2,lwd = 2) lines(y<-seq(0, 3000, length=100), dlnorm(y, mean(log(x)), sqrt(var(log(x)))), lty = 4, col = 4, lwd = 2) plot(vaal, which = "density", add = TRUE, lwd = 2) legend("topright", lty = c(1, 4, 2), col = c(1, 4, 2), bty = "n",lwd = 2, c(expression(paste("MABLE: ",hat(f)[B])), expression(paste("Log-Normal: ",hat(f)[P])), expression(paste("KDE: ",hat(f)[K])))) par(op) ``` In Figure \ref{fig:vaal-river-data-plot}, the unknown density $f$ is estimated by \texttt{MABLE} $\hat f_B$ using optimal degrees $m=$ `r vaal$m` selected using the exponential change-point method, the parametric estimate using \texttt{Log-Normal} model and the kernel density estimate \texttt{KDE}: $\hat f_K$. We can also look at the plots of AIC, BIC, and QHC, and likelihood ration(LR) of gamma change-point. ```{r Vaal-River-data-AIC-BIC-plot, fig.align='center', fig.cap="AIC and BIC based on Vaal River Data", fig.width=6, fig.height=4} M <- vaal$M[1]:vaal$M[2] aic <- vaal$ic$AIC bic <- vaal$ic$BIC qhc <- vaal$ic$QHC vaal.gcp <- optim.gcp(vaal) # choose m by gamma change-point model lr <- vaal.gcp$lr plot(M, aic, cex = 0.7, xlab = "m", ylab = "", main = "AIC, BIC, QHC, and LR", ylim = c(ymin<-min(aic,bic,qhc,lr), ymax<-max(aic,bic,qhc,lr)), col = 1) points(M, bic, pch = 2, cex = 0.7, col = 2) points(M, qhc, pch = 3, cex = 0.7, col = 3) points(M[-1], lr, pch = 4, cex = 0.7, col = 4) segments(which.max(aic)+M[1]-1->m1, ymin, m1, max(aic), lty = 2) segments(which.max(bic)+M[1]-1->m2, ymin, m2, max(bic), lty = 2, col = 2) segments(which.max(qhc)+M[1]-1->m3, ymin, m3, max(qhc), lty = 2, col = 3) segments(which.max(lr)+M[1]->m4, ymin, m4, max(lr), lty = 2, col = 4) axis(1, c(m1,m2, m3, m4), as.character(c(m1,m2,m3,m4)), col.axis = 4) legend("topright", pch=c(1,2,3,4), c("AIC", "BIC", "QHC", "LR"), bty="n", col=c(1,2,3,4)) ``` From Figure \ref{fig:Vaal-River-data-AIC-BIC-plot} we see that the gamma change-point method gives the same optimal degree $m=$ `r m4` as the exponential method; BIC and QHC give the same degree $m=$ `r m2`; the degree $m=$ `r m1` selected by AIC method is closer to the one selected by change-point methods. From this Figure we also see that unlike the LR plot the information criteria do not have clear peak points. For any given degree $m$, one can fit the data by specifying \texttt{M=m} in \texttt{mable()} to obtain an estimate of $f$. The \texttt{summary} method \texttt{summary.mable} prints and returns invisibly the summarized results. ```{r} summary(vaal) ``` The mixing proportions, the coefficients of the Bernstein polynomials, $p$ can be obtained as either \texttt{p<-vaal\$p} ```{r} p<-vaal$p p ``` or \texttt{summary(res)\$p}. The method of change-point for choosing model degree is computer-intensive especially for multimodal density. A better lower bound for model degree is $$m_b=\frac{\mu(1-\mu)-\sigma^2}{\sigma^2-\sum_{i=1}^k\lambda_i(\mu_i-\mu)^2}-2,$$ where $\lambda_i$ and $\mu_i$ are, respectively, the mixing proportion and mean value of teh $i$th component density. Less computer-intensive methods such as the method of moment and the method of mode implemented by \texttt{optimal()} can be used (see Figure \ref{fig:old-faithful-amrginal-data-plot}). ### Example: The Old Faithful data \label{old-faithful-data-mme} ```{r old-faithful-data-mme, echo=TRUE, message=FALSE, warning=FALSE, results = "hide"} x<-faithful x1<-faithful[,1] x2<-faithful[,2] a<-c(0, 40); b<-c(7, 110) mu<-(apply(x,2,mean)-a)/(b-a) s2<-apply(x,2,var)/(b-a)^2 # mixing proportions lambda<-c(mean(x1<3),mean(x2<65)) # guess component mean mu1<-(c(mean(x1[x1<3]), mean(x2[x2<65]))-a)/(b-a) mu2<-(c(mean(x1[x1>=3]), mean(x2[x2>=65]))-a)/(b-a) # estimate lower bound for m mb<-ceiling((mu*(1-mu)-s2)/(s2-lambda*(1-lambda)*(mu1-mu2)^2)-2) ``` ```{r old-faithful-data-optimable, echo=TRUE, message=FALSE, warning=FALSE, eval=FALSE, results = "hide"} m1<-optimable(x1, interval=c(a[1],b[1]), nmod=2, modes=c(2,4.5))$m m2<-optimable(x2, interval=c(a[2],b[2]), nmod=2, modes=c(52.5,80))$m erupt1<-mable(x1, M=mb[1], interval=c(a[1],b[1])) erupt2<-mable(x1, M=m1, interval=c(a[1],b[1])) wait1<-mable(x2, M=mb[2],interval=c(a[2],b[2])) wait2<-mable(x2, M=m2,interval=c(a[2],b[2])) ``` ```{r echo=FALSE} mb<-c(78,28) m1<-122 m2<-34 erupt1<-structure(list(p = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 1.24815896595563e-216, 3.38239376086288e-68, 0.257770688882492, 0.0979264711500575, 2.31634344879214e-53, 2.66253393626854e-138, 1.48544938482616e-237, 1.48219693752374e-323, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 1.48219693752374e-323, 1.72038175246372e-275, 1.94836109530888e-208, 1.79962617110949e-146, 1.30190873388111e-93, 2.01576739697645e-52, 5.69539554300385e-24, 8.58124703698123e-08, 0.0578639985565491, 0.00216726804817309, 2.45410656052947e-08, 3.8391464010493e-13, 5.14426500405859e-15, 9.58951405603689e-13, 3.35442886601859e-07, 0.154674097429257, 0.429597030135244, 4.57523201601482e-13, 4.80861465801093e-46, 1.60317842850254e-109, 1.85518266050181e-213, 9.88131291682493e-324, 9.88131291682493e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), m = c(eruptions = 78), mloglik = 249.26925075986, interval = c(0, 7), convergent = TRUE, del = 1e-07, xNames = "x1", data.type = "raw"), class = "mable") erupt2<-structure(list(p = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 4.30737643561733e-269, 1.62814062123862e-131, 4.90911843705386e-46, 2.83917053657873e-05, 0.348368171335392, 3.77694875446027e-25, 3.71830537130198e-68, 4.22179465798189e-121, 3.21584252902755e-175, 1.20894654106119e-222, 4.72975981679665e-257, 2.74709261002431e-274, 8.18149514208777e-273, 8.44589795695613e-254, 9.08617757941952e-221, 7.77296368749942e-179, 9.65734936105495e-134, 7.0036228159251e-91, 1.53780825637889e-54, 2.05562181625119e-27, 1.92411289705749e-10, 0.00638131265972255, 0.00741764733792521, 1.73042840569444e-07, 6.50478488180998e-14, 1.15942674061074e-19, 4.0249794950853e-23, 1.07928705660775e-23, 1.11787566542467e-21, 5.88049078257417e-18, 1.34318951955047e-13, 1.61134983988649e-09, 2.7405485571206e-06, 0.000433646477745988, 0.00826425048957231, 0.0339435343542211, 0.0527620183123396, 0.0434143468745409, 0.0205586948367354, 0.00540714866412622, 0.000831386742608205, 0.000101910891977919, 1.84126954946833e-05, 1.07457809214495e-05, 3.84061367255896e-05, 0.000881805715994463, 0.0460173726014059, 0.417170226473733, 0.00794765398033802, 5.37755739563085e-10, 3.26963511653428e-26, 8.53102373923467e-56, 1.20693339901693e-103, 4.53327375461247e-175, 1.2439128875044e-275, 9.88131291682493e-324, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), m = 122, mloglik = 259.459509914242, interval = c(0, 7), convergent = FALSE, del = 3.10216285015485e-06, xNames = "x1", data.type = "raw"), class = "mable") wait1<-structure(list(p = c(0, 9.88131291682493e-324, 6.6849130023567e-74, 0.0187548016948242, 0.15700514789747, 8.75954273815904e-05, 0.06545445419232, 0.119002640079203, 7.05858130071315e-13, 5.78009594099882e-38, 2.4885050618525e-70, 4.39810503775073e-96, 1.23641913319538e-101, 2.75300216831121e-83, 2.61907632865599e-50, 1.43811814261486e-18, 0.563846543514944, 0.0758488171931512, 7.39718007440611e-18, 1.76520675346944e-45, 2.17723696161449e-92, 2.91762259483708e-186, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0), m = c(waiting = 28), mloglik = 123.361607234209, interval = c(40, 110), convergent = TRUE, del = 1e-07, xNames = "x2", data.type = "raw"), class = "mable") wait2<-structure(list(p = c(0, 4.94065645841247e-324, 1.14264859149383e-232, 1.30720378382645e-40, 0.0838405408737573, 0.083619306273632, 0.00136187033955907, 0.00677135346637448, 0.100407492563723, 0.0895524269017653, 1.81869219281019e-06, 7.72322130777564e-20, 7.52703369686382e-42, 8.49919563980898e-66, 3.79021098912416e-81, 4.78224136880441e-80, 1.47865073410181e-62, 1.12999728173389e-36, 1.63428527813762e-13, 0.283502336517345, 0.334317016669562, 3.78931951581921e-07, 1.30305547325536e-09, 0.000245257916097455, 0.016380199550821, 1.18598828606837e-36, 1.59695089211477e-158, 9.88131291682493e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0), m = 34, mloglik = 124.360653955066, interval = c(40, 110), convergent = FALSE, del = 4.86893426909774e-07, xNames = "x2", data.type = "raw"), class = "mable") ``` For this data set, we obtain lower bounds for the marginal densities $m_b=(`r mb[1]`, `r mb[2]`)$ and $\hat m=(`r m1`, `r m2`)$. ```{r old-faithful-amrginal-data-plot, echo=TRUE, message=FALSE, warning=FALSE, results = "hide", fig.align='center', fig.cap="The Old Faithful data", fig.width=6, fig.height=4} op<-par(mfrow=c(1,2), cex=0.8) hist(x1, probability = TRUE, col="grey", border="white", main="", xlab="Eruptions", ylim=c(0,.65), las=1) plot(erupt1, add=TRUE,"density") plot(erupt2, add=TRUE,"density",lty=2,col=2) legend("topleft", lty=c(1,2),col=1:2, bty="n", cex=.7, c(expression(paste("m = ", m[b])),expression(paste("m = ", hat(m))))) hist(x2, probability = TRUE, col="grey", border="white", main="", xlab="Waiting", las=1) plot(wait1, add=TRUE,"density") plot(wait2, add=TRUE,"density",lty=2,col=2) legend("topleft", lty=c(1,2),col=1:2, bty="n", cex=.7, c(expression(paste("m = ", m[b])),expression(paste("m = ", hat(m))))) par(op) ``` ## Grouped Data With a grouped dataset, a frequency table of data from a continuous population, one can estimate the density from histogram using \texttt{mable.group()} with an optimal degree $m$ chosen from \texttt{M[1]:M[2]} or with a given degree $m$ using \texttt{M=m} [@Guan-2017-jns]. ### Example: The Chicken Embryo Data Consider the chicken embryo data contain the number of hatched eggs on each day during the 21 days of incubation period. The times of hatching ($n=43$) are treated as grouped by intervals with equal width of one day. The data were studied first by @Jassim-et-al-1996. @Kuurman-et-al-2003 and @Lin-and-He-2006-bka also analyzed the data using the minimum Hellinger distance estimator, in addition to other methods assuming some parametric mixture models including Weibull model. ```{r} data(chicken.embryo) head(chicken.embryo, 2) a <- 0 b <- 21 day <- chicken.embryo$day nT <- chicken.embryo$nT ``` ```{r results = "hide", warning = FALSE, eval=FALSE} embryo<-mable.group(x = nT, breaks = a:b, M=c(2,100), interval = c(a, b), IC = "aic", controls = mable.ctrl(sig.level = 1e-6, maxit = 2000, eps = 1.0e-7)) ``` ```{r echo=FALSE} embryo<-structure(list(p = c(0.114376931420402, 0.562791099356775, 5.32839664631301e-19, 3.38622777062756e-51, 2.84638292790797e-62, 5.80521798084961e-42, 1.3375757219419e-10, 0.0221221274888444, 2.32016100251404e-19, 5.82119897096528e-33, 5.79803303109582e-25, 0.0146874328491674, 0.286022408751054, 6.42237557241518e-50), mloglik = -107.317973971343, m = 13L, lk = c(-120.863259740031, -115.817965358717, -112.801354117708, -111.339648775425, -110.837265980592, -110.538493633193, -110.020689312112, -109.298178441611, -108.600045047698, -108.02557800418, -107.585638457218, -107.317973971343, -107.148545461053, -107.055966733679, -106.962452125559, -106.827342660057, -106.675514124012, -106.507298389744, -106.382851894322, -106.343553530154, -106.32440630411, -106.276574667465, -106.234080500234, -106.173604014522, -106.117202622847, -106.051240207298, -106.04113558772, -106.01265002765, -105.993737313517, -105.952654698097, -105.923203482832, -105.891479080936, -105.863368579659, -105.813513347504, -105.77608117454, -105.70208187719, -105.63713301691, -105.574326688937, -105.535462089958, -105.478437815024, -105.421349867517, -105.369197146256, -105.329451739024, -105.249852983251, -105.179117104329, -105.108501370057, -104.995452636252, -104.87999862755, -104.792658763962, -104.735615142702, -104.70061898287, -104.642623541647, -104.57707993091, -104.514966078791, -104.421327923093, -104.323686475727, -104.233815881355, -104.175575968485), lr = c(16.3230095456385, 29.0861065901675, 35.5917280060074, 36.2244662618225, 35.5961320336328, 36.9092000352466, 40.3857805702116, 44.3938517294349, 47.8959924401617, 50.4479709734884, 51.1408843739782, 50.6491244341255, 49.1547046764929, 47.7418558334538, 47.0014095938608, 46.5952326351995, 46.5470277981241, 45.8989731404594, 43.9161078582619, 41.6622622072695, 39.9148221305205, 38.1267477312251, 36.657180483749, 35.1605842594247, 33.8449803360561, 31.7168974676527, 29.9135101092751, 28.0178986716924, 26.4763347569636, 24.8114629046842, 23.2168258039656, 21.6148538745695, 20.3229395307165, 18.9068623738782, 17.9610145023234, 16.9259642996796, 15.8855583728488, 14.586933447357, 13.5239029400836, 12.4828494629328, 11.4095964315732, 10.2315111028659, 9.4880709572956, 8.66844603009012, 7.86276435685585, 7.50384102297772, 7.20640364067839, 6.62416451871473, 5.70882918770703, 4.57964608639039, 3.72343300913795, 2.96039616405586, 2.18684390983753, 1.70251083343507, 1.27496068347981, 0.819409224233579, 0), M = c(2L, 59L), interval = c(0, 21), convergence = 1L, del = 8.20010057264204e-06, ic = list(BIC = c(-124.624459855725, -121.459765532257, -118.443154291249, -115.100848891119, -116.479066154132, -116.180293806734, -119.423689601346, -118.701178730845, -116.122445279085, -117.428578293413, -116.988638746452, -116.720974260577, -118.432145808134, -118.339567080759, -118.24605247264, -119.991543064985, -117.959114471092, -117.790898736824, -119.547052299249, -125.149554108622, -121.369206766885, -119.440775072393, -123.159481020855, -123.099004535143, -126.803803259161, -124.857240785766, -126.727736224035, -128.579850721812, -128.560938007678, -124.758655276565, -128.490404176993, -128.458679775098, -130.311169331667, -130.261314099513, -130.223881926548, -135.791682802739, -137.607334000306, -135.663927614486, -133.74446295766, -131.806838624879, -129.869150619525, -137.339398129652, -129.777252491032, -125.936453619565, -125.865717740644, -125.795102006372, -127.562653330413, -138.730799668792, -131.121059573817, -125.422215779017, -121.626019503492, -125.329224177962, -125.263680567225, -125.201566715106, -125.107928559407, -128.771487227735, -123.039816459823, -126.742776662646), AIC = c(-122.863259740031, -118.817965358717, -115.801354117708, -113.339648775425, -113.837265980592, -113.538493633193, -115.020689312112, -114.298178441611, -112.600045047698, -113.02557800418, -112.585638457218, -112.317973971343, -113.148545461053, -113.055966733679, -112.962452125559, -113.827342660057, -112.675514124012, -112.507298389744, -113.382851894322, -116.343553530154, -114.32440630411, -113.276574667465, -115.234080500234, -115.173604014522, -117.117202622847, -116.051240207298, -117.04113558772, -118.01265002765, -117.993737313517, -115.952654698097, -117.923203482832, -117.891479080936, -118.863368579659, -118.813513347504, -118.77608117454, -121.70208187719, -122.63713301691, -121.574326688937, -120.535462089958, -119.478437815024, -118.421349867517, -122.369197146256, -118.329451739024, -116.249852983251, -116.179117104329, -116.108501370057, -116.995452636252, -122.87999862755, -118.792658763962, -115.735615142702, -113.70061898287, -115.642623541647, -115.57707993091, -115.514966078791, -115.421327923093, -117.323686475727, -114.233815881355, -116.175575968485)), pval = c(1, 1, 1, 0.258321430724618, 0.338416305131659, 0.250214806435255, 0.214437581083489, 0.20868755191322, 0.197315929355982, 0.175642869925156, 0.14821286270032, 0.116640782330745, 0.0884706406042653, 0.0651144119898133, 0.0483829749456973, 0.0371837728546789, 0.0291722924685416, 0.0233577113779639, 0.0184075792524701, 0.0138979785335623, 0.00869865992416541, 0.00557764814071127, 0.00361092339215108, 0.00247156882260313, 0.00170644473233761, 0.0012203212088614, 0.000770968855261067, 0.000514040664339466, 0.000339161727812232, 0.000239344568407684, 0.000166006143030195, 0.000117049776315992, 8.26142912390138e-05, 6.21789404039452e-05, 4.57675691831749e-05, 3.71939755646755e-05, 2.97547025729372e-05, 2.38360917681479e-05, 1.81563941096252e-05, 1.45346636180044e-05, 1.17058388606761e-05, 9.37289203561953e-06, 7.333857852454e-06, 6.32000816869205e-06, 5.35795229139602e-06, 4.5588365490401e-06, 4.26544467413414e-06, 4.00992976457015e-06, 3.56012781232984e-06, 2.98124930497856e-06, 2.39718559091884e-06, 2.02683100458678e-06, 1.74562533694633e-06, 1.49792137804639e-06, 1.37074402284387e-06, 1.26568991221099e-06, 1.15315511173275e-06, 9.92595902138405e-07 ), chpts = c(2L, 2L, 2L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 12L, 12L, 12L, 12L, 12L, 12L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L), xNames = "nT", data.type = "grp"), class = "mable") ``` ```{r chicken-embryo-data-plot, fig.align='center', fig.cap="Chicken Embryo Data", fig.width=6, fig.height=7} Day <- rep(day,nT) op <- par(mfrow = c(1,2), lwd = 2) layout(rbind(c(1, 2), c(3, 3))) plot(embryo, which = "likelihood") plot(embryo, which = "change-point") fk <- density(x = rep((0:20)+.5, nT), bw = "sj", n = 101, from = a, to = b) hist(Day, breaks = seq(a,b, length = 12), freq = FALSE, col = "grey", border = "white", main = "Histogram and Density Estimates") plot(embryo, which = "density", cex = 0.7, add = TRUE) lines(fk, lty = 2, col = 2) legend("top", lty = c(1:2), c("MABLE", "Kernel"), bty = "n", col = c(1:2)) par(op) ``` We see in Figure \ref{fig:chicken-embryo-data-AIC-BIC-plot} that AIC and gamma change-point method give the same optimal degree as the one, $m=$ `r embryo$m`, given by the exponential change-point method. However, BIC fails in choosing a useful model degree. ```{r chicken-embryo-data-AIC-BIC-plot, fig.align='center', fig.cap="AIC and BIC based on Chicken Embryo Data", fig.width=6, fig.height=4} M <- embryo$M[1]:embryo$M[2] aic <- embryo$ic$AIC bic <- embryo$ic$BIC res.gcp <- optim.gcp(embryo) # choose m by gamma change-point model lr <- res.gcp$lr plot(M, aic, cex = 0.7, col = 1, xlab = "m", ylab = "", ylim = c(ymin<-min(aic,bic,lr), ymax<-max(aic,bic,lr)), main = "AIC, BIC, and LR") points(M, bic, pch = 2, cex = 0.7, col = 2) points(M[-1], lr, pch = 3, cex = 0.7, col = 4) segments(which.max(aic)+M[1]-1->m1, ymin, m1, max(aic), lty = 2) segments(which.max(bic)+M[1]-1->m2, ymin, m2, max(bic), lty = 2, col = 2) segments(which.max(lr)+M[1]->m3, ymin, m3, max(lr), lty = 2, col = 4) axis(1, c(m1,m2, m3), as.character(c(m1,m2,m3)), col.axis = 4) legend("right", pch = 1:3, c("AIC", "BIC", "LR"), bty = "n", col = c(1,2,4)) ``` The results are summarized as follows. ```{r} summary(embryo) ``` ## Contaminated Data--Density Deconvolution Consider the additive measurement error model $Y = X + \epsilon$, where $X$ has an unknown distribution $F$, $\epsilon$ has a known distribution $G$, and $X$ and $\epsilon$ are independent. We want to estimate density $f=F'$ based on independent observations, $y_i = x_i + \epsilon_i$, $i=1,\ldots,n$, of $Y$, where $x_i$'s and $\epsilon_i$'s are unobservable. \texttt{mable.decon()} implements the method of @Guan-2019-mable-deconvolution and gives an estimate of the density $f$ using the approximate Bernstein polynomial model. ### Example: A Simulated Normal Dataset ```{r results = "hide", warning = FALSE} set.seed(123) mu <- 1; sig <- 2; a <- mu - sig*5; b <- mu + sig*5; gn <- function(x) dnorm(x, 0, 1) n <- 50; x <- rnorm(n, mu, sig); e <- rnorm(n); y <- x + e; ``` ```{r results = "hide", warning = FALSE, eval=FALSE} decn <- mable.decon(y, gn, interval = c(a,b), M = c(5, 50)) ``` ```{r echo=FALSE} decn<-structure(list(lk = c(26.4579423107072, 31.4140509648368, 31.6998197144108, 34.9370176142077, 35.2218689572872, 37.4273654796253, 37.7114555490011, 39.2389160569194, 39.5223722599646, 40.576827290447, 40.859759218697, 41.5717041718022, 41.8525657731218, 42.330823537089, 42.5867219176599, 42.9072301430247, 43.1220744872046, 43.3339784277991, 43.5017673605904, 43.6379069244143, 43.757776779825, 43.8406765363159, 43.9143423876334, 43.9597771358505, 43.9901711510228, 44.0132941175817), lr = c(5.02746478365044, 3.60541223095752, 7.52343356261833, 6.31033456027949, 9.44964523254241, 8.35690371851028, 10.8356879380868, 9.86810277277983, 11.7197952837592, 10.907275632404, 12.116971203553, 11.5012874953979, 12.1421805957699, 11.6492334826244, 11.796369369754, 11.322857648053, 11.0153806264742, 10.458442517396, 9.70462338391716, 8.93229705718044, 7.70390485228706, 6.4866761765979, 4.68721784830983, 2.46531643620066, 0), p = c(0, 0, 0, 0, 3.4620755997222e-241, 3.93510147238985e-145, 9.81419110226799e-79, 8.98907407872981e-36, 3.71844957194364e-11, 0.909202372386772, 0.0907976275756359, 4.07903846905113e-13, 5.27217693724487e-36, 4.73536199157314e-74, 2.34291815334911e-133, 6.88953241741248e-223, 0, 0, 0), m = 18L, pval = c(1, 1, 1, 0.248430944519551, 0.28874165936426, 0.471250353765544, 0.426508202618338, 0.488584206014697, 0.448941261229855, 0.455207804767377, 0.396335228868547, 0.37409521001946, 0.312571031616345, 0.272552270565568, 0.218467208799487, 0.181621363026878, 0.142831162959649, 0.110117441998513, 0.0835422496832733, 0.0627933371332644, 0.0458559664816752, 0.0322750310721527, 0.0228206225189981, 0.0156474674491344, 0.0102810737264621, 0.0067272913372286), chpts = c(0L, 0L, 0L, 1L, 3L, 1L, 5L, 1L, 1L, 1L, 3L, 3L, 5L, 5L, 5L, 5L, 7L, 7L, 7L, 7L, 9L, 9L, 9L, 11L, 11L, 13L), M = c(5L, 30L), interval = c(-9, 11), mloglik = 42.330823537089, convergence = 0, ic = list(BIC = c(24.5019308079931, 27.5020279594087, 29.7438082116967, 31.0249946087796, 33.2658574545731, 33.5153424741971, 35.755444046287, 37.2829045542053, 35.6103492545364, 38.6208157877329, 34.9917247105547, 39.6156926690881, 35.9845312649796, 36.4627890289468, 36.7186874095177, 37.0391956348825, 37.2540399790624, 37.4659439196569, 37.6337328524481, 35.813860913558, 35.9337307689687, 36.0166305254596, 36.0902963767771, 36.1357311249942, 36.1661251401665, 34.2332366040113 )), xNames = "y", data.type = "noisy"), class = "mable") ``` ```{r simulated-normal-data-plot, fig.align='center', fig.cap="Simulated Normal Data", fig.width=7, fig.height=7} op <- par(mfrow = c(2,2), lwd = 2) plot(decn, which = "likelihood") plot(decn, which = "change-point", lgd.x = "right") plot(xx<-seq(a, b, length = 100), yy<-dnorm(xx, mu, sig), type = "l", xlab = "x", ylab = "Density", ylim = c(0, max(yy)*1.1)) plot(decn, which = "density", add = TRUE, lty = 2, col = 2) # kernel density based on pure data lines(density(x), lty = 5, col = 4) legend("topright", bty = "n", lty = c(1,2,5), col = c(1,2,4), c(expression(f), expression(hat(f)), expression(tilde(f)[K]))) plot(xx, yy<-pnorm(xx, mu, sig), type = "l", xlab = "x", ylab = "Distribution Function") plot(decn, which = "cumulative", add = TRUE, lty = 2, col = 2) legend("bottomright", bty = "n", lty = c(1,2), col = c(1,2), c(expression(F), expression(hat(F)))) par(op) ``` ## Interval Censored Data When data are interval censored, the ``\texttt{interval2}'' type observations are of the form $(l,u)$ which is the interval containing the event time. Data is uncensored if $l = u$, right censored if $u =$ \texttt{Inf} or $u =$ \texttt{NA}, and left censored data if $l =0$. Let $f(t)$ and $F(t)=1-S(t)$ be the density and cumulative distribution functions of the event time, respectively, on $(0, \tau)$, where $\tau\le \infty$. If $\tau$ is unknown or $\tau=\infty$, then $f(t)$ on $[0,\tau_n]$ can be approximated by $fm(t; p)=\tau_n^{-1}\sum_{i=0}^m p_i\beta_{mi}(t/\tau_n)$, where $p_i\ge 0$, $i=0,\ldots,m$, $\sum_{i=0}^mp_i=1-p_{m+1}$, and $\tau_n$ is the largest observation, either uncensored time, or right endpoint of interval/left censored, or left endpoint of right censored time. So we can approximate $S(t)$ on $[0,\tau]$ by $Sm(t; p)=\sum_{i=0}^{m+1} p_i \bar B_{mi}(t/\tau)$, where $\bar B_{mi}(u)=1-\int_0^u\beta_{mj}(t)dt$, $i=0,\ldots,m$, is the beta survival function with shapes $(i+1,m-i+1)$, $\bar B_{m,m+1}(t)=1$, $p_{m+1}=1-\pi$, and $\pi=F(\tau_n)$. For data without right-censored time, $p_{m+1}=1-\pi=0$. The search for optimal degree $m$ among \texttt{M=c(m0,m1)} using the method of change-point is stopped if either \texttt{m1} is reached or the test for change-point results in a p-value \texttt{pval} smaller than \texttt{sig.level}. @Guan-SIM-2020 proposed a method, as a special case of a semiparametric regression model, for estimating $p$ with an optimal degree $m$. The method is implemented in R function \texttt{mable.ic()}. ### Example: The Breast Cosmesis Data Consider the breast cosmesis data as described in @Finkelstein-and-Wolfe-1985-biom is used to study the cosmetic effects of cancer therapy. The time-to-breast-retractions in months ($T$) were subject to interval censoring and were measured for 94 women among them 46 received radiation only ($X=0$) (25 right-censored, 3 left-censored and 18 interval censored) and 48 received radiation plus chemotherapy ($X=1$) (13 right-censored, 2 left-censored and 33 interval censored). The right-censored event times were for those women who did not experienced cosmetic deterioration. We fit the Breast Cosmesis Data as two-sample data separately. ```{r} head(cosmesis, 3) ``` ```{r results = "hide", eval=FALSE} bc.res0 <- mable.ic(cosmesis[cosmesis$treat == "RT",1:2], M = c(1,50), IC = "none") bc.res1 <- mable.ic(cosmesis[cosmesis$treat == "RCT",1:2], M = c(1,50), IC = "none") ``` ```{r echo=FALSE} bc.res0<-structure(list(m = 6L, mloglik = -63.8031655794969, tau.n = 48, tau = Inf, interval = c(0, 48), convergence = 1L, delta = 0.0452318935094813, lr = c(1.797653214647, 2.33522547947038, 3.85049561635444, 5.32264477270326, 6.76686805682158, 6.68839261253333, 5.54710898977501, 4.7171220648023, 4.31270548240009, 4.44524613565427, 4.71747217723365, 5.47775864353639, 5.77061880461072, 5.44348787393976, 5.13688286730415, 5.10918996332308, 5.33986236662574, 5.66263871075579, 5.53916267665075, 5.40126854231033, 5.22522481613409, 5.0955945137348, 4.85895738595319, 4.69084338490006, 4.16548894463274, 3.89016586588323, 3.66415533926587, 3.50039695698725, 3.16064116685934, 2.79489267555537, 2.41724483179608, 2.25781251231455, 2.10039544998944, 1.90524753484281, 1.67568418400231, 1.41612335718722, 1.19031626406561, 1.02897013202351, 0.9291442462955, 0.843169950636643, 0.712951620830879, 0.554875886926446, 0.407673558404831, 0.321960751223394, 0.310097411617388, 0.232784116364762, 0.182364773305355, 0.113090679275252, 0), p = c(0.000359923711055115, 0.309858448466314, 3.04750494256107e-23, 1.67251493580178e-19, 0.00440764730467407, 0.214093057356895, 0.0560988992762524, 0.41518202388481), M = c(1L, 50L), lk = c(-64.709663640355, -64.4760177529401, -64.3423824219475, -64.1490081036323, -63.9699645956084, -63.8031655794969, -63.7287943272788, -63.7210176268191, -63.7020575928541, -63.6600081663922, -63.5833523020403, -63.499108395227, -63.3861227472242, -63.305831793135, -63.2645442805545, -63.2235942678219, -63.1665501685051, -63.0949662585894, -63.0199399845413, -62.9719807002143, -62.9257685847454, -62.8825907837, -62.8376030291276, -62.7994602333017, -62.7581624187024, -62.7384399330872, -62.7050059082487, -62.6692768597209, -62.630317285939, -62.603054471946, -62.5785811371243, -62.5563089115923, -62.5194186484327, -62.4828773366503, -62.449632311125, -62.419648890746, -62.3929903499263, -62.3646569094906, -62.3314657506243, -62.2930799198645, -62.2537543395296, -62.2189602221154, -62.1875384628887, -62.1560051280178, -62.1181770988759, -62.072160952893, -62.0338361744811, -61.9930454059307, -61.9547124975021, -61.9206119896223), pval = c(1, 1, 1, 0.307356049843645, 0.624990055507922, 0.7768073291445, 0.713643373129939, 0.290003078374176, 0.200745867154769, 0.214326538313812, 0.286450421808513, 0.322331989826414, 0.396140269396927, 0.407063833797846, 0.361541062930716, 0.322797768530098, 0.308117590358187, 0.310595662434251, 0.31611107178967, 0.295271492563875, 0.275417683102587, 0.255484818519026, 0.239234102938603, 0.22014737812776, 0.205280735444488, 0.179784287430119, 0.164798695753463, 0.152641355691455, 0.143226670963006, 0.1300658023749, 0.117459054776074, 0.105661383363785, 0.0996882272692384, 0.0941723085803453, 0.0882863846496835, 0.0821590628590383, 0.0758842296607785, 0.0705988933080037, 0.0667677530176096, 0.0642189317893982, 0.0620212898998007, 0.0592187679356891, 0.056100749362483, 0.0532501402304651, 0.0515011137282283, 0.0509722803933899, 0.0494565951419158, 0.0483423585511763, 0.0469920744710006, 0.0452318935094813), ic = list(BIC = c(-68.5383050368441, -70.2189798476737, -71.9996652149257, -73.720611594855, -75.4558887850757, -73.3747690707196, -77.1290392149906, -77.1212625145309, -80.930943877055, -78.9745737523486, -82.7265592844857, -84.556636075917, -82.5293297296697, -80.534718077336, -78.5791098665108, -76.6238391555337, -86.1383985474396, -89.8954560340131, -78.3345055704977, -76.3722255879262, -78.2403341707017, -80.1114770679009, -85.8094514080622, -89.5999500087254, -85.730010797637, -87.6246090102663, -91.419816381917, -93.2984080316337, -91.3451277596072, -85.5749028508805, -91.2933916107925, -91.2711193852605, -87.4055877256119, -93.112008508563, -93.0787634830378, -91.1344593644142, -93.022121521839, -92.9937880814033, -96.7892383190262, -94.8365317900219, -94.7972062096869, -96.6767327905172, -96.6453110312905, -98.5280983946642, -96.5759496672777, -100.358574917784, -98.4059294411275, -98.3651386725771, -98.3268057641485, -94.4640638597796)), chpts = c(1L, 1L, 1L, 2L, 2L, 2L, 6L, 7L, 7L, 7L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L), xNames = "cosmesis[cosmesis$treat == \"RT\", 1:2]", data.type = "icen"), class = "mable") bc.res1<-structure(list(m = 3L, mloglik = -75.9639126133468, tau.n = 60, tau = Inf, interval = c(0, 60), convergence = 0L, delta = 0.00938045116286756, lr = c(3.47912110850533, 12.1343351064705, 9.37442015676128, 8.82385795132178, 9.91405632343216, 7.69613487735494, 5.20402369287645, 4.63360804198853, 2.33928524265641, 0), p = c(2.10858395908529e-60, 0.956034481977071, 2.8535377405612e-12, 0.0436760732213005, 0.00028944479877519), M = c(1L, 11L), lk = c(-80.867749746675, -78.1489599763581, -75.9639126133468, -75.9313482278459, -75.7285907062759, -75.4865347047459, -75.4455005550738, -75.4234956072416, -75.3339517958129, -75.3089250460338, -75.2873449104489), pval = c(1, 1, 1, 0.131715395076138, 0.136324840837406, 0.116536285631219, 0.0680914564604432, 0.0383575442495347, 0.0258641953898849, 0.0154622244746406, 0.00938045116286756), ic = list(BIC = c(-84.7389507575829, -82.020160987266, -81.7707141297087, -83.6737502496617, -83.4709927280916, -85.1645372320156, -87.0591035877975, -88.9726991454192, -90.8187558394444, -92.7293295951194, -94.6433499649883)), chpts = c(1L, 1L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), xNames = "cosmesis[cosmesis$treat == \"RCT\", 1:2]", data.type = "icen"), class = "mable") ``` As the warning message suggested, we check the \texttt{pval}. The \texttt{pval} when the search stopped is `r round(bc.res0$pval[length(bc.res0$pval)],4)`. ```{r Breast-Cosmesis-Data-plot, message=FALSE, warning=FALSE, fig.align='center', fig.cap="Breast Cosmesis Data", fig.width=7, fig.height=7, warning = FALSE} op <- par(mfrow = c(2,2),lwd = 2) plot(bc.res0, which = "change-point", lgd.x = "right") plot(bc.res1, which = "change-point", lgd.x = "right") plot(bc.res0, which = "survival", xlab = "Months", ylim = c(0,1), main = "Radiation Only") plot(bc.res1, which = "survival", xlab = "Months", main = "Radiation and Chemotherapy") par(op) ``` ## Multivariate Data A $d$-variate density $f$ on a hyperrectangle $[a,b]=[a_1,b_1]\times \cdots\times[a_d,b_d]$ can be approximated by a mixture of $d$-variate beta densities on $[a,b]$, $\beta_{mj}(x; a, b)=\prod_{i=1}^d\beta_{m_i,j_i}[(x_i-a_i)/(b_i-a_i)]/(b_i-a_i)$, with proportions $p(j_1,\ldots,j_d)$, $0\le j_i\le m_i, i=1,\ldots,d$. Because all the marginal densities can be approximated by Bernstein polynomials, we can choose optimal degree $m_i$ based on observations of the $i$-th component of $x$. For the $i$-th marginal density, an optimal degree is selected using \texttt{mable()}. Then fit the data using EM algorithm with the selected optimal degrees $m=(m_1,\ldots,m_d)$ to obtain a vector $p$ of the mixture proportions $p(j_1,\ldots,j_d)$, arranged in the column-major order of $j = (j_1,\ldots,j_d)$, $(p_0,\ldots,p_{K-1})$, where $K=\prod_{i=1}^d (m_i+1)$. The proposed method of @Wang-and-Guan-2019 is implemented by function \texttt{mable.mvar()}. ### Example: The Old Faithful Data ```{r} data(faithful) head(faithful, 3) a <- c(0, 40); b <- c(7, 110) ``` ```{r message=FALSE, warning=FALSE, results = "hide", eval=FALSE} #faith2 <- mable.mvar(faithful, M = c(60,30), interval = cbind(a,b)) faith2 <- mable.mvar(faithful, M = c(46,19), search =FALSE, interval = cbind(a,b)) ``` ```{r echo=FALSE} faith2<-structure(list(p = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 1.3221316673478e-228, 5.24126655208957e-148, 3.58593698500509e-138, 1.34605085143491e-193, 5.4302909694972e-309, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.65803359736336e-286, 1.5574570054808e-128, 9.34817911624198e-45, 8.01540120088229e-28, 5.26517363414655e-71, 2.09116125201502e-167, 1.21629215860908e-309, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5.8767900380512e-289, 2.93570388763052e-125, 5.20146812955496e-32, 0.223856895143672, 2.20844683438706e-25, 1.58828949496376e-95, 1.71725305697173e-202, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 9.98285084653229e-157, 4.725809947529e-50, 0.119372660378794, 0.0142457607851216, 3.78037263424733e-43, 8.47826694876534e-113, 2.37470405049902e-200, 3.47277882785513e-295, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 7.33787768636394e-223, 6.16623630470217e-100, 5.52250130206306e-31, 1.3935825599307e-06, 1.39279201030188e-16, 4.99328836992708e-50, 8.94698563326245e-96, 3.89289103167831e-144, 3.42901760653081e-189, 1.96873952262891e-230, 1.98602187270129e-272, 1.23516411460312e-322, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.06969132033503e-197, 2.47852153716239e-105, 1.95620185688318e-54, 1.59515745840986e-34, 5.49273871838908e-35, 2.676952692219e-45, 6.149732208264e-57, 9.71733040321565e-66, 4.48166700816958e-73, 9.20331769709128e-85, 2.07276059657562e-108, 7.17496400486674e-151, 8.66853598708211e-217, 3.4091381380599e-308, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 1.36006626817657e-238, 1.68979699872748e-161, 4.22395618684829e-114, 2.4811717663345e-86, 3.52780833110995e-68, 8.92352058203179e-52, 9.0306283940599e-34, 3.36297011213198e-16, 5.49982875865874e-05, 4.1138286720304e-07, 3.37028054592537e-28, 2.57877431475673e-71, 7.46827625032142e-137, 4.86364378634656e-223, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 2.42346124083514e-268, 2.69888542439505e-217, 1.99878363283798e-176, 9.88144660897455e-138, 7.79351187445989e-98, 9.95461799993446e-59, 5.3038070021248e-26, 7.32188081979319e-06, 0.00846812111984434, 1.21000679930579e-18, 9.00408076869284e-53, 5.4441892189117e-103, 8.88422004225226e-167, 1.3492664841149e-242, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 3.91289684308764e-315, 7.08060982664866e-253, 5.60072914786952e-190, 6.52801504042434e-132, 1.7985239610607e-84, 1.64281828456228e-51, 2.66614062400093e-34, 4.18433164506145e-32, 2.28246639499047e-43, 4.95112888119984e-67, 1.36475854164726e-103, 1.30240143751911e-155, 1.22099944059073e-227, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 1.12337969030237e-299, 3.09088515614429e-218, 4.85061963181019e-150, 2.13618681986012e-96, 5.22768992540722e-57, 2.83928001823328e-31, 1.91271890043528e-19, 2.45524418554909e-23, 1.76937033857996e-46, 7.65457956185778e-94, 2.26571155179032e-171, 1.51288392741403e-285, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 1.63947791121016e-281, 3.33876378615635e-193, 1.00430666006635e-120, 7.88360960574434e-64, 2.99358723109032e-23, 0.352356624260163, 0.28163581317857, 2.28943210390895e-28, 4.76633509105543e-88, 2.25135548688222e-186, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.15528238411785e-228, 5.89854706046106e-146, 9.31989092444444e-82, 1.37192733132859e-37, 6.62271837384485e-17, 3.1481632965328e-24, 8.19456495253776e-65, 2.14712758086369e-144, 4.65475330076341e-269, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 6.83703434018381e-287, 2.11118013509095e-200, 5.70064642117329e-134, 9.61090441342611e-91, 6.30368493838083e-75, 1.57078948224484e-91, 3.45933287162858e-146, 4.75311401934995e-245, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 2.78585214636477e-296, 8.31249127174436e-226, 4.51840266046455e-181, 3.5050315749115e-167, 3.66856307484763e-190, 8.92519060057451e-257, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 1.8065773079144e-310, 2.27865799488389e-298, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), mloglik = 507.561455546254, interval = structure(c(0, 40, 7, 110), dim = c(2L, 2L), dimnames = list(NULL, c("a", "b" ))), m = c(46L, 19L), xNames = c("eruptions", "waiting"), convergence = 0L, D = 3.47435467146505e-05, data.type = "mvar"), class = "mable") ``` ```{r Old-Faithful-Data-plot, fig.align='center', fig.cap="Density Estimate based on the Old Faithful Data", fig.width=6, fig.height=6, warning = FALSE} plot(faith2, which = "density") ``` The density surface for two-dimensional data can be plot using the \texttt{plot} method (see Figure \ref{fig:Old-Faithful-Data-plot}). The summarized results are given below. ```{r} summary(faith2) ``` For multivariate density, the model degree can be chosen based on marginal data. As in Section \ref{old-faithful-data-mme} we obtained $m_b= (78,28)$ and $\hat m=(122,34)$. ```{r old-faithful-joint-data, message=FALSE, warning=FALSE, results = "hide"} x<-faithful x1<-faithful[,1] x2<-faithful[,2] a<-c(0, 40); b<-c(7, 110) mb<-c(78,28) m1<-122 m2<-34 ``` ```{r old-faithful-joint-data-mable, message=FALSE, warning=FALSE, results = "hide", eval=FALSE} oldfaith1<- mable.mvar(faithful, M = mb, search =FALSE, interval = cbind(a,b)) oldfaith2<- mable.mvar(faithful, M = c(m1,m2), search =FALSE, interval = cbind(a,b)) ``` ```{r echo=FALSE} oldfaith1<-structure(list(p = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 2.25122442130665e-196, 3.90683685902832e-90, 2.54092315635501e-93, 3.42724540986185e-198, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 1.48219693752374e-323, 4.16222211609906e-114, 1.61909964781621e-07, 0.0109329192225244, 9.75780699952314e-90, 4.66646702043307e-258, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 1.48219693752374e-323, 2.0725039913165e-143, 4.173322528328e-25, 0.165862350440829, 8.44439974154824e-61, 5.27351200131856e-190, 1.48219693752374e-323, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 9.67664817298388e-185, 2.44540533358448e-49, 0.02230155174135, 1.05610947776941e-30, 7.61351637797493e-120, 2.34131329980738e-254, 9.88131291682493e-324, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 1.28704603865639e-234, 3.16742445724399e-79, 2.81874836487009e-06, 0.0456010555976217, 6.98910327994437e-52, 1.61913294830031e-140, 2.03440091350932e-253, 1.48219693752374e-323, 9.88131291682493e-324, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 1.48219693752374e-323, 6.90710771844596e-146, 2.00426774829538e-42, 0.110282080838145, 5.24489305215529e-10, 6.16653885458302e-54, 4.19737102148436e-121, 2.30689651180735e-201, 7.05770866651787e-288, 1.48219693752374e-323, 9.88131291682493e-324, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 3.69845745899559e-280, 1.79690860841452e-139, 4.24739727448326e-56, 1.24928914853546e-18, 1.89765114191323e-16, 5.21601723520415e-40, 3.68161291945324e-81, 1.83678330942014e-133, 2.07501622367666e-192, 1.99912033696786e-255, 1.44267168585644e-321, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 4.28695820239991e-319, 4.86997534483237e-187, 3.56816585072276e-100, 3.3374242394245e-51, 2.61947044385416e-33, 1.00512438336718e-39, 1.50183579934528e-63, 6.36939906472253e-98, 1.52431625973566e-136, 2.10124693822751e-174, 1.79023303874337e-208, 3.96379095261762e-238, 3.41465990887902e-265, 1.35164069965782e-293, 1.48219693752374e-323, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 9.88131291682493e-324, 1.14075132887648e-274, 5.850881003518e-181, 6.48551863564442e-126, 3.32736019748718e-103, 1.25399530859678e-104, 1.90687875737783e-120, 1.71803306459021e-140, 6.81039704049918e-156, 6.5408639540663e-161, 3.58849086946926e-154, 8.8990886872618e-139, 2.90909747984955e-121, 6.52570829059147e-110, 4.40736489767937e-113, 5.44214547696352e-138, 4.34617171663046e-190, 6.06951811915136e-273, 9.88131291682493e-324, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 9.88131291682493e-324, 1.48219693752374e-323, 3.45129504911654e-295, 1.52622256726772e-280, 1.70306114105489e-283, 6.25377006505271e-290, 2.96379946930982e-287, 2.46484693114734e-267, 1.88062968730828e-228, 2.88678306200082e-175, 7.39280386061108e-117, 6.06060302869177e-64, 1.82840826410893e-26, 3.76667620476766e-12, 4.6374760866578e-26, 9.09442405458105e-71, 6.94668399746503e-147, 1.14467710432066e-253, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 4.94065645841247e-324, 9.88131291682493e-324, 9.88131291682493e-324, 1.48219693752374e-323, 1.61220891316991e-263, 7.11235346820593e-167, 2.61465944512125e-86, 7.26481408663434e-30, 0.0381247552286833, 7.64086669649036e-05, 1.54909466342051e-37, 1.83837269192759e-98, 3.88853308133044e-185, 9.49567870592988e-295, 9.88131291682493e-324, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 1.48219693752374e-323, 5.93146698515005e-244, 1.51011002176599e-140, 1.46252788465513e-67, 7.57243704425736e-26, 4.75580975881424e-14, 1.7068674141294e-29, 7.24200920803359e-69, 1.16773809296744e-128, 1.87155995286579e-205, 3.84640618043906e-296, 1.48219693752374e-323, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 1.73445075792113e-289, 4.15140901721722e-172, 4.42533696531428e-89, 7.6330337662206e-38, 1.51410867550641e-14, 9.31345480653008e-15, 5.30318961664863e-34, 3.45750867422115e-68, 5.31280436730918e-114, 2.48827104826936e-169, 1.02672757405359e-233, 9.71581002814381e-309, 1.48219693752374e-323, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 2.64600034350161e-284, 1.34763914398222e-169, 9.55825442511616e-90, 2.67096434581941e-39, 2.46277991900535e-12, 0.00423634116886228, 1.2945834907389e-06, 2.98588418171207e-19, 2.96659267192499e-39, 3.92711801957411e-67, 7.70662579052916e-106, 1.12693383155376e-160, 6.44980003868106e-239, 1.48219693752374e-323, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 9.88131291682493e-324, 8.61169153620938e-273, 2.633461680471e-176, 3.37887277897486e-110, 6.61825666523089e-67, 2.57611752301548e-39, 1.50665868282675e-21, 1.06057905529395e-09, 0.0207219791641175, 0.396999708413068, 4.80455089145829e-10, 1.05874499092136e-35, 2.12076707820655e-86, 1.45182154698227e-171, 3.52931968049712e-301, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 9.88131291682493e-324, 1.48219693752374e-323, 1.49202076072269e-249, 6.10306230178103e-195, 4.4945406364702e-153, 2.07725623843256e-117, 5.54366902943125e-84, 6.0512633679732e-52, 1.18759671448733e-23, 0.000246993505272245, 0.180295114580285, 4.51040654905674e-24, 5.08519503578743e-82, 1.4924240502191e-185, 1.48219693752374e-323, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 1.48219693752374e-323, 5.7901614315831e-300, 3.56460689422475e-241, 6.33702419455204e-181, 2.44269592002106e-122, 2.67254812206156e-71, 1.0211201367627e-35, 8.16539626394499e-25, 1.53140543806553e-48, 3.75468982794533e-117, 6.40101732957270e-241, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 9.88131291682493e-324, 1.48219693752374e-323, 4.74544015866018e-291, 7.56912059721118e-201, 2.48183115683234e-124, 2.43041052285199e-70, 1.57621403969013e-48, 4.76537908647714e-69, 2.70818953777212e-142, 1.06215981155948e-278, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 9.88131291682493e-324, 1.49299935577756e-272, 1.6120734370589e-163, 7.50658673091742e-83, 5.21775339823355e-41, 3.97416253057112e-49, 1.20400629152363e-118, 5.31129157723333e-261, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 1.48219693752374e-323, 2.92388761733832e-186, 1.29302195181297e-71, 3.58733643458477e-06, 0.00431087678219732, 2.68432455526347e-75, 3.8554398558441e-235, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 2.86439089112478e-236, 8.52705457475289e-98, 2.94090147851487e-26, 1.19285944950346e-35, 8.39550008254342e-139, 1.48219693752374e-323, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 2.78078248787996e-300, 9.34066221270549e-259, 5.77404638981748e-318, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), mloglik = 544.851771022873, interval = structure(c(0, 40, 7, 110), dim = c(2L, 2L), dimnames = list(NULL, c("a", "b" ))), m = c(78L, 28L), xNames = c("eruptions", "waiting"), convergence = 1L, D = 1.81234600575398e-05, data.type = "mvar"), class = "mable") oldfaith2<-structure(list(p = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 1.48219693752374e-323, 4.17710051078279e-262, 6.91273686559664e-251, 1.77722342031501e-299, 9.88131291682493e-324, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 2.12503906235987e-244, 5.13586752871597e-115, 3.23368953795656e-49, 8.52369663520542e-42, 2.66811139584931e-87, 4.27133280158359e-180, 1.86299284142594e-314, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 1.48219693752374e-323, 7.67946031413103e-208, 2.41645088562135e-81, 1.25776903349016e-14, 0.117042409402987, 5.97034780114412e-36, 1.26681930390804e-110, 1.65628111825882e-218, 1.48219693752374e-323, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 5.99137559133206e-231, 3.74751259413953e-101, 3.95213120159438e-28, 2.58423554124475e-05, 1.75614012171923e-25, 9.06869697859756e-81, 9.1860513932555e-163, 1.05941395604907e-262, 1.48219693752374e-323, 9.88131291682493e-324, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 1.53604104950486e-255, 9.32497266695284e-121, 2.37434616263598e-40, 1.97746175228732e-07, 2.51986732431296e-14, 1.4287230118795e-52, 5.38526344950946e-113, 5.86783803915934e-186, 7.61482709727077e-262, 1.48219693752374e-323, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 1.16227943918473e-278, 1.07589121456217e-138, 1.0842709414888e-50, 1.39391541288016e-07, 0.143068877766985, 3.99185594490976e-24, 4.53217361768229e-66, 7.48870058477017e-118, 4.65644133905843e-170, 3.07973609211234e-214, 5.86143218534082e-244, 7.21137392860016e-256, 2.11786711807421e-250, 8.45209256156721e-232, 2.09060892386754e-207, 1.95084766995428e-186, 1.19953411139748e-178, 3.1429194006432e-193, 4.62910162385541e-238, 2.55303481832006e-319, 9.88131291682493e-324, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 3.96734713610521e-321, 2.38228321180123e-175, 2.24355490415988e-78, 4.3846460869131e-23, 0.0702859696738787, 9.67888791958767e-07, 2.40997958674504e-29, 2.69632510361597e-61, 1.59446325805512e-94, 6.92698137651745e-122, 6.02910979897504e-138, 8.71477109319053e-140, 1.37267910653887e-127, 9.6944592183305e-105, 2.43854267701181e-77, 3.02293777768737e-53, 6.15567925175383e-41, 1.82437771331487e-48, 1.11563716507946e-82, 9.87671113849043e-149, 8.08130007586177e-250, 9.88131291682493e-324, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 9.88131291682493e-324, 2.27260992644436e-255, 1.43481018338604e-146, 3.23114487353983e-75, 3.17972933845855e-34, 1.32881421624519e-16, 1.33715688132915e-15, 7.92227435647028e-25, 2.03761051305441e-38, 3.00589758604894e-51, 2.78929427490442e-59, 5.56710727041841e-60, 9.82591330073139e-53, 4.40516398799699e-39, 1.24871876558728e-22, 1.42047634602433e-08, 0.00592228525142569, 8.53319134591583e-12, 2.52255222245755e-40, 2.31504906819116e-92, 6.34969895878352e-170, 5.0332361741877e-273, 9.88131291682493e-324, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 9.88131291682493e-324, 1.57239558136003e-277, 8.04030989406859e-185, 3.6495070994978e-118, 4.60359749137387e-72, 7.68759955169194e-42, 1.40183442986619e-23, 5.68408093556966e-14, 3.81412965467973e-10, 1.40827665672717e-09, 3.63496727019572e-10, 1.89537964219262e-10, 6.58760313343344e-10, 3.38751862915194e-09, 3.68108597167384e-10, 1.41287895739583e-15, 1.18855681143491e-28, 2.11118615671052e-52, 4.63244909858007e-89, 1.22316997286364e-139, 1.00938837369137e-203, 1.55419588001709e-279, 1.48219693752374e-323, 9.88131291682493e-324, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 1.48219693752374e-323, 4.69996445570091e-271, 1.69743834279182e-190, 6.9311488624554e-126, 6.92861348848531e-76, 1.25100722293932e-39, 3.11827717371941e-16, 0.000196501130100752, 0.0196624571870383, 5.54431544116407e-08, 1.13274900201698e-18, 2.4895278961492e-32, 2.14565933146319e-47, 4.02187549962531e-63, 1.49378208236526e-79, 2.9838492047046e-97, 6.20099015693503e-117, 4.16375567078688e-139, 7.11025892182123e-164, 5.87977761985216e-191, 6.52077495723774e-220, 1.82122871875058e-250, 7.94970721055522e-283, 7.4656925765827e-318, 1.48219693752374e-323, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 1.48219693752374e-323, 6.47782782040771e-283, 1.07070157495328e-198, 3.71781106337979e-134, 7.69767306380937e-90, 3.01641226188984e-65, 1.37197166520161e-58, 8.44434558180855e-67, 8.88634883273715e-86, 7.80219878479903e-111, 2.83190328623085e-137, 5.42132814513969e-161, 6.05033907318413e-179, 2.16299418512248e-189, 7.05937673303753e-192, 5.05884491993285e-187, 3.1074911056016e-176, 1.84144111734259e-161, 5.04374853933844e-145, 1.6772498646085e-129, 1.24043564783444e-117, 3.63577119121158e-112, 1.14205661713341e-115, 2.5625464823286e-130, 1.00883685424895e-157, 7.67521283791138e-199, 5.102221844999e-254, 3.95252516672997e-323, 9.88131291682493e-324, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 9.88131291682493e-324, 1.14810081233224e-313, 4.27657805879465e-249, 4.4741655431411e-212, 1.63093986040163e-200, 1.87306576980595e-210, 3.02586254981935e-236, 5.83656920892198e-271, 4.49454983152284e-307, 1.48219693752374e-323, 1.48219693752374e-323, 1.48219693752374e-323, 1.48219693752374e-323, 4.40892571806051e-316, 2.05602821371351e-274, 1.06285271449812e-224, 7.93971606974809e-172, 1.27337412273315e-120, 1.86540905492213e-75, 5.50981310041468e-40, 5.38804296140319e-17, 2.55275816318401e-08, 7.678743338787e-15, 1.41278369166402e-36, 7.37409624047523e-73, 1.28785167775691e-122, 1.03178540402656e-184, 2.73754260049356e-258, 1.48219693752374e-323, 9.88131291682493e-324, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 7.07021540826054e-315, 1.58411136145934e-232, 2.9234734610678e-158, 4.0164073487785e-96, 9.49326191116509e-49, 1.85204781427906e-17, 0.0204927288534873, 0.00757955881882005, 2.44054714126782e-17, 3.53324599751617e-44, 9.8680647320182e-82, 8.47758482245494e-129, 7.52085885941726e-185, 4.05429579998359e-250, 1.48219693752374e-323, 9.88131291682493e-324, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 9.88131291682493e-324, 1.48219693752374e-323, 3.65366374167142e-254, 2.2196773098846e-179, 6.59262805838445e-122, 2.3511827615665e-81, 2.05032153219157e-56, 3.23335948412036e-45, 1.02889412188111e-45, 4.67372489253079e-56, 6.53642003061765e-75, 1.17932053650374e-101, 2.10962993008091e-136, 7.54163718195583e-180, 5.66026343216991e-233, 1.14117410364164e-296, 1.48219693752374e-323, 9.88131291682493e-324, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 9.88131291682493e-324, 1.48219693752374e-323, 2.26122559254531e-271, 3.95407442861869e-201, 1.17549892483025e-147, 7.46844044353035e-109, 1.65455547470637e-82, 1.23942801732281e-66, 9.81710784496418e-60, 5.92329245850707e-61, 4.56031091342159e-70, 2.47447789459766e-87, 3.137945124776e-113, 3.71984799440745e-148, 3.7223414431115e-192, 9.64299979705359e-245, 7.37121781182566e-305, 1.48219693752374e-323, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 9.88131291682493e-324, 7.06603626203638e-310, 6.97433425745339e-229, 3.00750744811762e-165, 1.60578275084925e-116, 2.10416913641613e-80, 5.64853357803269e-55, 8.1736992912681e-39, 5.14642288280168e-31, 4.38975907709066e-31, 9.24902332741767e-39, 8.39702738479118e-54, 8.43446517805984e-76, 4.56013248827827e-104, 1.2200285017748e-137, 2.12700034692696e-175, 2.74749403295783e-216, 1.3746465825351e-259, 3.47366347304267e-305, 1.48219693752374e-323, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 1.48219693752374e-323, 5.65899787659358e-272, 3.00811693268112e-191, 6.35223984652064e-129, 4.00467707187403e-82, 3.32850398560086e-48, 6.93253077015153e-25, 2.57280578983042e-10, 0.00471974764196129, 0.0562251154626174, 3.56253211141097e-06, 7.96825951479344e-15, 3.96021645000909e-27, 2.58308576613207e-42, 1.00747775060822e-59, 5.85094417159565e-79, 4.55694987967676e-100, 1.0621723534981e-123, 3.18581366350767e-151, 9.47653782470652e-185, 4.29181859475405e-227, 1.14379346497535e-281, 1.48219693752374e-323, 9.88131291682493e-324, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 9.88131291682493e-324, 9.59845488549511e-283, 2.19760475297773e-194, 6.56514565931101e-128, 1.00541222512914e-79, 2.48971594693105e-46, 1.74768488490829e-24, 2.54643345211426e-11, 0.00019915433886375, 0.0736052812652973, 0.0414854413686152, 0.000481553515000026, 7.35928734850515e-07, 4.72274585689514e-10, 1.92389265840523e-13, 3.06770949225973e-17, 3.91336525616327e-22, 2.20822612681487e-29, 7.17827539291604e-41, 3.98353522840695e-59, 2.74900101552637e-87, 4.99802008885799e-129, 1.84212925845006e-188, 4.87445587912573e-270, 1.48219693752374e-323, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 9.88131291682493e-324, 2.38721842890848e-298, 3.48450748134973e-209, 1.30642248804749e-144, 2.45311141312381e-100, 4.71712526243218e-72, 1.27309445255325e-55, 2.83854249212496e-47, 9.05757506983457e-44, 1.64350524613334e-42, 1.41617422755314e-41, 1.0838936735032e-39, 3.62518324309656e-36, 8.21690550221954e-31, 7.02453924259889e-24, 4.69298121348534e-16, 1.78505263737913e-08, 0.00881246744050754, 0.356385535332593, 1.91294974037059e-06, 5.70925130388257e-22, 1.13530768196752e-51, 6.1923720057585e-100, 1.66224419817663e-171, 2.22667926075833e-271, 9.88131291682493e-324, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 9.88131291682493e-324, 4.51715754323207e-282, 2.99231439884979e-220, 1.17041181013887e-179, 2.77441967501879e-155, 2.89702255708503e-142, 3.67600408893915e-136, 3.34225736375995e-133, 1.82678460912064e-130, 6.13742763416153e-126, 1.76474947474656e-118, 1.12039776466326e-107, 1.09838311099446e-93, 4.34565449512892e-77, 8.0253844042204e-59, 3.50481557066237e-40, 6.89036776108713e-23, 3.50951606837611e-09, 0.0665804518398999, 0.000154467783126929, 3.45377911606487e-21, 1.69784854586951e-56, 1.61872909472278e-115, 1.43329439800229e-203, 1.48219693752374e-323, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 1.48219693752374e-323, 1.48219693752374e-323, 1.52074718717217e-307, 1.78663474451378e-295, 1.0927701936996e-284, 1.86687364048039e-272, 3.61064850639161e-257, 3.44077380717756e-238, 1.37778557064771e-215, 7.33220215332177e-190, 9.76254779193691e-162, 4.39013109547058e-132, 5.7466890603338e-102, 8.57700478139646e-73, 1.69387029054522e-46, 1.05227623313241e-25, 8.23420162353439e-14, 5.64877199355068e-15, 5.24808217533763e-34, 3.35308279774547e-76, 3.86944329674407e-147, 1.68916570064933e-252, 9.88131291682493e-324, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 1.48219693752374e-323, 3.7015892252072e-319, 1.26111527879745e-276, 1.02705918422774e-233, 8.58269555348366e-191, 1.6282797321953e-148, 3.86558557026525e-108, 8.80689795751501e-72, 1.53792648259897e-42, 1.81002202417554e-24, 1.97884438307384e-22, 7.42109034381049e-42, 1.69528857105352e-88, 3.40412933770554e-168, 1.05189574683234e-286, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 1.48219693752374e-323, 1.48219693752374e-323, 1.65768968452242e-266, 1.17333477848427e-209, 1.76960952978511e-154, 1.287973106036e-103, 6.16492963737777e-61, 2.62152580618684e-31, 2.67628947557508e-20, 7.28780930298222e-34, 4.43461424068779e-78, 5.87240987823475e-159, 2.48368056922358e-282, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 9.88131291682493e-324, 9.88131291682493e-324, 1.79980124364772e-308, 5.86543551944327e-230, 1.58259651850187e-155, 3.03201961712785e-90, 4.58578605898647e-40, 1.91689252449237e-11, 5.0150947194884e-11, 1.98900240525269e-45, 4.21054386177112e-121, 2.86427340751391e-244, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 9.88131291682493e-324, 1.48219693752374e-323, 2.4686993253813e-245, 5.47290224361043e-148, 4.88944900910601e-70, 3.83658710891909e-19, 0.00706651300836731, 1.49991266242615e-28, 4.14003118992495e-103, 4.26869500326167e-233, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 9.88131291682493e-324, 3.48118663712807e-274, 8.81819257928783e-172, 4.6883510097538e-106, 4.7742668489057e-85, 3.50137323945985e-116, 2.94056640236557e-206, 1.48219693752374e-323, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 9.88131291682493e-324, 9.88131291682493e-324, 1.48219693752374e-323, 1.48219693752374e-323, 9.88131291682493e-324, 9.88131291682493e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 4.94065645841247e-324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), mloglik = 562.649123135864, interval = structure(c(0, 40, 7, 110), dim = c(2L, 2L), dimnames = list( NULL, c("a", "b"))), m = c(122L, 34L), xNames = c("eruptions", "waiting"), convergence = 1L, D = 8.19263040034263e-06, data.type = "mvar"), class = "mable") ``` ```{r old-faithful-joint-data-plot, message=FALSE, warning=FALSE, results = "hide", fig.align='center', fig.cap="The Old Faithful data", fig.width=6, fig.height=3} op<-par(mfrow=c(1,2), cex=0.7) plot(oldfaith1, which="density", contour=TRUE) plot(oldfaith2, which="density", contour=TRUE, add=TRUE, lty=2, col=2) plot(oldfaith1, which="cumulative", contour=TRUE) plot(oldfaith2, which="cumulative", contour=TRUE, add=TRUE, lty=2, col=2) par(op) ``` # Event Time Data with Covariates Let $T$ be an event time and $X$ be an associated $d$-dimensional covariate with distribution $H(x)$ on $\cal{X}$. We denote the marginal and the conditional survival functions of $T$, respectively, by $S(t) = \bar F(t)= 1- F(t)=\Pr(T > t)$ and $S(t | x) = \bar F(t | x)= 1- F(t | x)=\Pr(T > t | X=x).$ Let $f(t | x)$ denote the conditional density of a continuous $T$ given $X=x$. The conditional cumulative hazard function and odds ratio, respectively, $\Lambda(t | x)=-\log S(t | x)$ and $O(t | x)={F(t | x)}/S(t | x)$. We will consider the general situation where the event time is subject to interval censoring. The observed data are $(X, Y)$, where $Y=(Y_1,Y_2]$, $0\le Y_1\le Y_2\le\infty$. The reader is referred to @Huang-and-Wellner-1997 for a review and more references about interval censoring. Special cases are right-censoring $Y_2=\infty$, left-censoring $Y_1=0$, and current status data. Let $f_0(\cdot)=f(\cdot\mid x_0)$, $F_0(\cdot)=F(\cdot\mid x_0)$, and $S_0(\cdot)=S(\cdot\mid x_0)$ be the baseline density, cumulative distribution, and survival functions, respectively. If the baseline density $f_0$ has support or can be truncated by $[0,\tau]$ then $$f_0(t)\approx f_m(t; \bm p)=\frac{1}{\tau}\sum_{j=0}^m p_j\beta_{mj}\Big(\frac{t}{\tau}\Big),\quad t\in[0,\tau],$$ where $\bm p=\bm p(\bm x_0)=(p_0,\ldots, p_{m})^\top$, $\sum_{j=0}^{m}p_j=1$, $p_j\ge 0$, $\beta_{mj}(t)=(m+1){m\choose j}t^j(1-t)^{m-j}$. Thus we have approximations $$F_0(t)\approx F_m(t; \bm p)= \sum_{j=0}^{m} p_j B_{mj}\Big(\frac{t}{\tau}\Big),\quad t\in[0,\tau],$$ $$S_0(t)\approx S_m(t; \bm p)= \sum_{j=0}^{m} p_j \bar B_{mj}\Big(\frac{t}{\tau}\Big),\quad t\in[0,\tau],$$ where $B_{mj}(t)=\int_0^t\beta_{mj}(u)du$, and $\bar B_{mj}(t)=1-B_{mj}(t)$, $j=0,\dots,m$. ## Accelerated Failure Time Model The accelerated failure time (AFT) model can be specified as \begin{equation}\label{eq: AFT model in terms of density} f(t\mid x)=f(t\mid x;\gamma)=e^{\gamma^T x}f(te^{\gamma^T x}\mid 0),\quad t\in[0,\infty), \end{equation} where $\gamma\in \mathbb{G}\subset \mathbb{R}^d$. Let $\gamma_0\in \mathbb{G}$ be the true value of $\gamma$. The AFT model (\ref{eq: AFT model in terms of density}) is equivalent to $$S(t\mid x;\gamma)=S(te^{\gamma^T x}\mid 0),\quad t\in[0,\infty).$$ Thus this is actually a scale regression model. The AFT model can also be written as linear regression $\log(T)=-\gamma^T x +\varepsilon$. It is clear that one can choose any $x_0$ in $\mathcal{X}$ as baseline by transform $\tilde{x}=x-x_0$. If $f_0(t)=f(t\mid x_0)$ has support $[0,\tau)$, $\tau\le\infty$, then $f(t\mid x)$ has support $[0,\tau e^{-\gamma_0^T \tilde x})$. The above AFT model can also be written as $$f(t\mid x;\gamma)=e^{\gamma^T \tilde {x}} f_0(te^{\gamma^T \tilde {x}}), \quad S(t\mid x;\gamma)=S_0(te^{\gamma^T \tilde {x}}),$$ where $f_0(t)=f(t\mid x_0)$ and $S_0(t)=S(t\mid x_0)=\int_t^\infty f_0(u)du$. We choose $\tau>y_{(n)}= \max\{y_{i1}, y_{j2}: y_{j2}<\infty;\, i,j=1,\dots,n\}$ so that $S(\tau)$ and $\max_{x\in\mathcal{X}} S(\tau\mid x)$ are believed very small. @Guan-SIM-2023 proposed to approximate $f_0(t)$ and $S_0(t)$ on $[0,\tau]$, respectively, by \begin{align*} f_m(t; p)&=\frac{1}{\tau}\sum_{j=0}^m p_j\beta_{mj}\Big(\frac{t}{\tau}\Big),\quad t\in[0,\tau];\\ S_m(t; p)&= \sum_{j=0}^{m} p_j \bar B_{mj}\Big(\frac{t}{\tau}\Big),\quad t\in[0,\tau]. \end{align*} Then $f(t\mid x; \gamma)$ and $S(t\mid x; \gamma)$ can be approximated, respectively, by \begin{align}\nonumber f_m(t\mid x; \gamma, p)&=e^{\gamma^T {x}}f_m\Big(te^{\gamma^T {x}}; p\Big)\\ &= \frac{e^{\gamma^T {x}}}{\tau_0}\sum_{j=0}^m p_j\beta_{mj}\Big(e^{\gamma^T {x}}\frac{t}{\tau_0}\Big),\quad t\in[0,\tau_0e^{-\gamma^T {x}}];\\\nonumber S_m(t\mid x;\gamma,p)&=S_m\Big(te^{\gamma^T {x}}; p\Big)\\ &= \sum_{j=0}^{m} p_j\bar B_{mj}\Big(e^{\gamma^T { x}}\frac{t}{\tau_0}\Big),\quad t\in[0,\tau_0e^{-\gamma^T {x}}]. \end{align} @Guan-SIM-2023's proposal is implemented by functions \texttt{mable.aft()} and \texttt{maple.aft()}. ### Example: Breast Cosmesis Data ```{r results = "hide", warning=FALSE, message=FALSE, eval=FALSE} library(mable) g <- 0.41 # Hanson and Johnson 2004, JCGS aft.res<-mable.aft(cbind(left, right) ~ treat, data = cosmesis, M =c(1, 30), g=.41, tau =100, x0=data.frame(treat = "RCT")) ``` ```{r echo=FALSE} aft.res<- structure(list(tau.n = 60, tau = 100, xNames = "treatRCT", x0 = 1, egx0 = 1.78233232862897, coefficients = 0.577922803596607, SE = 0.138888991045412, z = 4.16104112533762, M = c(1L, 13L ), lk = c(-149.195918793891, -147.122312893587, -146.805263205135, -144.818686096867, -143.395154013088, -143.152935116853, -143.089245666885, -142.963504206661, -142.948600114667, -142.918672815138, -142.830835168889, -142.795663720243, -142.774833692847), lr = c(1.97816658449391, 1.22644620428721, 4.70424542755593, 11.4650706877448, 11.9783865073313, 10.0832933323407, 9.69514494410889, 7.01982970449579, 4.85182104972307, 4.16479172050901, 2.32465708600656, 0), m = 6L, mloglik = -143.152935116853, p = c(0.105019596347297, 0.81462133922546, 0.0803590644272423, 2.93130761545064e-17, 7.82383246269338e-28, 1.75136396871527e-34, 5.74593171631315e-37), pval = c(1, 1, 1, 0.286290085033634, 0.574140880339093, 0.424953874506396, 0.23121084068413, 0.150788821414346, 0.0778775395330669, 0.0410203813856562, 0.0283966810218848, 0.0161430173880066, 0.00896652722735714), chpts = c(1L, 1L, 1L, 2L, 2L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L), convergence = 0L, delta = 9.96236906303238e-08, model = "aft", callText = "cbind(left, right) ~ treat", data.name = "cosmesis", allvars = c("left", "right", "treat")), class = "mable_reg") ``` ```{r} summary(aft.res) ``` ```{r Breast-Cosmesis-Data-aft-plot, fig.align='center', fig.cap="AFT Model Fit for Breast Cosmesis Data", fig.width=7, fig.height=4, warning = FALSE} op <- par(mfrow = c(1,2), lwd = 1.5) plot(x = aft.res, which = "likelihood") plot(x = aft.res, y = data.frame(treat = "RT"), which = "survival", type = "l", col = 1, main = "Survival Function") plot(x = aft.res, y = data.frame(treat = "RCT"), which = "survival", lty = 2, col = 1, add = TRUE) legend("topright", bty = "n", lty = 1:2, col = 1, c("Radiation Only", "Radiation & Chemotherapy"), cex = .7) par(op) ``` Alternatively we can use \texttt{mable.reg()}. ```{r results = "hide", eval=FALSE} aft.res1 <- mable.reg(cbind(left, right) ~ treat, data = cosmesis, model='aft', M = c(1, 30), g=.41, tau=100, x0=data.frame(treat = "RCT")) ``` ## Proportional Hazards Model Consider the Cox proportional hazard regression model [@Cox1972] \begin{equation}\label{eq: Cox PH model-conditional survival function} S(t | x) =S(t | x; \gamma, f_0)=S_0(t)^{\exp(\gamma^T\tilde{x})}, \end{equation} where $\gamma\in \mathbb{G} \subset \mathbb{R}^d$, $\tilde{x} =x -x_0$, $x_0$ is any fixed covariate value, $f_0(\cdot)=f(\cdot | x_0)$ is the unknown baseline density and $S_0(t)=\int_t^\infty f_0(s)ds$. Define $\tau=\inf\{t: F(t | x_{0})=1\}$. It is true that $\tau$ is independent of $x_0$, $0<\tau\le\infty$, and $f(t | x)$ have the same support $[0,\tau]$ for all $x\in \cal{X}$. Let $(x_i, y_i)=(x_i, (y_{i1}, y_{j2}])$, $i=1,\dots,n$, be independent observations of $(X, Y)$ and $\tau_n\ge y_{(n)}=\max\{y_{i1}, y_{j2}: y_{j2}<\infty;\, i,j=1,\dots,n\}$. Given any $x_0$, denote $\pi=\pi(x_0)=1-S_0(\tau_n)$. For integer $m\ge 1$ we define $\mathbb{S}_m\equiv \{(u_0,\ldots,u_m)^T\in \mathbb{R}^{m+1}: u_i\ge 0, \sum_{i=0}^m u_i=1.\}.$ @Guan-SIM-2020 propose to approximate $f_0(t)$ on $[0,\tau_n]$ by $f_m(t; p) = \tau_n^{-1}\sum_{i=0}^{m} p_i\beta_{mi}(t/\tau_n)$, where $p=p(x_0)=(p_0,\ldots,p_{m+1})$ are subject to constraints $p\in \mathbb{S}_{m+1}$ and $p_{m+1}=1-\pi$. Here the dependence of $\pi$ and $p$ on $x_0$ will be suppressed. If $\pi< 1$, although we cannot estimate the values of $f_0(t)$ on $(\tau_n,\infty)$, we can put an arbitrary guess on them such as $f_m(t; p)= p_{m+1} \alpha(t-\tau_n)$, $t\in (\tau_n,\infty)$, where $\alpha(\cdot)$ is a density on $[0,\infty)$ such that $(1-\pi)\alpha(0)=(m+1)p_m/\tau_n$ so that $f_m(t; p)$ is continuous at $t=\tau_n$, e.g., $\alpha(t)=\alpha(0)\exp[-\alpha(0)t]$. If $\tau$ is finite and known we choose $\tau_n=\tau$ and specify $p_{m+1}=0$. Otherwise, we choose $\tau_n=y_{(n)}$. For data without right-censoring or covariate we also have to specify $p_{m+1}=0$ due to its unidentifiability. The above method is implemented in function \texttt{mable.ph()} which returns maximum approximate Bernstein likelihood estimates of $(\gamma, p)$ with an optimal model degree $m$ and a prespecified $m$, respectively. With a efficient estimate of $\gamma$ obtained using other method, \texttt{maple.ph()} can be used to get an optimal degree $m$ and a mable of $p$. The \texttt{plot} method \texttt{plot.mable\_reg()} for class \texttt{mable\_reg} object returned by all the above functions produces graphs of the loglikelihoods at $m$ in a set \texttt{M[1]:M[2]} of consecutive candidate model degrees, the likelihood ratios of change-point at $m$ in \texttt{(M[1]+1):M[2]}, estimated density and survival function on the truncated \texttt{support}=$[0,\tau_n]$. ### Example: Ovarian Cancer Survival Data The Ovarian Cancer Survival Data is included in package \texttt{survival}. ```{r} library(survival) futime2 <- ovarian$futime futime2[ovarian$fustat==0] <- Inf ovarian2 <- data.frame(age = ovarian$age, futime1 = ovarian$futime, futime2 = futime2) head(ovarian2, 3) ``` ```{r results = "hide", warning = FALSE, message=FALSE, eval=FALSE} ova<-mable.ph(cbind(futime1, futime2) ~ age, data = ovarian2, M = c(2,35), g = .16, x0=data.frame(age=35)) ``` ```{r echo=FALSE} ova<-structure(list(M = c(2L, 35L), lk = c(-88.7099064570546, -88.6620193993037, -88.2063522797698, -88.0131362460812, -87.8779536441607, -87.5655730107827, -87.5246642110846, -87.1569385722069, -87.0686242248424, -86.915670957794, -86.6788875253044, -86.666889215658, -86.5005315387765, -86.4065185043726, -86.3597169716833, -86.1821891001746, -86.1446034760707, -86.0258604593683, -85.9006410893274, -85.8657313162381, -85.7858791729902, -85.7358640386468, -85.7192529335192, -85.6912022391792, -85.6793708330369, -85.6675085454649, -85.653167565389, -85.6361928984706, -85.6152024909356, -85.5948501128815, -85.569569754034, -85.5370922983257, -85.5119851122809, -85.4940530541106 ), lr = c(0.205895446573663, 1.44263977655423, 1.86111822391648, 1.9015830678408, 3.44520273496437, 2.75568700548905, 5.19060446171997, 4.95634337665228, 5.50631946111153, 7.31828283106786, 6.17604096480808, 7.26803778310007, 7.40644622131155, 6.83142369309167, 8.62784339961351, 7.92405270252482, 9.00740911688008, 10.7164642427666, 10.0494398288182, 10.9728242154829, 11.062985391912, 9.83970731301719, 9.13143872371833, 7.77202866646937, 6.46826391705687, 5.31004650066532, 4.28819289941039, 3.43349694885334, 2.58648273185088, 1.91437372126593, 1.50030813623333, 0.886993260948897, 0), m = 23L, tau.n = 1227, tau = Inf, mloglik = -85.7358640386468, p = c(2.56693371155884e-77, 4.17866353417925e-18, 0.00113629152423077, 2.02023368088419e-09, 1.49592034685534e-12, 9.62714381851158e-13, 8.08250811847392e-11, 9.16635704056953e-09, 2.27033359829197e-07, 0.0139227433708654, 0.00162007955037951, 4.38457006844612e-07, 7.17550846126675e-08, 1.18432071956013e-08, 2.116007943926e-09, 5.84101246782658e-10, 3.03802613707742e-10, 2.68652002753205e-10, 3.13034257781426e-10, 3.97618558289567e-10, 5.36293202680539e-10, 8.31949616488662e-10, 1.52738578415643e-09, 2.89001554215835e-09, 0.983320115427131), x0 = c(`27` = 35), egx0 = 484.940411000966, SE = 0.0105170971976927, z = 16.7999262893574, coefficients = 0.176686457699245, convergence = 0L, xNames = "age", pval = c(1, 1, 1, 0.207546359169874, 0.362766381048207, 0.433455588982329, 0.494949581051942, 0.556381521983088, 0.620890131052961, 0.665350546115876, 0.689845769796123, 0.380072973092538, 0.767935580798054, 0.782057828497582, 0.63532367457914, 0.759226925738293, 0.628923385821675, 0.639927577497975, 0.65638732444678, 0.551679856793764, 0.515057274506269, 0.448194300584407, 0.356269076540544, 0.291260323725976, 0.189822051638857, 0.11149881774139, 0.0642389260124165, 0.0411652770436294, 0.0295197991622843, 0.0213389970123636, 0.0170326796719841, 0.0152088066513323, 0.0122032848360507, 0.00899958750363183), chpts = c(2L, 2L, 2L, 3L, 3L, 3L, 7L, 3L, 3L, 3L, 3L, 12L, 3L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 20L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L), delta = 0.00899958750363183, model = "ph", callText = "cbind(futime1, futime2) ~ age", data.name = "ovarian2", allvars = c("futime1", "futime2", "age")), class = "mable_reg") ``` ```{r} summary(ova) ``` ```{r ovarian-Data-plot, fig.align='center', fig.cap="Ovarian Cancer Data", fig.width=7, fig.height=7, warning = FALSE, message=FALSE} op <- par(mfrow = c(2,2)) plot(ova, which = "likelihood") plot(ova, which = "change-point") plot(ova, y=data.frame(age=60), which="survival", type="l", xlab="Days", main="Age = 60") plot(ova, y=data.frame(age=65), which="survival", type="l", xlab="Days", main="Age = 65") par(op) ``` Alternatively we can use \texttt{mable.reg()}. ```{r results = "hide", warning = FALSE, eval=FALSE} ova1 <- mable.reg(cbind(futime1, futime2) ~ age, data = ovarian2, M = c(2,35)) ``` ## Proportional Odds Model As an important alternative of the Cox proportional hazards regression model [@Cox1972] the proportional odds (PO) model [@McCullagh-JRSSB-1980; @Bennett-stats-Med-1983; @Bennett-app-stats-1983; @Pettitt-AS-1984] is defined by \begin{equation}\label{eq: proportional odds model} \frac{O(t|\bm x)}{O(t|\bm x_0)}= e^{\bm\gamma^\top\tilde{\bm x}} \end{equation} where $\tilde{\bm x}=\bm x-\bm x_0$ and $O(t|\bm x_0)$ is an unknown increasing baseline odds function on $(0, \infty)$ and $\bm\gamma$ are unknown regression coefficients. If the baseline density $f_0$ has support or can be truncated by $[0,\tau]$ then $$f_0(t)\approx f_m(t; \bm p)=\frac{1}{\tau}\sum_{j=0}^m p_j\beta_{mj}\Big(\frac{t}{\tau}\Big),\quad t\in[0,\tau].$$ Therefore $f(t\mid\bm x; \bm\gamma, f_0)$, $F(t\mid\bm x; \bm\gamma, f_0)$, and $S(t\mid\bm x; \bm\gamma, f_0)$ can be approximated, respectively, by \begin{align}\nonumber f_m(t\mid\bm x; \bm\gamma, \bm p) &= \frac{e^{\bm\gamma^\top \tilde{\bm x}}f_m(t;\bm p)} {[e^{\bm\gamma^\top\tilde{\bm x}}+(1-e^{\bm\gamma^\top\tilde{\bm x}})S_m(t;\bm p)]^{2}}= \frac{e^{\bm\gamma^\top\tilde{\bm x}}f_m(t;\bm p)} {[1+(e^{\bm\gamma^\top\tilde{\bm x}}-1)F_m(t;\bm p)]^{2}},\\\nonumber F_m(t|\bm x; \bm\gamma, \bm p) =& \frac{e^{\bm\gamma^\top\tilde{\bm x}} F_m(t;\bm p)} {e^{\bm\gamma^\top\tilde{\bm x}}+(1-e^{\bm\gamma^\top\tilde{\bm x}})S_m(t;\bm p)}=\frac{e^{\bm\gamma^\top\tilde{\bm x}} F_m(t;\bm p)} {1+(e^{\bm\gamma^\top\tilde{\bm x}}-1)F_m(t;\bm p)},\\\nonumber S_m(t\mid\bm x; \bm\gamma, \bm p)&=\frac{S_m(t;\bm p)} {e^{\bm\gamma^\top\tilde{\bm x}}+(1-e^{\bm\gamma^\top\tilde{\bm x}})S_m(t;\bm p)}=\frac{S_m(t;\bm p)} {1+(e^{\bm\gamma^\top\tilde{\bm x}}-1)F_m(t;\bm p)}. \end{align} ### Example: HIV Infection time Data The data set is contained in R package \texttt{Epi}. This is an interval-censored data. The semiparametric methods (package \texttt{icenReg}) are available for PH and PO (but not AFT) models based on interval-censored data. The estimated survival functions are plotted in Fig \ref{fig:HIV-infection-time-data-po}. ```{r hiv-data, message=FALSE, warning=FALSE} ## Interval-censored data on times to HIV infection. Carstensen (1996) ## see also Lawless and Babineau (2006) # inspection times(months): started in Dec. 31, 1980, it<-c(0, 12, 16, 27, 45, 76, 101, Inf) # observed frequencies d[l,r] in interval (it[l], it[r]], 0<=l