Это было интересно, и я не совсем доволен открытием, но я расскажу вам, что я нашел, и посмотрю, поможет ли это.
При вызове функции fitdist по умолчанию она хочет использовать mledist из того же пакета. Это само по себе приводит к вызову stats::optim, который является общей функцией оптимизации. В возвращаемом значении он дает код ошибки сходимости, подробности см. в разделе ?optim. 100, которое вы видите, не является одним из тех, что возвращает optim. Поэтому я разобрал код для mledist и fitdist, чтобы найти, откуда берется этот код ошибки. К сожалению, он определен более чем в одном случае и является общим кодом ошибки ловушки. Если вы разберете весь код, то fitdist попытается сделать следующее, с учетом различных предварительных проверок и т. д.
fnobj <- function(par, fix.arg, obs, ddistnam) {
-sum(do.call(ddistnam, c(list(obs), as.list(par),
as.list(fix.arg), log = TRUE)))
}
vstart = list(xmin=5,alpha=5)
fnobj <- function(par, fix.arg obs, ddistnam) {
-sum(do.call(ddistnam, c(list(obs), as.list(par),
as.list(fix.arg), log = TRUE)))
}
ddistname=dplcon
fix.arg = NULL
meth = "Nelder-Mead"
lower = -Inf
upper = Inf
optim(par = vstart, fn = fnobj,
fix.arg = fix.arg, obs = data, ddistnam = ddistname,
hessian = TRUE, method = meth, lower = lower,
upper = upper)
Если мы запустим этот код, мы обнаружим более полезную ошибку «функция не может быть оценена при начальных параметрах». Что имеет смысл, если мы посмотрим на определение функции. Наличие xmin=0 или alpha=1 даст логарифмическую вероятность -Inf. Итак, подумайте, попробуйте разные начальные значения, я попробовал несколько случайных вариантов, но все они вернули новую ошибку: «неконечное конечно-разностное значение 1".
При дальнейшем поиске источника этих двух ошибок в источнике optim они не являются частью самого источника R, однако есть вызов .External2, поэтому я могу только предположить, что ошибки исходят оттуда. Неконечная ошибка подразумевает, что одно из вычислений функции где-то дает нечисловой результат. Функция dplcon сделает это, когда alpha <= 1 или xmin <= 0. fitdist позволяет вам указать дополнительные аргументы, которые передаются mledist или другому (в зависимости от выбранного вами метода, по умолчанию используется mle), из которых lower — один для управления нижними границами оптимизируемых параметров. Поэтому я попытался ввести эти ограничения и повторить попытку:
fitpl <- fitdist(data,"plcon",start = list(xmin=1,alpha=2), lower = c(xmin = 0, alpha = 1))
Досадно, что это по-прежнему дает код ошибки 100. Отслеживание этого приводит к ошибке «L-BFGS-B нужны конечные значения 'fn'». Метод оптимизации Нелдера-Мида по умолчанию изменился по мере того, как вы указываете границу, и где-то во внешнем вызове кода C возникает эта ошибка, предположительно близко к пределам либо xmin, либо alpha, где стабильность численного расчета по мере приближения к бесконечности равна важный.
Я решил сделать сопоставление квантилей, а не максимальную вероятность, чтобы попытаться узнать больше.
fitpl <- fitdist(data,"plcon",start = list(xmin=1,alpha=2),
method= "qme",probs = c(1/3,2/3))
fitpl
## Fitting of the distribution ' plcon ' by matching quantiles
## Parameters:
## estimate
## xmin 0.02135157
## alpha 46.65914353
что говорит о том, что оптимальное значение xmin близко к 0, это пределы. Причина, по которой я не удовлетворен, заключается в том, что я не могу получить соответствие распределения с максимальным правдоподобием, используя fitdist, однако, надеюсь, это объяснение поможет, а сопоставление квантилей даст альтернативу.
Изменить:
Узнав немного больше о распределении по степенному закону в целом, становится понятно, что это не работает так, как вы ожидаете. Параметр мощности параметра имеет функцию правдоподобия, которая может быть максимизирована при заданном xmin. Однако такого выражения для xmin не существует, поскольку функция правдоподобия возрастает по xmin. Обычно оценка xmin берется из статистики Колмогорова--Смирнова, см. этот вопрос mathoverflow и виньетка d_jss_paper пакета poweRlaw для получения дополнительной информации и соответствующих ссылок.
В самом пакете poweRlaw есть функционал для оценки параметров степенного распределения.
m = conpl$new(data)
xminhat = estimate_xmin(m)$xmin
m$setXmin(xminhat)
alphahat = estimate_pars(m)$pars
c(xmin = xminhat, alpha = alphahat)
person
jamieRowen
schedule
11.05.2016
