==46881== Memcheck, a memory error detector ==46881== Copyright (C) 2002-2024, and GNU GPL'd, by Julian Seward et al. ==46881== Using Valgrind-3.24.0 and LibVEX; rerun with -h for copyright info ==46881== Command: /data/blackswan/ripley/R/R-devel-vg/bin/exec/R --vanilla ==46881== R Under development (unstable) (2025-12-17 r89192) -- "Unsuffered Consequences" Copyright (C) 2025 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. Natural language support but running in an English locale R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > pkgname <- "Dimodal" > source(file.path(R.home("share"), "R", "examples-header.R")) > options(warn = 1) > library('Dimodal') Loading required package: statmod > > base::assign(".oldSearch", base::search(), pos = 'CheckExEnv') > base::assign(".old_wd", base::getwd(), pos = 'CheckExEnv') > cleanEx() > nameEx("Dicpt") > ### * Dicpt > > flush(stderr()); flush(stdout()) > > ### Name: Dicpt > ### Title: Class methods for Dicpt Objects > ### Aliases: Dicpt print.Dicpt summary.Dicpt mark.Dicpt plot.Dicpt > ### Keywords: Dimodal Dicpt spacing changepoints > > ### ** Examples > > ## Not run because detectors may not be available. > ## Not run: m <- Dimodal(faithful$eruptions, Diopt.local(analysis='cpt')) > ## Not run: m$cpt > ## Not run: summary(m$cpt) > ## You want a tall plot window for this. > ## Not run: dev.new(width=4, height=12) ; plot(m$cpt) > > > > cleanEx() > nameEx("Didata") > ### * Didata > > flush(stderr()); flush(stdout()) > > ### Name: Didata > ### Title: Class methods for Didata Objects > ### Aliases: Didata print.Didata summary.Didata > ### Keywords: Didata Dimodal > > ### ** Examples > > ## We override the analysis option to avoid changepoints, > ## which may not be available. > m <- Dimodal(faithful$eruptions, Diopt.local(analysis=c('diw','lp'), diw.window=16)) ==46881== Invalid read of size 8 ==46881== at 0x123454D8: segment_midrun (packages/tests-vg/Dimodal/src/utility.c:389) ==46881== by 0x123454D8: impl_midq (packages/tests-vg/Dimodal/src/utility.c:213) ==46881== by 0x123454D8: C_midq (packages/tests-vg/Dimodal/src/utility.c:98) ==46881== by 0x4A714D: R_doDotCall (svn/R-devel/src/main/dotcode.c:760) ==46881== by 0x4A76A3: do_dotcall (svn/R-devel/src/main/dotcode.c:1437) ==46881== by 0x4E330C: bcEval_loop (svn/R-devel/src/main/eval.c:8132) ==46881== by 0x4F1317: bcEval (svn/R-devel/src/main/eval.c:7515) ==46881== by 0x4F1317: bcEval (svn/R-devel/src/main/eval.c:7500) ==46881== by 0x4F164A: Rf_eval (svn/R-devel/src/main/eval.c:1167) ==46881== by 0x4F33CD: R_execClosure (svn/R-devel/src/main/eval.c:2389) ==46881== by 0x4F4086: applyClosure_core (svn/R-devel/src/main/eval.c:2302) ==46881== by 0x4F1755: Rf_applyClosure (svn/R-devel/src/main/eval.c:2324) ==46881== by 0x4F1755: Rf_eval (svn/R-devel/src/main/eval.c:1280) ==46881== by 0x4F641C: do_set (svn/R-devel/src/main/eval.c:3581) ==46881== by 0x4F19B2: Rf_eval (svn/R-devel/src/main/eval.c:1232) ==46881== by 0x52810B: Rf_ReplIteration (svn/R-devel/src/main/main.c:264) ==46881== Address 0xf435160 is 0 bytes after a block of size 2,224 alloc'd ==46881== at 0x4843866: malloc (/builddir/build/BUILD/valgrind-3.24.0/coregrind/m_replacemalloc/vg_replace_malloc.c:446) ==46881== by 0x535B81: Rf_allocVector3 (svn/R-devel/src/main/memory.c:2894) ==46881== by 0x5B56FD: Rf_allocVector (svn/R-devel/src/include/Rinlinedfuns.h:609) ==46881== by 0x5B56FD: Rf_ExtractSubset (svn/R-devel/src/main/subset.c:134) ==46881== by 0x5B7B50: VectorSubset (svn/R-devel/src/main/subset.c:216) ==46881== by 0x5B7B50: do_subset_dflt (svn/R-devel/src/main/subset.c:861) ==46881== by 0x4D89BB: VECSUBSET_PTR (svn/R-devel/src/main/eval.c:6367) ==46881== by 0x4E581E: bcEval_loop (svn/R-devel/src/main/eval.c:8357) ==46881== by 0x4F1317: bcEval (svn/R-devel/src/main/eval.c:7515) ==46881== by 0x4F1317: bcEval (svn/R-devel/src/main/eval.c:7500) ==46881== by 0x4F164A: Rf_eval (svn/R-devel/src/main/eval.c:1167) ==46881== by 0x4F33CD: R_execClosure (svn/R-devel/src/main/eval.c:2389) ==46881== by 0x4F4086: applyClosure_core (svn/R-devel/src/main/eval.c:2302) ==46881== by 0x4F4A78: Rf_applyClosure (svn/R-devel/src/main/eval.c:2324) ==46881== by 0x53C32B: dispatchMethod (svn/R-devel/src/main/objects.c:473) ==46881== > m$data data source low-pass spacing with 41 (0.150) kaiser filter interval spacing with 16 (0.059) interval positions at interval end; shift by -8 vs. low-pass Information row valid range sd x 1 - 272 3.5 1.141 Di 2 - 272 0.25 0.02574 LP Di 22 - 252 0.25 0.02686 Diw 17 - 272 1.116 0.2214 > > > > cleanEx() > nameEx("Diflat") > ### * Diflat > > flush(stderr()); flush(stdout()) > > ### Name: Diflat > ### Title: Class methods for Diflat Objects > ### Aliases: Diflat print.Diflat summary.Diflat mark.Diflat > ### Keywords: Diflat Dimodal > > ### ** Examples > > ## We override the analysis option to avoid changepoints, > ## which may not be available. > m <- Dimodal(faithful$eruptions, Diopt.local(analysis=c('diw','lp'), diw.window=16)) > m$lp.flats statistics endpoints raw len hexcur 148 - 238 (4.083 4.65) 91 0.002212 probabilities endpoints plen pexcur pflat pass 148 - 238 0.0213 0.0000 0.0000 T (2) accept at 0.0500 0.0100 > summary(m$diw.flats) summary of flat endpoints raw pflat passing tests 211 - 240 (4.45 4.603) 0.0004 excursion > > > > cleanEx() > nameEx("Dimodal") > ### * Dimodal > > flush(stderr()); flush(stdout()) > > ### Name: Dimodal > ### Title: Detect modality in the spacing of data. > ### Aliases: Dimodal print.Dimodal summary.Dimodal plot.Dimodal > ### Keywords: Dimodal modality spacing > > ### ** Examples > > ## Not running cpt analysis because local setup is not known. > ## The interval spacing is noisy with the default options, so require a > ## larger peak height with a temporary value to Diopt. > oldopt <- Diopt(analysis=c('lp','diw'), diw.param=list(peak.fht=0.125)) > ## Run the analysis. > m <- Dimodal(faithful$waiting) > ## If printing the results, the interval spacing peaks have a probability > ## just under 0.05 but fail the acceptance levels. > summary(m) Setup data source low-pass spacing with 41 (0.150) kaiser filter interval spacing with 27 (0.100) interval Peaks in Low-Pass Spacing feature positions at middle of filter summary of maxima pos raw ppeak passing tests 98 (65.80) 0.0000 height model, excursion 185 (80.05) 0.6735 none Flats in Low-Pass Spacing feature positions at middle of filter summary of flat endpoints raw pflat passing tests 170 - 200 (78.56 81.4) 0.0031 excursion Peaks in Interval Spacing feature positions at end of interval summary of maxima pos raw ppeak passing tests 104 (63.57) 0.0452 none 114 (68.00) 0.0211 excursion Flats in Interval Spacing feature positions at end of interval no flats found > ## Details about the peaks in both spacings. > print(m, feature="peaks") Setup data source low-pass spacing with 41 (0.150) kaiser filter interval spacing with 27 (0.100) interval positions at interval end; shift by -13 vs. low-pass Information row valid range sd x 1 - 272 53 13.59 Di 2 - 272 2 0.4244 LP Di 22 - 252 2 0.3868 Diw 28 - 272 9 2.455 Peaks in Low-Pass Spacing feature positions at middle of filter location of maxima pos raw minima pos (raw) support pos 98 (65.80) 43 - 160 (53.25 - 77.78) 53 - 141 185 (80.05) 160 - 210 (77.78 - 82.19) 166 - 205 statistics pos ht hexcur 98 3.4890 0.26500 185 0.3193 0.02438 probabilities pos pht pexcur ppeak pass 98 0.0002 0.0000 0.0000 T (2) 185 0.6735 0.7026 0.6735 accept at 0.0100 0.0500 Peaks in Interval Spacing feature positions at end of interval location of maxima pos raw minima pos (raw) support pos 104 (63.57) 44 - 112 (51.27 - 66.33) 69 - 109 114 (68.00) 112 - 155 (66.33 - 76.29) 113 - 151 statistics pos hexcur nrun runlen runht 104 6 33 14 6 114 8 22 8 9 probabilities pos pexcur pnrun prunlen prunht ppeak pass 104 0.1350 0.4879 0.0452 0.3061 0.0452 114 0.0330 0.4097 0.1429 0.0211 0.0211 T (1) accept at 0.0500 0.0100 0.0100 0.0050 > ## We find one peak in both spacings, but only the low-pass is significant. > match.features(m) matching peaks within 10 points after centering intervals low-pass interval pos raw #pass pos raw #pass 98 (65.8) 2 114 ( 68) 1 matching flats each overlapping by 0.70 after centering intervals none $peak.lp2diw [1] NA NA 5 NA NA NA NA $peak.diw2lp [1] NA NA NA NA 3 NA NA $flat.lp2diw [1] NA $flat.diw2lp integer(0) > ## Three plots side by side. The limited resolution of the data is clear > ## in the interval spacing. > dev.new(width=12, height=4) ; plot(m) > ## Restore the old option values. Diopt(NULL) returns to defaults. > oldopt <- Diopt(oldopt) > > > > cleanEx() > nameEx("Diopt") > ### * Diopt > > flush(stderr()); flush(stdout()) > > ### Name: Diopt > ### Title: Options for Dimodal package > ### Aliases: Diopt Diopt.local > ### Keywords: Dimodal modality spacing > > ### ** Examples > > ## Use the triangular kernel spanning 10% of the data for the > ## low-pass analysis. Set the interval to 30 and override the > ## ripple spec in the interval spacing but not for the low-pass. > Diopt(lp.kernel='bartlett', lp.window=0.10, diw.window=30, diw.param=list(flat.fripple=0.02)) $lp.kernel [1] "kaiser" $lp.window [1] 0.15 $diw.window [1] 0.1 $diw.param $diw.param$peak.fht [1] 0.125 > ## Diopt()$lp.kernel also accesses the new value. > Diopt('lp.kernel', 'diw.window') $lp.kernel [1] "bartlett" $diw.window [1] 30 > ## Temporarily override the kernel, without changing default. > Diopt.local(list(lp.kernel='Hanning'))$lp.kernel [1] "hanning" > Diopt()$lp.kernel [1] "bartlett" > ## Reset. > Diopt(NULL) Spacing to Analyze: analysis lp, diw, cpt General Data Preparation (data.*): midq 0 Low-Pass Filter Options (lp.*): kernel bartlett window 0.1 tests ht, pkexcur, len, ftexcur param -none- Interval Spacing Options (diw.*): window 30 tests pkexcur, runht, nrun, runlen, ftexcur param flat.fripple=0.02 Peak Detector Options (peak.*): fht 0.05 frelht 0.15 fhtie 0.001 fhsupp 0.9 Flat Detector Options (flat.*): fripple 0.05 minlen 30 fminlen 0.05 noutlier 1 distrib logistic Excursion Test Options (excur.*): nrep 15000 ntop 8 seed 0 Permutation (runht) Test Options (perm.*): nrep 5000 seed 0 Significance Thresholds (alpha.*): ht 0.01 pkexcur.lp 0.05 pkexcur.diw 0.05 len 0.05 ftexcur.lp 0.01 ftexcur.diw 0.01 runht 0.005 nrun 0.01 runlen 0.01 alpha TRUE Changepoint Detector Options (cpt.*): fncpt.max 0.05 qvote (0.1, 0.9) timeout 2 fsep 0.02 sep 10 libsep 2 libs -astsa, -ecp Tracking Options (track.*): maxwindow 0.4 Display Options: palette Dark 2 digits 4 Display Options - Color (colID.*): peak 8 flat 5 cpt 2 data 7 filter 4 hist 6 cdf 1 > > > > cleanEx() > nameEx("Dipeak") > ### * Dipeak > > flush(stderr()); flush(stdout()) > > ### Name: Dipeak > ### Title: Class methods for Dipeak Objects > ### Aliases: Dipeak print.Dipeak summary.Dipeak mark.Dipeak > ### Keywords: Dipeak Dimodal > > ### ** Examples > > ## We override the analysis option to avoid changepoints, > ## which may not be available. > m <- Dimodal(faithful$eruptions, Diopt.local(analysis=c('diw','lp'), diw.window=16)) > summary(m$lp.peaks) summary of maximum pos raw ppeak passing tests 99 (3.317) 0.0000 height model, excursion > m$diw.peaks location of maxima pos raw minima pos (raw) support pos 107 (3.317) 36 - 126 (1.852 - 3.784) 68 - 120 128 (3.822) 126 - 160 (3.784 - 4.104) 127 - 150 181 (4.258) 160 - 221 (4.104 - 4.506) 165 - 219 253 (4.707) 221 - 257 (4.506 - 4.750) 234 - 256 statistics pos hexcur nrun runlen runht 107 1.000 59 5 17 128 0.201 23 4 9 181 0.067 41 6 4 253 0.083 21 5 5 probabilities pos pexcur pnrun prunlen prunht ppeak pass 107 0.0005 0.3753 0.3233 0.0164 0.0005 T (1) 128 0.2922 0.5090 0.3475 0.0956 0.0956 181 0.9463 0.5477 0.1548 0.2351 0.1548 253 0.7274 0.1442 0.2113 0.2219 0.1442 accept at 0.0500 0.0100 0.0100 0.0050 > > > > cleanEx() > nameEx("Ditest") > ### * Ditest > > flush(stderr()); flush(stdout()) > > ### Name: Ditest > ### Title: Class Methods for Ditest Objects > ### Aliases: Ditest print.Ditest summary.Ditest > ### Keywords: Ditest > > ### ** Examples > > ptst <- Dipeak.test(1:4, 250, 0.15, 'bartlett') > ptst low-pass spacing height model evaluated with pinvgauss(ht, corrht, mu, lambda) ht corrht p.value 1 1.079774 0.5748570 2 1.790633 0.8970272 3 2.969477 0.9917696 4 4.924401 0.9998737 with mu=1.10256, lambda=4.661764 for features less than or equal to ht model parameters n=250, flp=0.15, filter=bartlett > ## Summary results. > summary(ptst) low-pass spacing height model ht 1.000000 2.0000000 3.0000000 4.0000000 p.value 0.574857 0.8970272 0.9917696 0.9998737 > > > > cleanEx() > nameEx("Ditrack") > ### * Ditrack > > flush(stderr()); flush(stdout()) > > ### Name: Ditrack > ### Title: Class methods for Ditrack Objects > ### Aliases: Ditrack print.Ditrack summary.Ditrack plot.Ditrack > ### Keywords: Dimodal spacing > > ### ** Examples > > ## Not run because of the run time. > ## Two plots side by side. > > > > cleanEx() > nameEx("excursion_tests") > ### * excursion_tests > > flush(stderr()); flush(stdout()) > > ### Name: Dimodal Excursion and Permutation Tests > ### Title: Feature significance test using permutation or bootstrap > ### (excursion) sampling. > ### Aliases: Diexcurht.test Dipermht.test > ### Keywords: Dimodal modality low-pass spacing > > ### ** Examples > > ## Recommended number of excursions/permutations is larger. > ## 2000 reduces the run-time. > set.seed(2) ; x <- round(runif(50,-1,1) * 10) / 10 > Diexcurht.test(6.5, 30, x, 2000, TRUE, lower.tail=FALSE, seed=3) excursion height test evaluated with Diexcurht.test(ht, ndraw, nexcur, seed) ht ndraw nexcur seed p.value 6.5 30 2000 3 0.07825 for features greater than or equal to ht > xrun <- rle(sign(diff(x)))$values * rle(sign(diff(x)))$lengths > Dipermht.test(6.5, xrun, 2000, FALSE, 3) run height permutation test evaluated with Dipermht.test(ht, nperm, seed) ht nperm seed p.value 6.5 2000 3 0.437 for features greater than or equal to ht > > > > cleanEx() > nameEx("kirkwood") > ### * kirkwood > > flush(stderr()); flush(stdout()) > > ### Name: kirkwood > ### Title: Asteroid Orbital Axis Showing Kirkwood Gaps > ### Aliases: kirkwood > ### Keywords: datasets > > ### ** Examples > > > ## Load data. > data(kirkwood, package="Dimodal") > > ## The Ditrack and changepoints can take a while to run, so to pass > ## the CRAN time limit in examples they are marked don't run. To see > ## the complete example call example(kirkwood, run.dontrun=TRUE). > > ## Start from scratch. > ## Normally you would not need to save the results of the call, but > ## this clutters up the screen. > opt <- Diopt(NULL) > > ## Set up the analysis. > ## This is a large data set so spacing in features is small but > ## smooth, which allows tighter detector parameters. Values depend > ## on the range of the spacing and are found by trial and error. > ## Leave off changepoints for repeatability in this example. > ## Set RNG seeds for repeatability. > ## Using bars to mark flats is easier to read with dense data. > opt1 <- Diopt(peak.fht=0.015, flat.fripple=0.0075, analysis=c("lp", "diw"), + excur.seed=3, perm.seed=5, mark.flat="bar") > > ## Use Ditrack to see where passing features (filled in dots) appear. > ## Interval spacing is the same cut-off. Wrapped in \donttest because of > ## the run time. > opt1 <- c(opt1, Diopt(lp.window=0.05, diw.window=0.05)) > > ## Run analysis. > m <- Dimodal(kirkwood) Warning in check.modelargs(n, flp, lower.tail) : models used outside recommended range (flp: 0.05-0.30, n 60-500) Warning in check.modelargs(n, flp, lower.tail) : models used outside recommended range (flp: 0.05-0.30, n 60-500) > > ## Summarize the data and spacing. > m$data data source low-pass spacing with 105 (0.050) kaiser filter interval spacing with 105 (0.050) interval positions at interval end; shift by -52 vs. low-pass Information row valid range sd x 1 - 2093 3.293 0.361 Di 2 - 2093 0.3187 0.01331 LP Di 54 - 2041 0.09225 0.00235 Diw 106 - 2093 1.263 0.1143 > > ## Print features that exist in both low-pass and interval spacing. > ## Allow large distance between peaks because of data set size. > mtch <- match.features(m, near=30) matching peaks within 30 points after centering intervals low-pass interval pos raw #pass pos raw #pass 141 (2.296) 1 198 (2.304) 1 378 (2.487) 1 419 (2.471) 1 1021 (2.848) 1 1101 (2.864) 1 1162 (2.944) 1 1189 (2.921) 2 1385 (3.056) 1 1444 (3.061) 3 matching flats each overlapping by 0.70 after centering intervals low-pass interval pos raw #pass pos raw #pass 850 - 959 (2.755 2.789) 1 886 - 1008 (2.749 2.789) 1 1519 - 1793 (3.118 3.173) 1 1564 - 1921 (3.116 3.191) 1 > > ## Gap distance (in AU), dropping the extra LP peak. > ## There are gaps at 2.30, 2.48, 2.86, 2.93, and 3.06 AU. > a_gap <- (select.peaks(m$lp.peaks)$x[-3L] + select.peaks(m$diw.peaks)$x) / 2 > a_gap [1] 2.300249 2.478981 2.855980 2.932442 3.058457 > > ## Orbital resonance with Jupiter, which has orbital axis of > ## 5.201 AU, from Kepler's third law. The corresponding > ## resonance is 3.40 (ratio 7:2 or 10:3), 3.04 (3:1), 2.46 (5:2), > ## 2.36 (7:3), and 2.22 (9:4). Simulations confirm all but the > ## first and last. > resonance <- (5.201 / a_gap) ^ (3/2) > resonance [1] 3.399916 3.038928 2.457523 2.362034 2.217566 > > ## The flats include families of asteroids, although other > ## considerations must be taken into account to find true > ## family members (composition of asteroid, orbital > ## eccentricity and inclination). > ## The second range at 2.76 AU includes the Ceres family, > ## the third at 3.13 AU the Themis. The first range at > ## 2.64 AU includes the Eunomia and Prosperina families. > m$diw.flats statistics endpoints raw hexcur 541 - 767 (2.592 2.690) 0.009402 886 - 1008 (2.749 2.789) 0.009502 1564 - 1921 (3.116 3.191) 0.009624 probabilities endpoints pexcur pflat pass 541 - 767 0.0001 0.0001 T (1) 886 - 1008 0.0060 0.0060 T (1) 1564 - 1921 0.0000 0.0000 T (1) accept at 0.0100 > > ## 3 graphs side-by-side with results. > dev.new(width=12, height=4) ; plot(m) dev.new(): using pdf(file="Rplots1.pdf") > > ## The full test results of the low-pass peaks. > m$lp.peaks location of maxima pos raw minima pos (raw) support pos 141 (2.296) 71 - 230 (2.228 - 2.379) 87 - 208 378 (2.487) 230 - 663 (2.379 - 2.668) 256 - 500 747 (2.710) 663 - 937 (2.668 - 2.785) 679 - 902 1021 (2.848) 937 - 1076 (2.785 - 2.884) 960 - 1064 1162 (2.944) 1076 - 1290 (2.884 - 3.017) 1111 - 1243 1385 (3.056) 1290 - 1724 (3.017 - 3.160) 1307 - 1524 statistics pos ht hexcur 141 0.8660 0.0007238 378 1.0754 0.0009012 747 0.2723 0.0002269 1021 0.7334 0.0006079 1162 1.0858 0.0009114 1385 0.5757 0.0004802 probabilities pos pht pexcur ppeak pass 141 0.0006 0.1341 0.0006 T (1) 378 0.0002 0.3075 0.0002 T (1) 747 0.0056 0.9352 0.0056 T (1) 1021 0.0010 0.1613 0.0010 T (1) 1162 0.0002 0.0601 0.0002 T (1) 1385 0.0019 0.7057 0.0019 T (1) accept at 0.0100 0.0500 > > ## In the Ditrack plot there are two peaks in the very lower right > ## corner. We need smaller filter kernels/intervals to find these. > ## The parameters already set will still apply. > opt3 <- Diopt(lp.window=0.015, diw.window=0.015, peak.fht=0.025, peak.frelht=0.10) > ## Running Dimodal with such small windows will generate warnings > ## about the model tests being pushed out of bounds. > ## This adds gaps at 1.93 (2:1) and 1.74 (7:4), both expected > ## from simulations. The 2:1 gap position shifts because the > ## spacing is so large that there is data point to anchor the > ## value and a second, smaller increase pulls the peak to 1.93. > m2 <- Dimodal(kirkwood) ; mtch <- match.features(m2, near=30) Warning in check.modelargs(n, flp, lower.tail) : models used outside recommended range (flp: 0.05-0.30, n 60-500) Warning in check.modelargs(n, flp, lower.tail) : models used outside recommended range (flp: 0.05-0.30, n 60-500) matching peaks within 30 points after centering intervals low-pass interval pos raw #pass pos raw #pass 152 (2.319) 0 149 (2.292) 0 382 (2.515) 0 391 (2.477) 1 1012 (2.833) 0 1034 (2.846) 1 1166 (2.965) 0 1183 (2.965) 1 1995 (3.348) 1 2015 (3.368) 1 2056 (3.633) 2 2068 (3.565) 1 matching flats each overlapping by 0.70 after centering intervals low-pass interval pos raw #pass pos raw #pass 411 - 520 (2.548 2.609) 1 422 - 539 (2.544 2.610) 1 543 - 654 (2.614 2.666) 1 549 - 673 (2.613 2.667) 1 797 - 986 (2.735 2.799) 1 791 - 1003 (2.724 2.800) 1 1216 - 1338 (3.004 3.022) 1 1227 - 1352 (3.002 3.022) 1 1484 - 1907 (3.107 3.201) 1 1490 - 1960 (3.104 3.214) 1 > dev.new(width=12, height=4) ; plot(m2) dev.new(): using pdf(file="Rplots2.pdf") > > ## To run changepoints, do this. > ## Using all supported libraries they surround the 3:1, 2:1, > ## and 7:4 gaps to the side of the peak, but not directly at > ## the feature. They also mark one side of the 5:2, 7:3, and > ## 9:4 gaps. > ## Not run: mcpt <- Dimodal(kirkwood, Diopt.local(analysis="cpt")) ; mcpt > ## Not run: dev.new(width=8, height=4) ; plot(mcpt) > ## Show the voting from all libraries. This is a tall graph. > ## Not run: dev.new(height=10, width=6) ; plot(mcpt$cpt) > > ## Restore default. > opt <- Diopt(opt) > > > > cleanEx() > nameEx("model_tests") > ### * model_tests > > flush(stderr()); flush(stdout()) > > ### Name: Dimodal Model Tests > ### Title: Significance models of features in the low-pass spacing. > ### Aliases: Dipeak.test Dipeak.critval Diflat.test Diflat.critval > ### Keywords: Dimodal modality low-pass spacing > > ### ** Examples > > pval <- Dipeak.test(0.25*(1:16), 200, 0.15,'kaiser', lower.tail=FALSE) > pval$p.value [1] 7.733783e-01 6.790977e-01 5.728901e-01 4.619492e-01 3.542696e-01 [6] 2.570857e-01 1.755776e-01 1.121774e-01 6.659554e-02 3.645105e-02 [11] 1.823023e-02 8.244314e-03 3.330622e-03 1.185178e-03 3.654521e-04 [16] 9.581747e-05 > ## Recovers pval. > Dipeak.critval(pval$p.value, 200, 0.15,'kaiser') [1] 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 2.25 2.50 2.75 3.00 3.25 3.50 3.75 [16] 4.00 > > pval <- Diflat.test(10*(1:12), 200, 0.15,'kaiser', 'logistic', lower.tail=FALSE) > pval$p.value [1] 7.654237e-01 5.308473e-01 2.962710e-01 7.336286e-02 4.099450e-02 [6] 2.441491e-02 1.319057e-02 5.939218e-03 2.124478e-03 6.346244e-04 [11] 1.742913e-04 4.656318e-05 > Diflat.critval(pval$p.value, 200, 0.15,'kaiser', 'logistic') [1] 10 20 30 40 50 60 70 80 90 100 110 120 > > > > cleanEx() > nameEx("runs_tests") > ### * runs_tests > > flush(stderr()); flush(stdout()) > > ### Name: Dimodal Runs Tests > ### Title: Tests of runs within subsets of data > ### Aliases: Dinrun.test Dirunlen.test > ### Keywords: Dimodal runs Markov Chain > > ### ** Examples > > ## The inner diff generates the spacing, the outer the signed difference. > xrun <- sign(diff( diff(sort( iris$Petal.Width )) )) > ## No epsilon needed for signed values. > Dinrun.test(xrun, 1, length(xrun), 0) Kaplanksy-Riordan test on number of runs evaluated with pnorm(nrun, Erun, Vrun, feps) nrun Erun Vrun feps p.value 59 63.58108 14.79334 0 0.1168144 for features less than or equal to nrun > Dirunlen.test(xrun, 1, length(xrun), 0) Markov chain model of longest run evaluated with test.markov.runlen(runlen, featlen, tmat, wt) runlen featlen p.value 27 148 0.01845608 Markov chain variables tmat, wt stored with result > > > > cleanEx() > nameEx("utilities") > ### * utilities > > flush(stderr()); flush(stdout()) > > ### Name: Dimodal Utility Functions > ### Title: Dimodal Utility Functions > ### Aliases: midquantile runs.as.rle select.peaks center.diw match.features > ### shiftID.place > ### Keywords: Dimodal runs peaks flats interval spacing > > ### ** Examples > > > m <- Dimodal(faithful$eruptions, Diopt.local(analysis=c('lp','diw'))) > # How many peaks were found? Use print.data.frame to see the full structure. > nrow(select.peaks(m$lp.peaks)) [1] 1 > nrow(select.peaks(m$diw.peaks)) [1] 2 > # Compare to m$diw.peaks. > m$diw.peaks location of maxima pos raw minima pos (raw) support pos 117 (3.450) 40 - 163 (1.848 - 4.093) 63 - 144 167 (4.121) 163 - 229 (4.093 - 4.517) 165 - 228 statistics pos hexcur nrun runlen runht 117 1.134 77 8 24 167 0.083 45 4 5 probabilities pos pexcur pnrun prunlen prunht ppeak pass 117 0.0019 0.1522 0.0489 0.0018 0.0018 T (2) 167 0.9485 0.7368 0.3460 0.2006 0.2006 accept at 0.0500 0.0100 0.0100 0.0050 > center.diw(m)$diw.peaks location of maxima pos raw minima pos (raw) support pos 104 (3.450) 40 - 163 (1.848 - 4.093) 63 - 144 154 (4.121) 163 - 229 (4.093 - 4.517) 165 - 228 statistics pos hexcur nrun runlen runht 104 1.134 77 8 24 154 0.083 45 4 5 probabilities pos pexcur pnrun prunlen prunht ppeak pass 104 0.0019 0.1522 0.0489 0.0018 0.0018 T (2) 154 0.9485 0.7368 0.3460 0.2006 0.2006 accept at 0.0500 0.0100 0.0100 0.0050 > # Flats do not match because the Diw feature only covers 50% of the LP. > match.features(m) matching peaks within 10 points after centering intervals low-pass interval pos raw #pass pos raw #pass 99 (3.317) 2 117 ( 3.45) 2 matching flats each overlapping by 0.70 after centering intervals none $peak.lp2diw [1] NA NA 3 NA NA $peak.diw2lp [1] NA NA 3 NA NA NA NA $flat.lp2diw [1] NA $flat.diw2lp [1] NA > > plot(sort(iris$Petal.Length)) > lines(midquantile(iris$Petal.Length, type=1L), col='red') > lines(midquantile(iris$Petal.Length, type=2L), col='blue') > lines(midquantile(iris$Petal.Length, type=3L), col='green') > lines(midquantile(iris$Petal.Length, type=4L), col='orange') > > # See the Dimodal.R source code for the use of shiftID.place. > > # To simplify the runs in the signed difference of the interval spacing > # runs.as.rle(Dimodal:::find.runs(m$data['signed',], 0.01), m$data['signed',]) > > > > ### *