Глобальная оптимизация по сайтам для отрицательной экспоненциальной функции - PullRequest
0 голосов
/ 10 мая 2018

Попытка оценить лучшие параметры, которые подходят для всего исследования. Каждое исследование имеет различную длину наборов данных и должно видеть, каковы наилучшие параметры для отрицательной экспоненциальной модели, согласованной между переменными «Остальное» и «Время». Примерный набор данных выглядит следующим образом;

Study <- as.factor(c(1, 1, 1, 1, 2, 2, 2, 3, 3))
Time <- as.numeric(c(0, 0.08, 0.16, 0.24, 0, 0.05, 0.88, 0, 0.99))
Remaining <- as.numeric(c(100, 80, 69, 45, 100, 60, 35, 0, 25))
data_n <- cbind(Study, Time, Remaining)
head(data_n)

Необходимость оптимизации отрицательной экспоненциальной функции позволяет рассчитывать K и K2, которые подходят для всех сайтов. Функция

Predicted_remaining <- 41* exp(-k*Time) + (100-41) * exp(-k2*Time)

Значение 41 является константой, а -k и -k2 - параметры отрицательной экспоненциальной функции, которые необходимо оптимизировать. Попытка оптимизировать в R, так что любой совет с благодарностью.

1 Ответ

0 голосов
/ 10 мая 2018

Сначала я создаю фрейм данных из ваших данных. Затем я определяю функцию, которая будет соответствовать данным, которая принимает параметры подбора, p и time в качестве аргументов. Затем я создаю функцию стоимости, которая вычисляет сумму квадратов как показатель качества соответствия. После этого я использую optim, чтобы минимизировать функцию стоимости и, наконец, отобразить результаты.

# Create data frame
data_n <- data.frame(Study, Time, Remaining)

# Function to be fit
fit_function <- function(p, time){
  (41 * exp(-p[1] * time) + 59 * exp(-p[2] * time))
}

# Cost function using the sum of squares
cost_function <- function(p, data){
  sum((data$Remaining - fit_function(p, data$Time))^2)
}

# Using 'optim' to minimise the cost function
fit <- optim(c(1, 1), cost_function, data = data_n)

# Plot results
plot(data_n$Time, data_n$Remaining, xlab = "Time", ylab = "Remaining")
lines(seq(0, 1, by = 0.1), fit_function(fit$par, seq(0, 1, by = 0.1)))

enter image description here

...