Maximum Approximate Bernstein/Beta Likelihood Estimation in R package 'mable'"

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.

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,

data(Vaal.Flow)
head(Vaal.Flow, 3)

we want to estimate the density and the distribution functions of annual flow

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

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.

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.

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.

summary(vaal)

The mixing proportions, the coefficients of the Bernstein polynomials, $p$ can be obtained as either \texttt{p<-vaal\$p}

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}

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)
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]))
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)$.

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.

data(chicken.embryo)
head(chicken.embryo, 2)
a <- 0
b <- 21
day <- chicken.embryo$day
nT <- chicken.embryo$nT
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))
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")
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.

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.

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

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;
decn <- mable.decon(y, gn, interval = c(a,b), M = c(5, 50))
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")
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.

bcos=cosmesis
head(bcos, 3)
bc.res0 <- mable.ic(bcos[bcos$treat == "RT",1:2], M = c(1,50), IC = "none")
bc.res1 <- mable.ic(bcos[bcos$treat == "RCT",1:2], M = c(1,50), IC = "none")
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 = "bcos[bcos$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 = "bcos[bcos$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).

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

data(faithful)
head(faithful, 3)
a <- c(0, 40); b <- c(7, 110)
#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))
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")
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.

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)$.

x<-faithful                                                                                                                 
x1<-faithful[,1]                                                                                                            
x2<-faithful[,2]                                                                                                            
a<-c(0, 40); b<-c(7, 110)                                                                                                   
mb<-c(78,28)                                                              
m1<-122                                                            
m2<-34                                                          
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))
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")
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)={S(t | x)}/[{1-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.

Cox 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}.

library(survival)
futime2 <- ovarian$futime
futime2[ovarian$fustat==0] <- Inf
ovarian2 <- data.frame(age = ovarian$age, futime1 = ovarian$futime, futime2 = futime2)
head(ovarian2, 3)
ova<-mable.ph(cbind(futime1, futime2) ~ age, data = ovarian2, M = c(2,35), g = .16, x0=35)
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(5.00578691357198e-77, 4.63890599398755e-18, 0.00113628909769538, 
    2.06489209178884e-09, 1.57518116307075e-12, 1.01557329781422e-12, 
    8.37164666051139e-11, 9.31136041022891e-09, 2.27989567840045e-07, 
    0.0139227354457905, 0.00162008026644403, 4.39635387607504e-07, 
    7.23159357443982e-08, 1.20181755306809e-08, 2.16238912626703e-09, 
    6.00078203501842e-10, 3.12955915494075e-10, 2.76886839233189e-10, 
    3.2242592200995e-10, 4.09143725977586e-10, 5.51156911645021e-10, 
    8.5346035201065e-10, 1.56296029304163e-09, 2.94958649779085e-09, 
    0.9833201217674), x0 = 35, egx0 = 484.940948828593, SE = 0.0105224510461618, 
    z = 16.7913814577528, coefficients = 0.176686489386633, convergence = 1L, 
    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"), class = "mable_reg")
summary(ova)
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()}.

ova1 <- mable.reg(cbind(futime1, futime2) ~ age, data = ovarian2, M = c(2,35))

Accelerated Failure Time Model

The 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(t\mid 0)$ has support $[0,\tau_0)$, $\tau_0\le\infty$, then $f(t\mid x)$ has support $[0,\tau_0 e^{\gamma_0^T x})$. We define $\tau=\max{\tau_0 e^{\gamma_0^T x}: x\in\mathcal{X}}$ if $\tau_0<\infty$ and $\tau=\infty$ otherwise. The above AFT model can also be written as $$f(t\mid x;\gamma)=e^{-\gamma^T {x}} f_0(te^{-\gamma^T {x}}), \quad S(t\mid x;\gamma)=S_0(te^{-\gamma^T {x}}),$$ where $f_0(t)=f(t\mid 0)$ and $S_0(t)=S(t\mid 0)=\int_t^\infty f_0(u)du$.

