LS-регрессия с ограничением порядка

Я новичок в R и буду признателен за вашу помощь. У меня проблема оптимизации с ограничением. Хотя в R есть несколько способов решить задачу оптимизации, я не смог правильно выразить свою проблему с ограничением, которое мне нужно применить.

Предположим, у меня есть следующие данные в трех категориях:

A<-c(99.1,  96.5,   94.4,   92.7,   91.5,   91.3,   91.4,   90.1,   87.1,   82.6,   76.4)
B<-c(146.4, 140.2,  133.6,  126.5,  118.7,  109.4,  101.2,  101.8,  103.7,  102.5,  98.3)
C<-c(237.5, 213.9,  191,    168.9,  147.4,  124.9,  108.3,  95.7,   84.4,   73.5,   63)
t<-seq(1:11)
DT<-cbind.data.frame(t,A,B,C)

Я хотел бы подогнать экспоненциальную функцию y(t) к точкам данных в каждой из категорий (минимизируя квадраты ошибок), чтобы y(t)_c > y(t)_b > y(t)_a > 0 для выбранного t [1;15]


r optimization regression least-squares simplex
person Walle    schedule 17.10.2016    source источник

Ответы (2)


arrow_upward
0
arrow_downward

Я бы не стал пытаться включить ограничение в регрессию. Просто постройте три отдельные регрессии:

fit_loga <- lm(y ~ log(A) + t, data=DT) 
fit_logb <- lm(y ~ log(B) + t, data=DT) 
fit_logc <- lm(y ~ log(C) + t, data=DT)

fit_a <- exp(A)
fit_b <- exp(B)
fit_c <- exp(C)

затем убедитесь, что они удовлетворяют ограничению везде в этом диапазоне (или, по крайней мере, в каждой целочисленной точке данных): (fit_c > fit_b) & (fit_b > fit_a). Только если нет, то мы беспокоимся об этом. например, добавить некоторые другие термины в t в модель, возможно, exp(t), I(t^2), poly(t, <order>)...

РЕДАКТИРОВАТЬ: я не знал о пакете solnp.

person smci    schedule 20.10.2016
comment
smci спасибо за ответ! Я думал сделать это в два этапа, как вы предлагаете, но как бы вы это сделали в этом случае, когда C явно пересекает B и A? Я боюсь, что мне нужно сохранить одно и то же предположение о форме кривой (здесь: экспоненциальная) для всех кривых. Следовательно, параметры одного из них будут влиять на параметры другого, что требует выполнения оптимизации для них всех сразу. Я полностью согласен с тем, что лучше искать параметры линейных кривых в логарифмическом пространстве, чем пытаться подобрать экспоненциальные кривые. - person Walle; 21.10.2016
comment
О, дерьмо. Хорошо, если они пересекаются, то как мы применяем это ограничение? Не делая чего-то неуклюжего, например, зажимая, используя min/max? - person smci; 21.10.2016
comment
Я думал, что проблема может быть решена с помощью функции оптимизации, такой как {solnp()}, которая может обрабатывать ограничения неравенства. Однако он не сходится или приводит к недопустимому решению при непосредственном применении для минимизации квадратичных ошибок с заданными ограничениями. - person Walle; 23.10.2016
comment
@Walle: я не знаю solnp. Предлагаем вам добавить свой код и недопустимую ошибку к деталям и заголовку вопроса. Кроме того, это было бы более актуально на CrossValidated.SE. - person smci; 23.10.2016
comment
Оказалось, что это проблемы форматирования с solnp(), поскольку функция ожидала, что все параметры будут оценены как входные данные в одном единственном векторе. Поэтому мне пришлось преобразовать имеющуюся у меня матрицу в вектор, и это сработало. Не совсем понятно из сообщений об ошибках, которые появляются при запуске функции, т.е. non-numeric argument to binary operator для невыполнимых ограничений (слишком жестких). Не уверен, стоит ли мне публиковать код или просто закрыть эту тему. - person Walle; 25.10.2016
comment
@Walle: рад, что ты решил это. Пожалуйста, опубликуйте свой код и добавьте свой собственный ответ. Да, сообщения об ошибках R часто загадочны и плохо документированы. Если вы можете предложить авторам solnp, как улучшить это сообщение об ошибке и/или документацию, сделайте это во что бы то ни стало (как здесь, так и в docbug или в их списке), это единственный способ улучшить ситуацию. - person smci; 25.10.2016

arrow_upward
0
arrow_downward

Я понял, что сообщения об ошибках, которые появляются при использовании solnp, в основном относятся к неадекватным ограничениям. Также, как сказано в документации, требуется уложить все параметры в один вектор. После соответствующих корректировок кода я смог реализовать свои ограничения y(t)_c > y(t)_b > y(t)_a > 0 напрямую, без необходимости изменения проблемы. Удобнее всего использовать для решателя матричное формирование. Используя приведенные выше данные, я получил следующее: Результаты показаны здесь

library(data.table)
library(Rsolnp)

t<-as.vector(10:20)
DT<-cbind.data.frame(A,B,C)
tlogDT<-transpose(log(DT))

# min[log(y)'- Ax-b]
# Arr = [A1 A2 .. An b1 b2 .. bn]

gofn = function(arrin)
{
  arr = cbind(arrin[1:3],arrin[4:6])
  sum(
    (as.matrix(arr[,1])%*%t + arr[,2] - tlogDT)^2
  )
}

nocross=t #defines the times where the curves are not allowed to intersect
ineqfn2=function(arrin)
{
  #constrains:
  # 0<f_a(t)<f_b(t)<f_c(t), for some t,
  arr = cbind(arrin[1:3],arrin[4:6])
  nextarr=as.matrix(rbind(rep(0,2),arr[1:(length(arr[,1])-1),]))
  ineqmat=as.matrix(arr[,1])%*%nocross+arr[,2]-nextarr[,1]%*%nocross-nextarr[,2]
  as.vector(t(ineqmat))
}

#lines should be aligned with the first startvalue
eqfn = function(arrin)
{  
  arr = cbind(arrin[1:3],arrin[4:6])
  arr[,1]*t[1]+arr[,2]-tlogDT[,1]
}
#start values:
An=c(1,1,1)
bn=tlogDT[,1]
vstart=c(An,bn)

r <- solnp(
  vstart,  gofn,

  eqfun = eqfn, eqB= c(0,0,0),

  ineqfun = ineqfn2,
  ineqLB = rep(0,length(DT[1,])*length(nocross)),
  ineqUB = rep(5000,length(DT[1,])*length(nocross))
)

r$pars[1]
line1 = exp(r$pars[4]+r$pars[1]*t)
line2 = exp(r$pars[5]+r$pars[2]*t)
line3 = exp(r$pars[6]+r$pars[3]*t)

plot(t, DT[,3],log = "y")
points(t, DT[,2],col="green")
points(t, DT[,1],col="blue")

lines(t, line1,lwd=2, col = "blue",  xlab = "Time (s)", ylab = "Counts")
lines(t, line2,lwd=2, col = "green",  xlab = "Time (s)", ylab = "Counts")
lines(t, line3,lwd=2, col = "black",  xlab = "Time (s)", ylab = "Counts")
person Walle    schedule 26.10.2016