Это будет длинная история, возможно, больше подходит для https://stats.stackexchange.com/.
====== Часть 1 - Проблема ======
Это последовательность генерации ошибки:
library(fExtremes)
samp <- read.csv("optimdata.csv")[ ,2]
## does not converge
para <- gevFit(samp, type = "mle")
Мы сталкиваемся с типичной причиной lack- конвергенции при использовании optim() и друзей: неадекватные стартовые значения для оптимизации.
Чтобы узнать, что пойдет не так, дадим оценку PWM (http://arxiv.org/abs/1310.3222); это состоит из аналитической формулы, следовательно, она не несет в проблемы конвергенции, так как это не имеет никакого использования optim():
para <- gevFit(samp, type = "pwm")
fitpwm<- attr(para, "fit")
fitpwm$par.ests
Предполагаемый параметр хвоста xi является отрицательным, что соответствует ограниченному верхнему хвосту; на самом деле встроенные распределительные дисплеи еще более «верхний хвост ограниченности», чем данные выборки, как вы можете увидеть от «выравнивания» из квантиль-квантиль графике справа:
qqgevplot <- function(samp, params){
probs <- seq(0.1,0.99,by=0.01)
qqempir <- quantile(samp, probs)
qqtheor <- qgev(probs, xi=params["xi"], mu=params["mu"], beta=params["beta"])
rang <- range(qqempir,qqtheor)
plot(qqempir, qqtheor, xlim=rang, ylim=rang,
xlab="empirical", ylab="theoretical",
main="Quantile-quantile plot")
abline(a=0,b=1, col=2)
}
qqgevplot(samp, fitpwm$par.ests)
Для xi<0.5 ОМП оценивани не является регулярным (http://arxiv.org/abs/1301.5611): значение -0,46, оцененное PWM для xi, очень близко к этому. Теперь оценки ШИМ используются внутри gevFit() в качестве исходного значения для optim(): вы можете увидеть это, если вы распечатать код для функции gevFit():
print(gevFit)
print(.gevFit)
print(.gevmleFit)
Отправной значение для Оптим является theta, полученный PWM. Для конкретных данных это исходное значение не является адекватным, поскольку оно приводит к неконвергенции optim().
====== Часть 2 - решения? ======
Решение 1 заключается в использовании para <- gevFit(samp, type = "pwm"), как указано выше. Если вы хотите использовать ML, тогда вы должны указать хорошие начальные значения для optim(). К сожалению, пакет fExtremes не позволяет это сделать. Затем вы можете переопределить свою собственную версию .gevmleFit, чтобы включить ее, например.
.gevmleFit <- function (data, block = NA, start.param, ...)
{
data = as.numeric(data)
n = length(data)
if(missing(start.param)){
theta = .gevpwmFit(data)$par.ests
}else{
theta = start.param
}
fit = optim(theta, .gevLLH, hessian = TRUE, ..., tmp = data)
if (fit$convergence)
warning("optimization may not have succeeded")
par.ests = fit$par
varcov = solve(fit$hessian)
par.ses = sqrt(diag(varcov))
ans = list(n = n, data = data, par.ests = par.ests, par.ses = par.ses,
varcov = varcov, converged = fit$convergence, nllh.final = fit$value)
class(ans) = "gev"
ans
}
## diverges, just as above
.gevmleFit(samp)
## diverges, just as above
startp <- fitpwm$par.ests
.gevmleFit(samp, start.param=startp)
## converges
startp <- structure(c(-0.1, 1, 1), names=names(fitpwm$par.ests))
.gevmleFit(samp, start.param=startp)$par.ests
Теперь проверьте это: beta оценивается ШИМ 0.1245; изменив это крошечное количество, ОМП производится сходиться:
startp <- fitpwm$par.ests
startp["beta"]
startp["beta"] <- 0.13
.gevmleFit(samp, start.param=startp)$par.ests
Это, надеюсь, ясно показывает, что слепо optim() ISE работает до тех пор, пока не могут и затем превратиться в весьма деликатное начинание действительно. По этой причине было бы полезно оставить этот ответ здесь, а не перейти на CrossValidated.
Трудно сказать, что пошло не так, потому что ваш пример не воспроизводится: см. Здесь. http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example Вероятность отказа связана с особенностями вашего набора данных 'sample' (например, повторных значений); Я не могу сделать много, если вы не публикуете хотя бы подмножество данных или, может быть, другой набор данных, где возникает проблема. –
Привет, я обновил свои данные здесь! –