Блог

Проще говоря, у меня есть функция f(x, t2) , и я хочу найти значение x , которое максимизирует интеграл f(x, t2) по отношению к t2 . Я выбираю алгоритм pso для оптимизации. Извиняемый код выглядит следующим образом

 library(pso)
xl=0; xu=2000; n=1; t2l=100; t2u=2000; t1=1
g<-function(x, t2) t1*x/(t2 x)
h<-function(z) 1/z^n
gdot<-function(x, t2){
  c(x/(t2 x),-t1*x/(t2 x)^2)
}
logdetHinv<-function(dp, dw, t2){
  gmat=mapply(function(x) gdot(x,t2),dp)
  D0=gmat%*%diag(dw)%*%t(gmat)
  D1=gmat%*%diag(1/h(g(dp,t2)))%*%diag(dw)%*%t(gmat)
  2*log(det(D1))-log(det(D0))
}
obj<-function(x){
  dp=x[1:2]; dw=c(x[3],1-x[3])
  fitness_value=-integrate(Vectorize(function(t2) logdetHinv(dp, dw, t2)*1/(t2u-t2l)), t2l, t2u)$value
  return(ifelse(dw[2]>0, fitness_value, fitness_value 1e3))
}

x <- psoptim(rep(1,3), fn = obj, lower = c(rep(xl,2),0.1), upper = c(rep(xu,2), 0.9))$par
x
 

Поскольку глобальная оптимизация включает в себя некоторую случайную процедуру, она иногда выдает правильный результат

 > x
[1] 2000.0000  754.4146    0.5000
 

в других случаях он сообщает об ошибке

 Error in integrate(Vectorize(function(t2) logdetHinv(dp, dw, t2) * 1/(t2u -  : 
  non-finite function value
In addition: There were 11 warnings (use warnings() to see them)
> warnings()
Warning messages:
1: In log(det(D1)) : NaNs produced
2: In log(det(D0)) : NaNs produced
3: In log(det(D1)) : NaNs produced
4: In log(det(D0)) : NaNs produced
 

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

Я хочу избежать таких значений, возможно , с tryCatch помощью , например, если в функции есть предупреждение logdetHinv , она возвращает очень маленькое значение , но нет NaN , поэтому это не приведет к ошибке integrate , и psoptim маловероятно, что при максимизации целевой функции (минимизации -integrate(logdetHinv) ) будут выбраны такие значения . Я не знаком с tryCatch такой сложной ситуацией. Куда мне его положить tryCatch ? Спасибо.

Кроме того, я хотел бы знать, существуют ли в R какие-либо методы отладки, которые позволяют мне знать, какое случайное значение ( D0/D1 ) вызывает ошибку в этом случае. Я предполагаю, что это какое-то отрицательное значение в журнале, но оно не должно быть, так как внутри журнала находится определитель положительно определенной матрицы. В режиме обратной трассировки, в разделе обзор, если я введу D0 объект «D0», он не будет найден.

В этом случае я бы не стал использовать tryCatch, который обычно более подходит для тестирования, чем в вашем основном коде. Почему бы вам просто не проверить детерминанты в вашей функции ? Что-то в этом роде должно сработать:

 logdetHinv<-function(dp, dw, t2){
    gmat=mapply(function(x) gdot(x,t2),dp)
    D0=gmat%*%diag(dw)%*%t(gmat)
    D1=gmat%*%diag(1/h(g(dp,t2)))%*%diag(dw)%*%t(gmat)
    detD1 <- max(0.01, det(D1))
    detD0 <- max(0.01, det(D0)) 
    2*log(detD1)-log(detD0)
}