As in Cox PH model, we choose $\tau_n>y_{(n)}= \max{y_{i1}, y_{j2}: y_{j2}<\infty;\, i,j=1,\dots,n}$ so that $S(\tau_n)$ and $\max_{x\in\mathcal{X}} S(\tau_n\mid x)$ are believed very small. Then we approximate $f_0(t)$ and $S_0(t)$ on $[0,\tau_n]$, respectively, by \begin{align} f_0(t)&\approx f_m(t; p)=\frac{1}{\tau_n}\sum_{j=0}^m p_j\beta_{mj}\Big(\frac{t}{\tau_n}\Big),\quad t\in[0,\tau_n];\ S_0(t)&\approx S_m(t; p)= \sum_{j=0}^{m} p_j \bar B_{mj}\Big(\frac{t}{\tau_n}\Big),\quad t\in[0,\tau_n], \end{align} where $S_m(\infty; p)=0$, and $p= (p_0,\ldots,p_{m})^T\in \mathbb{S}{m}$. 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_n}\sum{j=0}^m p_j\beta_{mj}\Big(e^{-\gamma^T {x}}\frac{t}{\tau_n}\Big),\quad t\in[0,\tau_ne^{\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_n}\Big),\quad t\in[0,\tau_ne^{\gamma^T {x}}]. \end{align} @Guan-2019-mable-aft's proposal is implemented by functions \texttt{mable.aft()} and \texttt{maple.aft()}.

Example: Breast Cosmesis Data

bcos2 <- data.frame(bcos[,1:2], x = 1*(bcos$treat == "RCT"))
g <- 0.41 # Hanson and Johnson 2004, JCGS
aft.res<-mable.aft(cbind(left, right) ~ x, data = bcos2, M =c(1, 30), g=.41,  
                   tau =100, x0=1)
aft.res<-structure(list(tau.n = 60, tau = 100, xNames = "x", M = c(1L, 
13L), coefficients = 0.577727331151052, x0 = 1, egx0 = 1.78198396581868, 
    SE = 0.129305349189261, lk = c(-149.195918910649, -147.122313202153, 
    -146.80526335503, -144.818686219186, -143.395154072282, -143.152935393756, 
    -143.089245862901, -142.963504497459, -142.948600408076, 
    -142.918673160554, -142.830835466953, -142.795664061396, 
    -142.77483409929), lr = c(1.97816638326472, 1.22644632716039, 
    4.70424614111728, 11.4650745849899, 11.9783884966994, 10.083296884476, 
    9.69514766293329, 7.01983199129244, 4.85182210899123, 4.16479533453579, 
    2.32466006469601, 0), m = 6L, mloglik = -143.152935393756, 
    p = c(0.105122858016893, 0.813986484599243, 0.0808906573838619, 
    1.95536547047634e-15, 6.80359923674277e-24, 3.3884962605114e-29, 
    3.58166309889493e-31), z = -4.46793063684818, pval = c(1, 
    1, 1, 0.286290100499535, 0.5741409239889, 0.424953704342698, 
    0.231210762842995, 0.150788738403658, 0.0778774886493685, 
    0.0410203684068813, 0.0283966794539574, 0.016143013481566, 
    0.00896652233227913), chpts = c(1L, 1L, 1L, 2L, 2L, 5L, 5L, 
    5L, 5L, 6L, 6L, 6L, 6L), convergence = 0L, delta = 0.00896652233227913, 
    model = "aft", callText = "cbind(left, right) ~ x", data.name = "bcos2"), class = "mable_reg")
summary(aft.res)
op <- par(mfrow = c(1,2), lwd = 1.5)
plot(x = aft.res, which = "likelihood")
plot(x = aft.res, y = data.frame(x = 0), which = "survival", type = "l", col = 1, 
    main = "Survival Function")
plot(x = aft.res, y = data.frame(x = 1), which = "survival", lty = 2, col = 1, add = TRUE)
legend("bottomleft", bty = "n", lty = 1:2, col = 1, c("Radiation Only", 
    "Radiation and Chemotherapy"), cex = .7)
par(op)

