Выбор плоских частей кривой

Я думаю, что лучшим подходом было бы определить сегменты, в которых изменения минимальны. С учетом этих данных

Unsubtracted.Weight <- c(2.519903, 2.480581, 2.453806, 2.440516, 2.439226, 2.434226, 2.428516, 2.424839, 2.422839, 2.421839, 
                         2.420419, 2.419258, 2.418258, 2.41729, 2.416645, 2.415677, 2.415097, 2.414194, 2.413032, 2.412516, 
                         2.412806, 2.411839, 2.410677, 2.409935, 2.408355, 2.407452, 2.406323, 2.405419, 2.404355, 2.403355, 
                         2.40271, 2.401839, 2.401, 2.400419, 2.400032, 2.399452, 2.39871, 2.397581, 2.39671, 2.395806, 
                         2.395097, 2.394387, 2.393613, 2.386, 2.367258, 2.347581, 2.324, 2.287097, 2.230484, 2.144806, 
                         2.016871, 1.846968, 1.639097, 1.408452, 1.172484, 0.960161, 0.873258, 0.873065, 0.873226, 0.874194, 
                         0.87529, 0.875452, 0.876613, 0.876258, 0.877032, 0.877355, 0.878129, 0.878645, 0.879774, 0.880194, 
                         0.880452, 0.881419, 0.882226, 0.882935, 0.883806, 0.884419, 0.885032, 0.885581, 0.886387, 0.887, 
                         0.887645)
Sample.Temperature <- c(29.82, 29.95, 30, 30, 36.48, 53.15, 69.83, 86.5, 103.17, 119.82, 136.49, 153.16, 169.82, 186.48, 203.15, 
                        219.83, 236.49, 253.15, 269.82, 286.48, 303.16, 319.82, 336.49, 353.15, 369.82, 386.49, 403.15, 419.82, 
                        436.48, 453.15, 469.82, 486.48, 503.15, 519.81, 536.48, 553.15, 569.81, 586.48, 600, 600, 600, 600, 600, 
                        600, 600.21, 616.59, 633.26, 649.93, 666.6, 683.26, 699.92, 716.6, 733.26, 749.93, 766.6, 783.26, 
                        799.92, 816.59, 833.25, 849.92, 866.59, 883.25, 899.92, 916.59, 933.25, 949.91, 966.57, 983.23, 999.68, 
                        999.99, 1000.01, 1000.01, 999.99, 1000, 999.99, 1000, 1000, 999.99, 1000, 1000)

length_to_use <- min(length(Unsubtracted.Weight), length(Sample.Temperature))
Unsubtracted.Weight <- Unsubtracted.Weight[1:length_to_use]
Sample.Temperature <- Sample.Temperature[1:length_to_use]

tga.data <- data.frame(Sample.Temperature, Unsubtracted.Weight)

plot(Unsubtracted.Weight ~ Sample.Temperature, data = tga.data, pch=19, ylim = c(-1,3))

затем вы можете идентифицировать сегмент, а затем выполнить необходимое извлечение и средние вычисления.

diffs <- abs(diff(tga.data$Unsubtracted.Weight))

threshold <- 0.01
flat_sections <- which(diffs < threshold)

flat_start <- c()
flat_end <- c()

i <- 1
while (i < length(flat_sections)) {
  start <- flat_sections[i]
  while (i < length(flat_sections) && flat_sections[i + 1] == flat_sections[i] + 1) {
    i <- i + 1
  }
  end <- flat_sections[i]
  flat_start <- c(flat_start, start)
  flat_end <- c(flat_end, end)
  i <- i + 1
}

flat_means <- sapply(1:length(flat_start), function(i) {
  mean(tga.data$Unsubtracted.Weight[flat_start[i]:flat_end[i]])
})

print(flat_means)

который дает

> print(flat_means)
[1] 2.4108065 0.8791627

Основываясь на простой идее сравнения абсолютных значений разностей с допуском, мы можем сделать:

Редактировать

2024-05-26

Завернут в функцию. В комплект входит scale, что делает настройку tol более надежной.

const1D = \(x, tol = .01) {
  stopifnot(is.numeric(x), is.numeric(tol))
  b = abs(diff(scale(x))) <= tol
  o = b[c(TRUE, !b[-length(b)] == b[-1L])]
  i = 1L + which(diff(b) != 0L)
  split(x, cumsum(seq_along(x) %in% i))[o]
}
vapply(const1D(tga.data$Unsubtracted.Weight), mean, numeric(1L))

