интерполяция между двумя кривыми в R

У меня есть две кривые с одинаковым количеством точек, заполненные с использованием приблизительной функции, для значений x и y отдельно для каждой кривой. Значения осей x и y являются логарифмическими, поэтому при аппроксимации и интерполяции я возвращаюсь к обычной десятичной шкале. Черная и синяя линии являются исходными линиями, а красная интерполирована между ними. Как видите, красная линия не имитирует изгиб с правой стороны, так как интерполяция выполняется на основе предположения, что каждая пара x и y является ближайшей.

Есть ли способ выполнить интерполяцию между кривыми в R на основе реальных ближайших точек между ними? Может быть, существуют алгоритмы для этого? Что угодно было бы полезно, так как я не уверен, как это называется в математике.

    base="ftp://cdsarc.u-strasbg.fr/pub/cats/J/A+A/508/355/ms/"
    setwd("~/Desktop")
    file1=paste(base,"z001y23_1.60.dat",sep="")
    file2=paste(base,"z001y23_1.70.dat",sep="")

    cols=c("no","age","logL","logTef", "grav","stage")
    ncol <- length(count.fields(file=file1, sep = ","))
    second=read.table(file=file1,fill=T, blank.lines.skip=F, skip=2, header=F, strip.white=T, col.names = paste("V", seq_len(ncol)))
    second$V.6<-second$V.23
    colnames(second) <-cols
    second$logL=as.numeric(second$logL)
    #performing some filtering of data here
    pos1=which(second$stage == "trgb")[1]
    second=second[1:pos1,]

    ncol <- length(count.fields(file=file2, sep = ","))
    first=read.table(file=file2,fill=T, blank.lines.skip=F, skip=2, header=F, strip.white=T, col.names = paste("V", seq_len(ncol)))
    first$V.6<-first$V.23
    colnames(first) <-cols
    #performing some filtering of data here
    pos2=which(first$stage == "trgb")[1]
    first=first[1:pos2,]

    #plotting data
    len=max(c(min(first[[4]]),min(second[[4]])))
    first=first[first[[4]]>len,]
    second=second[second[[4]]>len,]

    plot(second[[4]],second[[3]],t="l",xlim=rev(range(second[[4]])),xlab="x",ylab="y")
    lines(first[[4]],first[[3]],t="l",col="blue")
    n=max(c(length(second[[4]]),length(first[[4]])))
    #approximating missing points
    xf1 <- approx(10^second[[4]],n=n)
    yf1 <- approx(10^second[[3]],n=n)

    xf2 <- approx(10^first[[4]],n=n)
    yf2 <- approx(10^first[[3]],n=n)

    #calculating interpolated line
    ratio=2
    s1<-log10((xf1$y-xf2$y)/ratio+xf2$y)
    s2<-log10((yf1$y-yf2$y)/ratio+yf2$y)
    lines(s1,s2, col ="red")

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


r interpolation approximation curves
person saganas    schedule 12.06.2014    source источник

Ответы (1)


arrow_upward
3
arrow_downward

Хотя это и не окончательный ответ, вот кое-что, адаптированное из того, что я сделал некоторое время назад для миграции потокового канала. Обратите внимание, что они обычно не пересекаются самостоятельно, поэтому ваш пробег может отличаться. Вся идея состоит в том, чтобы рассчитать кривизну и использовать динамическую деформацию времени. соответствовать экстремумам.

Грубо можно резюмировать так:

  1. Параметризуйте обе кривые так, чтобы L1 и L2 были векторами, представляющими длины от начала кривой до рассматриваемого индекса.
  2. Рассчитайте smooth.spline xsp1, ysp1, xsp2, ysp2, используя L1 и L2 для x и y каждой кривой. Обратите внимание на параметр сглаживания, так как ваши кривые временами выглядят резкими.
  3. Явно получить кривизна со знаком для каждой сглаженной линии
  4. Используйте dtw для сопоставления пиков кривизны каждой сглаженной линии.
  5. Используйте индексы, возвращенные dtw, чтобы установить сопоставление между кривыми
  6. ...
  7. ВЫГОДА!!!

Учтите, что dtw не творит чудес, и для этого потребуются некоторые эксперименты.

P.S. Чтобы сэкономить ваше время, я попытался использовать dtw непосредственно для x и y без кривизны, но это не очень хорошо, так как мы хотели бы отображать обе координаты одновременно.

ИЗМЕНИТЬ

library(dtw)
df1 <- data.frame(x=first[[4]], y=first[[3]])
df2 <- data.frame(x=second[[4]], y=second[[3]])
measure <- function(df)
  within(df, m <- c(0, cumsum(diff(x)^2 + diff(y)^2)))
df1 <- measure(df1)
df2 <- measure(df2)

curvify <- function(df) {
  xsp <- with(df, smooth.spline(m, x))
  ysp <- with(df, smooth.spline(m, y))
  xx <- predict(xsp, df$m)$y
  yy <- predict(ysp, df$m)$y
  xp <- predict(xsp, df$m, deriv=1)$y
  xpp <- predict(xsp, df$m, deriv=2)$y
  yp <- predict(ysp, df$m, deriv=1)$y
  ypp <- predict(ysp, df$m, deriv=2)$y
  # http://en.wikipedia.org/wiki/Curvature#Signed_curvature
  within(df, c <- (xp*ypp - yp*xpp)/(xp^2 + yp^2)^1.5)
}

df1 <- curvify(df1)
df2 <- curvify(df2)

d <- dtw(df1$c, df2$c, keep=TRUE)
# plot(d, type='three')

xx <- ( df1$x[d$index1] + df2$x[d$index2] ) /2
yy <- ( df1$y[d$index1] + df2$y[d$index2] ) /2

lines(xx, yy, col="green")

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

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

РЕДАКТИРОВАТЬ

Для интерполяции с весами, отличными от 1/2

fr <- 1/3
xx <- df1$x[d$index1] * fr + df2$x[d$index2] * (1-fr)
yy <- df1$y[d$index1] * fr + df2$y[d$index2] * (1-fr)
lines(xx, yy, col="yellow")

fr <- 2/3
xx <- df1$x[d$index1] * fr + df2$x[d$index2] * (1-fr)
yy <- df1$y[d$index1] * fr + df2$y[d$index2] * (1-fr)
lines(xx, yy, col="brown")
person mlt    schedule 13.06.2014
comment
Я только что понял, что забыл sqrt для евклидова расстояния, но, тем не менее, результат выглядит хорошо для меня и становится уродливым с этим... да ладно... лень смотреть на это. Возможно, это влияет на сглаживание дефолтов. - person mlt; 13.06.2014
comment
Большое тебе спасибо! Выглядит именно то, что мне было нужно. - person saganas; 16.06.2014
comment
Вы можете принять это как ответ, если он работает для вас. Также обратите внимание, что он соответствует кривизне, а не ближайшим точкам, поэтому кривые должны иметь одинаковые повороты. - person mlt; 16.06.2014
comment
Это именно то, что мне было нужно. Так что это достаточно хорошо, просто отрегулировал вес на 1,5, чтобы уменьшить его, и кажется, что все хорошо. - person saganas; 19.06.2014
comment
Большое спасибо за интерполяцию произвольной дроби. - person saganas; 19.06.2014