Alternatively we can use \texttt{mable.reg()}.

aft.res1 <- mable.reg(cbind(left, right) ~ x, data = bcos2, 'aft', M = c(1, 30), 
      g=.41, tau=100, x0=1)

Two-sample Data

In addition to the PH and AFT regression models with single binary covariate, other two-sample semiparametric models can be estimated using the maximum approximate Bernstein/Beta likelihhod method.

Density Ratio Model

Two-sample density ratio(DR) model (see for example, @QinZhang1997 and @Qin-and-Zhang-2005) is useful for fitting case-control data. Suppose that the densities $f_0$ and $f_1$ of $X_0$ and $X_1$, respectively, satisfy the following density ratio model (see @CHENG-and-CHU-Bernoulli-2004 and @Qin-and-Zhang-2005, for example) \begin{equation}\label{eq: exponential tilting model} f_1(x)=f(x;\bm\alpha)=f_0(x)\exp{\bm\alpha^T \tilde r(x)}, \end{equation} where $\bm\alpha=(\alpha_0,\ldots,\alpha_d)^T \in\mathcal{A}\subset R^{d+1}$, and $\tilde r(x)=(1, r^T(x))^T$ with linearly independent components. @Guan-arXiv-2021-mable-dr-model proposed to approximate the \emph{baseline} density $f_0$ by Bernstein polynomial and to estimate $f_0$ and $\bm\alpha$ using the maximum approximate Bernstein/Beta likelihood estimation. This method is implemented by function \texttt{mable.dr()} and \texttt{maple.dr()} in which $\bm\alpha$ is a given estimated value such as the one obtained by the logistic regression as described in @QinZhang1997 and @Qin-and-Zhang-2005.

Example: Coronary Heart Disease Data

# Hosmer and Lemeshow (1989):                                           
# ages and the status of coronary disease (CHD) of 100 subjects         
x<-c(20, 23, 24, 25, 26, 26, 28, 28, 29, 30, 30, 30, 30, 30, 32,        
32, 33, 33, 34, 34, 34, 34, 35, 35, 36, 36, 37, 37, 38, 38, 39,         
40, 41, 41, 42, 42, 42, 43, 43, 44, 44, 45, 46, 47, 47, 48, 49,         
49, 50, 51, 52, 55, 57, 57, 58, 60, 64)                                 
y<-c(25, 30, 34, 36, 37, 39, 40, 42, 43, 44, 44, 45, 46, 47, 48,        
48, 49, 50, 52, 53, 53, 54, 55, 55, 56, 56, 56, 57, 57, 57, 57,         
58, 58, 59, 59, 60, 61, 62, 62, 63, 64, 65, 69) 
a<-20; b<-70
regr<-function(x) cbind(1,x)  
chd.mable<-mable.dr(x, y, M=c(1, 15), regr, interval = c(a,b))   
chd.mable<-structure(list(regr = function (x) 
regr(a + (b - a) * x, ...), x = c(20, 23, 24, 25, 26, 26, 28, 
28, 29, 30, 30, 30, 30, 30, 32, 32, 33, 33, 34, 34, 34, 34, 35, 
35, 36, 36, 37, 37, 38, 38, 39, 40, 41, 41, 42, 42, 42, 43, 43, 
44, 44, 45, 46, 47, 47, 48, 49, 49, 50, 51, 52, 55, 57, 57, 58, 
60, 64), y = c(25, 30, 34, 36, 37, 39, 40, 42, 43, 44, 44, 45, 
46, 47, 48, 48, 49, 50, 52, 53, 53, 54, 55, 55, 56, 56, 56, 57, 
57, 57, 57, 58, 58, 59, 59, 60, 61, 62, 62, 63, 64, 65, 69), 
    interval = c(20, 70), mloglik = 24.412633993085, baseline = "Control", 
    lk = c(17.8495273773309, 20.0479815122578, 24.412633993085, 
    24.4214627038492, 24.461994343828, 24.5649515929184, 25.0704010073581, 
    25.6201041864968, 25.6774975123695, 25.7255171394501, 25.8355982818961, 
    25.9784597331011, 26.1533441209478, 26.224961081822), lr = c(1.46576799829151, 
    11.7440810573308, 9.06057505676728, 6.9459846642635, 5.39069589437123, 
    5.78863205111739, 7.32069295174204, 5.51798837590901, 3.80746410992697, 
    2.65902564975591, 1.79878802150627, 1.3393156127635, 0, 0
    ), p = c(0.0968615482255904, 0.898341880557055, 7.81704620921449e-07, 
    0.00479578951273396), wt = c(0.313928877027077, 1.04237811208442, 
    2.84669511165798, 6.91819858110639), m = 3L, alpha = c(-5.0400603627388, 
    0.111168522000774), chpts = c(1L, 1L, 1L, 3L, 3L, 3L, 3L, 
    3L, 3L, 3L, 3L, 3L, 3L, 3L), pval = c(1, 1, 1, 0.104528285150261, 
    0.0579937377184714, 0.0395265612490706, 0.064387242799094, 
    0.0741133539078582, 0.0471604715140564, 0.0301377062253065, 
    0.0212983347025477, 0.016077600053771, 0.0129017186124526, 
    0.00926163010372738), M = c(1L, 14L)), class = "mable_dr")