Код

tol = .01
b = abs(diff(tga.data$Unsubtracted.Weight)) <= tol
o = b[c(TRUE, !b[-length(b)] == b[-1L])]
i = 1L + which(diff(b) != 0L)
l = with(tga.data, split(Unsubtracted.Weight, 
                         cumsum(seq_along(Unsubtracted.Weight) %in% i)))
vapply(l[o], mean, numeric(1L))
#>         1         3 
#> 2.4108065 0.8798155

Сюжет

with(tga.data, plot(x = Sample.Temperature, y = Unsubtracted.Weight,
                    col = c('red','blue')[b+1L], pch = ".", cex = 3L))

Примечание

Более сложные подходы будут включать в себя алгоритмы обнаружения точек изменения. Байесовские версии могут быть хорошей идеей. Для справки, взгляните на beast() или beast.irreg() из {Rbeast} . Вероятно, на первый взгляд список аргументов (и файл справки) покажется вам ошеломляющим. Если все ваши ряды данных представляют собой одиночные S-образные кривые ( сигмовидная(x) или, точнее, -сигмовидная(x) функции), это может быть стрельба по воробьям из пушек. Тем не менее, вам нужно будет найти способ исключить точки перегиба, поскольку у этих функций ровно одна.

Я думаю, что хороший подход к поиску точек останова — это сегментированная линейная модель, поэтому вот пример.

tga.data = data.frame(y_weight = c(2.519903, 2.480581, 2.453806, 2.440516, 2.439226, 2.434226, 2.428516, 2.424839, 2.422839, 2.421839, 2.420419, 2.419258, 2.418258, 2.41729, 2.416645, 2.415677, 2.415097, 2.414194, 2.413032, 2.412516, 2.412806, 2.411839, 2.410677, 2.409935, 2.408355, 2.407452, 2.406323, 2.405419, 2.404355, 2.403355, 2.40271, 2.401839, 2.401, 2.400419, 2.400032, 2.399452, 2.39871, 2.397581, 2.39671, 2.395806, 2.395097, 2.394387, 2.393613, 2.386, 2.367258, 2.347581, 2.324, 2.287097, 2.230484, 2.144806, 2.016871, 1.846968, 1.639097, 1.408452, 1.172484, 0.960161, 0.873258, 0.873065, 0.873226, 0.874194, 0.87529, 0.875452, 0.876613, 0.876258, 0.877032, 0.877355, 0.878129, 0.878645, 0.879774, 0.880194, 0.880452, 0.881419, 0.882226, 0.882935, 0.883806, 0.884419, 0.885032, 0.885581, 0.886387, 0.887, 0.887645), 
                      x_temp = c(29.82, 29.95, 30, 30, 36.48, 53.15, 69.83, 86.5, 103.17, 119.82, 136.49, 153.16, 169.82, 186.48, 203.15, 219.83, 236.49, 253.15, 269.82, 286.48, 303.16, 319.82, 336.49, 353.15, 369.82, 386.49, 403.15, 419.82, 436.48, 453.15, 469.82, 486.48, 503.15, 519.81, 536.48, 553.15, 569.81, 586.48, 600, 600, 600, 600, 600, 600, 600, 600.21, 616.59, 633.26, 649.93, 666.6, 683.26, 699.92, 716.6, 733.26, 749.93, 766.6, 783.26, 799.92, 816.59, 833.25, 849.92, 866.59, 883.25, 899.92, 916.59, 933.25, 949.91, 966.57, 983.23, 999.68, 999.99, 1000.01, 1000.01, 999.99, 1000, 999.99, 1000, 1000, 999.99, 1000, 1000))

library(segmented)

plot(tga.data$x_temp,tga.data$y_weight)
(tga_lm <- lm(y_weight ~ x_temp,data=tga.data))

(seg_tga_lm <- segmented(tga_lm,npsi = 2))

(x_break_points_found <- seg_tga_lm$psi[,2])

plot(seg_tga_lm)
abline(v=x_break_points_found)

library(tidyverse)

(summary_info <- 
mutate(tga.data,
       grp=cut(x_temp,breaks=c(
         -Inf,
         x_break_points_found,
         Inf),
               include.lowest=TRUE)) |> 
  group_by(grp) |> summarise(mean_y = 
                               mean(y_weight),
                              sd = sd(y_weight),
                              count = length(y_weight)))

summary_info$mean_y
# 2.4046728 1.5984056 0.8798155

Поэтому интересующие нас значения составляют 2,4046728 и 0,8798155.