Я бы попросил немного помочь с моей проблемой в поиске процедуры оптимизации.Я должен объяснить некоторые этапы своей работы и надеюсь, что смогу четко объяснить свою проблему ... В настоящее время я провожу исследование отсутствующих данных в R, работая с имитированными данными и реальными данными.В смоделированных данных я истолковал отсутствующие шаблоны данных, аналогичные реальным данным некоторыми методами обучения.То, что я сейчас хочу сделать, это некоторые уроки от моделирования до реальных данных.
В моделировании у меня есть вектор истинных параметров с параметрами регрессионной модели, оцененными на основе полных данных без пропусков.Назовем этот вектор REF.При искусственном генерировании пропусков и множественного вменения (например, с MICE) пропущенных значений я, наконец, получаю параметры регрессионной модели на основе пропущенных данных.Назовем эти параметры регрессии MICE.Теперь я хочу смоделировать разницу между REF и MICE, давайте назовем их PA-REF и перенесем это обучение в реальные данные.Довольно хорошее моделирование для PA-REF возможно с экстремальным распределением значений.Это означает, что обычно MICE находится в основном рядом с REF, но в некоторых случаях нет.Таким образом, для построения распределения для оценочных параметров регрессионной модели в +++ реальных данных (!) С отсутствием и подключением вменений через Mice +++ я получил очень широкое распределение с экстремальным значением вокруг параметров регрессии.Для этих обстоятельств я оцениваю логит с некоторыми ковариатами, получающими вероятности, если параметр регрессии является выбросом или нет (с точки зрения разницы REF и MICE в моделировании).Я пытаюсь управлять EVD-параметрами с этими вероятностями, используя прибл.
В приведенном ниже коде вы найдете реализацию этой идеи.
В начале имеет место базовая функция зависимости evd-параметра от вероятности выброса: просто линейная функция с начальным и конечным значением.Моя отправная точка - это различие между низкими и высокими значениями PA.Ref.Для обеих групп я включил выброс.Для низких значений отсечка для выбросов (10), для высоких значений самое низкое значение из группы низких значений.Я делаю это, чтобы включить неуверенность в том, что оценочные вероятности отражают неопределенность оценки параметров регрессии, основанной на MICE, как выброса или нет.В симуляции модель для предсказания выброса параметра (высокий PA-Ref) работает очень хорошо, но есть некоторые недостатки, которые здесь не нужно объяснять.
На следующем шаге я использую этофункция для извлечения значений из EVD, управляемой таким образом, то есть 125 значений для PA-REF с простой линейной функцией для параметров evd в зависимости от вероятностей выбросов.Тиражи должны соответствовать требованию быть выше 0 и не выше 1,5 макс. PA-REF.Управляемые отрисовки evd дают распределение для каждого параметра регрессии.Здесь я получаю матричный формат 125x44, строки которого могут быть удвоены из-за неизвестного направления PA-REF.Я снова приближаюсь к этому evd-распределению и могу вычислить значение pdf для наблюдаемого PA-REF (разность REF и MICE) в моделировании, см. Следующий шаг.
Третий шаг в основном делает это для каждого элемента y.pa.ref: Приблизьте плотность отрисовок evd и интегрируйте эту приближенную функцию с наблюдаемым y.pa.ref, выполняя при этом некоторые требования (начало.val, end.val, max.val) для получения вероятностной массы покрытия для наблюдаемого PA-REF.Наконец, среднее - это показатель того, насколько хорошо evd покрывает значения PA-Ref в среднем.Здесь он равен 0,33, что означает, что при таком подходе это средняя вероятность того, что y.pa.ref будут охвачены в extr.Распределение значений указывается для каждого параметра регрессионной модели.(Если учесть направление, значение может быть удвоено.)
То, что я хочу оптимизировать, заключается в следующем: теперь я хочу максимизировать среднее значение pdf-значений для всех элементов в наблюдаемом векторе PA-REF.Итак, мой вопрос: возможно ли оптимизировать эти значения, указав оптимизированную функцию для параметров EVD?
Я попробовал несколько идей самостоятельно.Например, я использовал пары с самыми высокими значениями prob.outl и y.pa.ref, извлекал значения из многовариантной нормали и использовал эти значения для более высокой оценки evd-параметров.Затем я строю функцию с линейным увеличением параметров evd в outl.prob от 0 до 0,5 и с постоянными параметрами evd более высоких параметров EVD на основе MVN с более чем 0,5 до 1. К сожалению, это не улучшает среднее значение массы pdf привсе.В основном я считаю, что это созвездие можно решить с помощью оптимизации.Но может я ошибаюсь ...?: -))
THX за вашу поддержку.Если вы хотите, я буду ссылаться на вашу помощь в моей публикации, по крайней мере, я должен анонимно!
PA.REF <- c(0.43997, 1.99193, 2.36845, 1.37931, 24.32432, 0.42373, 2.77907,
2.49998, 18.18181, 0.52219, 0.37312, 5.24694, 4.21052, 0.59523,
5.02122, 4.78839, 3.11992, 4.42476, 16.09194, 3.48839, 3.24674,
0.39386, 3.74449, 6.20155, 2.02429, 2.02953, 2.77776, 2.22222,
7.34693, 5.15163, 1.57566, 0.16821, 3.20001, 1.17646, 0.83331,
0.55556, 1.88996, 0.38461, 4.76189, 0.91186, 1.0118, 0.29985,
1.34529, 13.58024)
Outl.prob <- c(0.15371, 0.02757, 0.01265, 0.05866, 0.89814, 0.01828, 0.00442,
0.03018, 0.63696, 0.00845, 0.01356, 0.00507, 0.07715, 0.07176,
0.02763, 0.00671, 0.00998, 0.03585, 0.01629, 0.16794, 0.09089,
0.00406, 0.01003, 0.03254, 0.01134, 0.00711, 0.00424, 0.00589,
0.09664, 0.01311, 0.00851, 0.01428, 0.01049, 0.25723, 0.0257,
0.06623, 0.00728, 0.00789, 1, 0.00587, 0.01292, 0.00889, 0.0179,
1)
my.data <- data.frame(Outl.prob, PA.REF)
y.pa.ref <- my.data[,c("PA.REF")]
outl.prob <- my.data[,c("Outl.prob")]
cut.outl <- 10
require(extRemes)
###########################
#### Step 1
fit.evd.min <- fevd( c( sort(y.pa.ref)[1:(length(y.pa.ref)/2)],
cut.outl))
loc.evd.fit.min = fit.evd.min$results$par[c("location")]
scale.evd.fit.min = fit.evd.min$results$par[c("scale")]
shape.evd.fit.min = fit.evd.min$results$par[c("shape")]
fit.evd.max <- fevd( c( sort(y.pa.ref)[(length(y.pa.ref)/2): length(y.pa.ref)],
min(y.pa.ref)))
loc.evd.fit.max = fit.evd.max$results$par[c("location")]
scale.evd.fit.max = fit.evd.max$results$par[c("scale")]
shape.evd.fit.max = fit.evd.max$results$par[c("shape")]
basic.loc.evd.fun <- approxfun(approx(x= c(0,1),
y= c(loc.evd.fit.min, loc.evd.fit.max), n=1e4))
basic.scale.evd.fun <- approxfun(approx(x= c(0,1),
y=c(scale.evd.fit.min, scale.evd.fit.max), n=1e4))
basic.shape.evd.fun <- approxfun(approx(x= c(0,1),
y=c(shape.evd.fit.min, shape.evd.fit.max), n=1e4))
###########################
###########################
#### Step 2
number125 <- 1:125
evd.val.max <- max(y.pa.ref)*1.5
count.draw.simu <- 0
evd.var.vert.simu <- NULL
matrix.evd.v4 <- matrix(nrow=length(number125),
ncol=length(y.pa.ref))
set.seed(71108)
for (i in c(1:length(y.pa.ref))){
while (count.draw.simu < max(number125)){
evd.simu.draw <- rgev(1, loc=basic.loc.evd.fun(outl.prob[i]),
scale= basic.scale.evd.fun(outl.prob[i]),
shape=basic.shape.evd.fun(outl.prob[i]))
evd.simu.draw <- ifelse(evd.simu.draw<evd.val.max
& evd.simu.draw >0, evd.simu.draw, NA)
evd.var.vert.simu <- c(evd.var.vert.simu, evd.simu.draw)
evd.var.vert.simu <- evd.var.vert.simu[!is.na(evd.var.vert.simu)]
count.draw.simu <- length(evd.var.vert.simu)
}
matrix.evd.v4[,i] <- evd.var.vert.simu
evd.var.vert.simu <- NULL
count.draw.simu<- 0
}
EVD2 <- rbind(matrix.evd.v4, - matrix.evd.v4)
###########################
###########################
#### Step 3
fun.list.EVD2 <- list()
est.val.fun.EVD2 <- list()
est.prob.fun.EVD2 <- NULL
for (i in 1:length(y.pa.ref)){
EVD2.val <- EVD2[,i]
assign(paste0(i, ".fun.EVD2" ), approxfun( density(EVD2.val) ))
fun.list.EVD2[[i]] <- get(paste0(i,".fun.EVD2" ))
start.val <- min( density(EVD2.val)$x)
end.val <- y.pa.ref[i]
max.val <- max( density(EVD2.val)$x)
ind.fun.ex <- c("Yes")
ifelse( start.val < end.val & max.val > end.val,
ifelse(start.val > 0,
prob.EVD2fun <- integrate( get(paste0(i,".fun.EVD2" )),
start.val,
end.val) ,
prob.EVD2fun <- integrate( get(paste0(i,".fun.EVD2" )),
0,
end.val) ),
ind.fun.ex <- c("NO") )
if(ind.fun.ex == c("NO")){
est.val.fun.EVD2[[i]] <- NULL
est.prob.fun.EVD2 <- c(est.prob.fun.EVD2, 0)
} else {
assign(paste0("value.EVD2.fun", i), prob.EVD2fun)
est.val.fun.EVD2[[i]] <- get(paste0("value.EVD2.fun", i))
est.prob.fun.EVD2 <- c(est.prob.fun.EVD2, prob.EVD2fun$value)
}
}
###########################
optim.val.prob.EVD2 <- mean(est.prob.fun.EVD2)