z<-seq(a,b,length=512)
f0hat<-dmixbeta(z, p=chd.mable$p, interval=c(a, b))
rf<-function(x) regr((x-a)/(b-a))
f1hat<-dtmixbeta(z, p=chd.mable$p, alpha=chd.mable$alpha, 
                 interval=c(a, b), regr=regr)
op<-par(mfrow=c(1,2),lwd=1.2, cex=.7, mar=c(5,4,1,1))
hist(x, freq=F, col = "light grey", border = "white", xlab="Age", 
  ylab="Density", xlim=c(a,b), ylim=c(0,.055), main="Control")
lines(z, f0hat, lty=1, col=1)
hist(y, freq=F, col = "light grey", border = "white", xlab="Age", 
  ylab="Density", xlim=c(a,b), ylim=c(0,.055), main="Case")
lines(z, f1hat, lty=1, col=1)
par(op)

Example: Pancreatic Cancer Biomarker Data

For the logarithmic levels of CA 19-9 of the Pancreatic Cancer Data (@Wieand-et-al-1989-bka), the control group is used as baseline because the optimal model degree is smaller that the one using case as baseline.

data(pancreas)
head(pancreas,3)
x<-log(pancreas$ca199[pancreas$status==0])
y<-log(pancreas$ca199[pancreas$status==1])
a<-min(x,y); b<-max(x,y)
M<-c(1,29)
regr<-function(x) cbind(1,x,x^2)                                            
m=maple.dr(x, y, M, regr=regr, interval=c(a,b), controls=mable.ctrl(sig.level=.001))$m
pc.mable<-mable.dr(x, y, M=m, regr=regr, interval=c(a,b),
                  controls=mable.ctrl(sig.level=1/length(c(x,y))))
