Нахождение точки изменения в данных из кусочно-линейной функции

10

Приветствую,

Я провожу исследование, которое поможет определить размер наблюдаемого пространства и время, прошедшее с момента Большого взрыва. Надеюсь, вы можете помочь!

У меня есть данные, соответствующие кусочно-линейной функции, для которой я хочу выполнить две линейные регрессии. Есть момент, когда наклон и точка пересечения меняются, и мне нужно (написать программу) найти эту точку.

Мысли?

rhombidodecahedron
источник
3
Какова политика кросс-постинга? Точно такой же вопрос был задан на math.stackexchange.com: math.stackexchange.com/questions/15214/…
mpiktas
Что плохого в выполнении простых нелинейных наименьших квадратов в этом случае? Я что-то упускаю из виду?
grg s
Я бы сказал, что производная целевой функции по параметру точки изменения довольно негладкая
Андре Хольцнер,
Наклон изменился бы настолько сильно, что нелинейные наименьшие квадраты не были бы краткими и точными. Что мы знаем, так это то, что у нас есть две или более линейные модели, поэтому мы должны попытаться извлечь эти две модели.
HelloWorld,

Ответы:

1

mcpПакет может сделать это. Скажите, что ваши данные

Во-первых, давайте смоделируем некоторые данные:

df = data.frame(x = 1:100,
                y = c(rnorm(40, 10 + (1:40)*0.5),
                      rnorm(60, 10 + 40*0.5 -8 + (1:60)*0.2)))

Теперь давайте посмотрим, сможем ли мы восстановить точку изменения на 40 (и значения параметров), используя mcp:

model = list(
  y ~ 1 + x,  # linear segment
  ~ 1 + x  # another linear segment
)
library(mcp)
fit = mcp(model, df)

Сюжет. Серые линии - это случайные ничьи из подгонки, показывающие, что она отражает тренд. Синяя кривая - приблизительное местоположение точки изменения:

введите описание изображения здесь

Давайте посмотрим оценки отдельных параметров. int_являются перехватами, x_являются наклонами на х и cp_точками изменения:

summary(fit)

Population-level parameters:
    name  mean lower upper Rhat n.eff
    cp_1 40.48 40.02 41.00    1  2888
   int_1 11.12  9.11 13.17    1   778
   int_2 21.72 20.09 23.49    1   717
 sigma_1  3.23  2.76  3.69    1  5343
     x_1  0.46  0.36  0.54    1   724
     x_2  0.21  0.16  0.26    1   754

Отказ от ответственности: я разработчик mcp.

Йонас Линделёв
источник
8

R пакет Strucchange может помочь вам. Посмотрите на виньетку, в ней есть хороший обзор, как решать подобные проблемы.

mpiktas
источник
6

Иксязнак равно(Икся,Yя)язнак равно1,,,,NJ2N-2{Икс1,,,,,ИксJ}{Икс(J+1),,,,,ИксN}J


источник
Я опубликовал ответ, основанный на вашем простом, но эффективном предложении.
HelloWorld
5

Это проблема обнаружения (в автономном режиме) точки изменения. Наше предыдущее обсуждение содержит ссылки на журнальные статьи и R-код. Сначала рассмотрим Барри и Хартиган «модель разделения продуктов», потому что она обрабатывает изменения наклона и имеет эффективные реализации.

Whuber
источник
3

Также сегментированный пакет помог мне с подобными проблемами в прошлом.

Миша
источник
К сожалению, пакету требуется начальное значение для точки останова.
HelloWorld,
Кроме того, segmentedне может моделировать изменения перехвата между сегментами - только перехват для первого сегмента.
Йонас Линделёв
2

Я основывался на ответе mbq, что искал все возможности. Кроме того, я делаю это:

  • Проверьте значимость двух кусочных моделей, чтобы убедиться, что коэффициенты значимы
  • Проверьте разницу с суммой квадратов невязок для полной модели
  • Подтвердите мою модель визуально (убедитесь, что это не чушь)

Зачем проверять на значимость? Это потому, что точка с минимальным SSE не имеет смысла, если любая из кусочных моделей очень плохо соответствует данным. Это может произойти для двух сильно коррелированных переменных без четкой точки останова, где меняются наклоны.

Давайте проверим этот простой подход с помощью простого тестового примера:

x <- c(-50:50)
y <- abs(x)
plot(x,y,pch=19)

введите описание изображения здесь

Точка останова, очевидно, равна нулю. Используйте следующий скрипт R:

f <- function(x, y)
{
    d <- data.frame(x=x, y=y)
    d <- d[order(x),]
    r <- data.frame(k=rep(0,length(x)-4), sums=rep(0,length(x)-4))

    plm <- function(i)
    {
        d1 <- head(d,i)
        d2 <- tail(d,-i)

        # Make sure we've divided the region perfectly        
        stopifnot(nrow(d1)+nrow(d2) == nrow(d))

        m1 <- lm(y~x, data=d1)
        m2 <- lm(y~x, data=d2)

        r <- list(m1, m2)
        r
    }

    lapply(2:(nrow(d)-3), function(i)
    {
        r$k[i-2] <<- d[i,]$x

        # Fit two piecewise linear models
        m <- plm(i)

        # Add up the sum of squares for residuals
        r$sums[i-2] <<- sum((m[[1]]$residuals)^2) + sum((m[[2]]$residuals)^2)
    })

    b <- r[which.min(r$sums),]    
    b
}

Установите кусочно-линейные модели для всех возможных комбинаций:

f(x,y)
   k sums
   0    0

Если мы проверим коэффициенты для двух оптимальных моделей, они будут очень значительными. Их R2 тоже будет очень высоким.

Привет мир
источник