Экстраполяция в R

У меня есть эталонный набор данных, p,

p=structure(list(v = 0:26, t = c(Inf, 1.016, 0.568, 0.418666666666667, 
0.344, 0.2992, 0.269333333333333, 0.248, 0.232, 0.219555555555556, 
0.2096, 0.201454545454545, 0.194666666666667, 0.188923076923077, 
0.184, 0.179733333333333, 0.176, 0.172705882352941, 0.169777777777778, 
0.167157894736842, 0.1648, 0.162666666666667, 0.160727272727273, 
0.15895652173913, 0.157333333333333, 0.15584, 0.154461538461538
)), .Names = c("v", "t"), row.names = c(NA, -27L), class = "data.frame")

Какие p$v и p$t имеют следующее отношение:

t=(0.16*(0.75*v + 5.6))/v

мой второй набор данных - это измеренные данные, w, содержит ту же переменную, что и:

w=structure(list(v = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 
13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26), t = c(0.235291176470588, 
0.354020375722543, 0.310974343434343, 0.25272725498699, 0.20351968240702, 
0.163155804025208, 0.132330740162655, 0.108593108108108, 0.0859813015873016, 
0.0655131899302683, 0.0492580103144236, 0.0368029846567365, 0.030538003253355, 
0.0300744415648525, 0.0347586421891237, 0.0451097744360902, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA)), .Names = c("v", "t"
), row.names = c(NA, -27L), class = "data.frame")

Мои эталонные данные p соответствуют степенному закону, я хотел бы сделать прогноз на основе степенного закона, подходящего для данных, чтобы заменить NA в моих измеренных данных. Как я мог сделать это в R ?


r regression extrapolation nls
person user9112767    schedule 14.03.2018    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Как насчет чего-то подобного?

# Non-linear power law fit
fit.nls <- nls(
    t ~ a * v^b,
    data = w[-1, ],
    start = list(a = 0.4, b = -0.7),
    na.action = na.exclude);
# Linear fit with log-log transformation
fit.lm <- lm(
    log(t) ~ log(v),
    data = w[-1, ])

# Plot
w %>%
    rename(t.src = t) %>%
    mutate(
        t.pred.nls = predict(fit.nls, data.frame(v = v)),
        t.NA.pred.nls = ifelse(
            is.na(t.src),
            predict(fit.nls, data.frame(v = v)),
            t.src),
        t.pred.lm = exp(predict(fit.lm, data.frame(v = v)))) %>%
gather(source, t, 3:5) %>%
ggplot(aes(v, t, colour = source)) +
    geom_line() +
    geom_point(aes(v, t.src), colour = "black") +
    scale_colour_manual(values = c(
        "t.NA.pred.nls" = "black",
        "t.pred.nls" = "blue",
        "t.pred.lm" = "red"));

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

Черные точки показывают фактические измерения. Синяя кривая представляет собой соответствие нелинейной модели степенного закона, красная кривая - линейное соответствие после логарифмического преобразования; черная кривая соответствует вашим исходным данным, где значения v = NA были заменены оценками соответствия нелинейной модели (поэтому для v > 15 черные и синие кривые перекрываются).

person Maurits Evers    schedule 14.03.2018
comment
Спасибо за ваш ответ. Могу ли я сделать прогноз, чтобы заменить NA в w условием сверху, что наклон или тренд должны быть аналогичны тому, что у меня есть в моих справочных данных, p ? - person user9112767; 14.03.2018
comment
с условием сверху, что наклон или тенденция должны быть аналогичны тому, что у меня есть в моих справочных данных, p Не уверен, что вы имеете в виду. Как определить подобное? Я подгоняю кривую LOESS, которая пытается подогнать гладкую кривую между двумя переменными. Все дело в том, что вы не делаете никаких предположений относительно глобальной тенденции (отсюда и локальная в локальной полиномиальной регрессии). Я все еще не уверен, что вы пытаетесь сделать. Нельзя ли использовать явную функциональную форму t = f(t) для замены NA значений в w$t? Возможно, вам нужна модель нелинейной регрессии? - person Maurits Evers; 14.03.2018
comment
Если вы посмотрите на набор данных p, он выглядит как степенной закон. Я хотел бы заменить NA на w таким образом, чтобы в конце я получил плавный степенной закон. Используя LOESS, он возвращает параболу, а другой f <- function(v) (0.16 * (0.75 * v + 5.6)) / v; w$t[is.na(w$t)] <- f(w$v[is.na(w$t)]); имеет разрыв - person user9112767; 14.03.2018
comment
Да, именно поэтому я попросил вас уточнить. Итак, вы хотите подогнать к вашим данным нелинейную функцию, основанную на степенном законе, и использовать модель для замены значений NA. Пожалуйста, отредактируйте свой пост, чтобы включить эти детали. - person Maurits Evers; 14.03.2018
comment
@user9112767 user9112767 Я пересмотрел свое решение, посмотрите, пожалуйста. - person Maurits Evers; 14.03.2018
comment
Я получаю сообщение об ошибке: Error in a * v^b : non-numeric argument to binary operator In addition: Warning message: In nls(t ~ a * v^b, data = w[-1, ], na.action = na.exclude) : No starting values specified for some parameters. Initializing ‘b’ to '1.'. Consider specifying 'start' or using a selfStart model - person user9112767; 14.03.2018
comment
@user9112767 user9112767 Да, к сожалению, nls может быть несколько сложно. Проблема заключается в недостающих начальных значениях. Я отредактировал свое решение, включив в него набор начальных параметров для a и b для примера данных, который вы предоставили. Вам придется настроить их вручную для ваших данных; обычно это означает методом проб и ошибок. Добро пожаловать в мир подбора нелинейных моделей. - person Maurits Evers; 14.03.2018
comment
@ user9112767 Возможно, альтернативой было бы преобразование ваших данных в журнал; тогда у вас есть log y = log a + b log x, и вы можете использовать простую линейную модель формы lm(log(t) ~ log(v), data = w[-1, ]). - person Maurits Evers; 14.03.2018
comment
@user9112767 user9112767 Я добавил пример с использованием линейной модели после логарифмического преобразования ваших данных. Взгляни, пожалуйста. - person Maurits Evers; 14.03.2018
comment
Спасибо за очень полезные решения - person user9112767; 14.03.2018
comment
@user9112767 user9112767 Модель в mod не имеет никаких параметров (у вас нет a или b), поэтому на самом деле вы не подходите к (нелинейной) модели; отсюда и ошибка. Если вы хотите построить график функции, взгляните на ?stat_function. - person Maurits Evers; 15.03.2018