#pc.mable   
m<-2
pc.mable<-structure(list(regr = function (x) 
regr(a + (b - a) * x, ...), x = c(3.3322045101752, 2.7408400239252, 
2.10413415427021, 1.22377543162212, 2.85070650150373, 2.72129542785223, 
3.49347265777133, 2.40694510831829, 4.47163879336357, 2.78501124223834, 
4.68120487226409, 1.7404661748405, 3.24259235148552, 3.44041809481544, 
3.07269331469012, 4.01818320125654, 2.17475172148416, 1.87180217690159, 
3.09557760852371, 2.66722820658195, 3.78872478908365, 1.30833281965018, 
2.05412373369555, 2.18605127673809, 2.89037175789616, 1.87180217690159, 
1.58923520511658, 2.34180580614733, 1.6094379124341, 1.66770682055808, 
1.87180217690159, 1.93152141160321, 2.10413415427021, 3.08190996979504, 
1.88706964903238, 2.02814824729229, 2.73436750941958, 4.08092154188996, 
1.62924053973028, 2.30258509299405, 1.66770682055808, 3.48431228837266, 
1.52605630349505, 1.93152141160321, 1.38629436111989, 1.2947271675944, 
2.05412373369555, 3.48124008933569, 2.4423470353692, 1.38629436111989, 
2.32238772029023), y = c(0.8754687373539, 6.57786135772105, 7.65286235670063, 
10.0858091093301, 7.44716835960004, 1.28093384546206, 6.25670927444083, 
7.37775890822787, 6.11809719804135, 4.69774936728118, 3.16547504814109, 
6.13988455222626, 9.19115755255941, 5.54126354515843, 4.07243972683405, 
5.41610040220442, 4.50092016461429, 3.91202300542815, 1.7227665977411, 
8.31139827843664, 6.38350663488401, 3.35340671782581, 8.72583205652757, 
6.99393297522319, 2.34180580614733, 3.30688670219091, 5.08759633523238, 
8.17751582384608, 2.68784749378469, 4.4224485491728, 5.8171111599632, 
4.01998014693324, 7.32646561384032, 1.3609765531356, 1.75785791755237, 
2.13416644136908, 5.88887795833288, 5.91079664404053, 9.01554129367111, 
3.67122451887522, 3.77276093809464, 5.88887795833288, 2.54944517092557, 
2.89037175789616, 9.16847616787748, 6.31896811374643, 4.09767235231478, 
3.08190996979504, 6.80239476332431, 1.88706964903238, 5.47646355193151, 
8.03915739047324, 8.09407314806935, 6.52502965784346, 4.44734610079452, 
9.23892782882809, 6.64639051484773, 5.51181454081044, 9.41897923708751, 
4.72827238312207, 6.98378996525813, 3.81990771652034, 7.39633529380081, 
4.37449836825309, 6.23048144757848, 8.06777619577889, 6.29526600143965, 
6.92853781816467, 5.45958551414416, 5.52545293913178, 8.05832730658096, 
6.17170059741091, 5.40267738187228, 2.75366071235426, 7.83991936001258, 
9.36134324551271, 7.50108212425987, 1.93152141160321, 1.41098697371026, 
2.74727091425549, 9.19217640134851, 7.30653139893951, 2.75366071235426, 
3.82428409112014, 2.05412373369555, 2.54944517092557, 4.61048901590065, 
5.4249500174814, 4.26127043353808, 7.82404601085629), interval = c(0.8754687373539, 
10.0858091093301), mloglik = 54.5988013637832, baseline = "Case", 
    lk = 54.5988013637832, lr = 0, p = c(0.13894796716688, 0.747617734547747, 
    0.113434298285373), wt = c(2.58081131881879, 0.841320451147648, 
    0.109450352149474), m = 2L, alpha = c(-0.134581496823654, 
    1.73897878071386, -0.443320912157217), chpts = 2L, pval = 1, 
    M = c(2L, 2L)), class = "mable_dr")
z<-seq(a,b,length=512)
# baseline is "case"
f1hat<-dmixbeta(z, p=pc.mable$p, interval=c(a, b))
rf<-function(x) regr((x-a)/(b-a))
f0hat<-dtmixbeta(z, p=pc.mable$p, alpha=pc.mable$alpha, 
                 interval=c(a, b), regr=regr)
op<-par(mfrow=c(1,2),lwd=1.2, cex=.7, mar=c(5,4,1,1))
hist(x, freq=F, col = "light grey", border = "white", xlab="log(CA19-9)", 
  ylab="Density", xlim=c(a,b),  main="Control")
lines(z, f0hat, lty=1, col=1)
hist(y, freq=F, col = "light grey", border = "white", xlab="log(CA19-9)", 
  ylab="Density", xlim=c(a,b), ylim=c(0,.5), main="Case")
lines(z, f1hat, lty=1, col=1)
par(op)

References



Try the mable package in your browser

Any scripts or data that you put into this service are public.

mable documentation built on Aug. 24, 2023, 5:10 p.m.