(при необходимости игнорируйте код R, так как мой главный вопрос не зависит от языка)
Если я хочу посмотреть на изменчивость простой статистики (например, среднее), я знаю, что могу сделать это с помощью теории, например:
x = rnorm(50)
# Estimate standard error from theory
summary(lm(x~1))
# same as...
sd(x) / sqrt(length(x))
или с помощью начальной загрузки, как:
library(boot)
# Estimate standard error from bootstrap
(x.bs = boot(x, function(x, inds) mean(x[inds]), 1000))
# which is simply the standard *deviation* of the bootstrap distribution...
sd(x.bs$t)
Однако, что мне интересно, это может быть полезно / правильно (?) Смотреть на стандартную ошибку дистрибутива начальной загрузки в определенных ситуациях? Ситуация, с которой я имею дело, является относительно шумной нелинейной функцией, такой как:
# Simulate dataset
set.seed(12345)
n = 100
x = runif(n, 0, 20)
y = SSasymp(x, 5, 1, -1) + rnorm(n, sd=2)
dat = data.frame(x, y)
Здесь модель даже не сходится, используя исходный набор данных,
> (fit = nls(y ~ SSasymp(x, Asym, R0, lrc), dat))
Error in numericDeriv(form[[3L]], names(ind), env) :
Missing value or an infinity produced when evaluating the model
поэтому интересующая меня статистика - это более устойчивые оценки этих параметров nls - возможно, их значения для ряда загрузочных репликаций.
# Obtain mean bootstrap nls parameter estimates
fit.bs = boot(dat, function(dat, inds)
tryCatch(coef(nls(y ~ SSasymp(x, Asym, R0, lrc), dat[inds, ])),
error=function(e) c(NA, NA, NA)), 100)
pars = colMeans(fit.bs$t, na.rm=T)
Вот они, действительно, в центре событий того, что я использовал для имитации исходных данных:
> pars
[1] 5.606190 1.859591 -1.390816
Верстка выглядит так:
# Plot
with(dat, plot(x, y))
newx = seq(min(x), max(x), len=100)
lines(newx, SSasymp(newx, pars[1], pars[2], pars[3]))
lines(newx, SSasymp(newx, 5, 1, -1), col='red')
legend('bottomright', c('Actual', 'Predicted'), bty='n', lty=1, col=2:1)
Теперь, если мне нужна изменчивость этих оценок стабилизированных параметров, я думаю, что при условии нормальности этого распределения при начальной загрузке я могу просто рассчитать их стандартные ошибки:
> apply(fit.bs$t, 2, function(x) sd(x, na.rm=T) / sqrt(length(na.omit(x))))
[1] 0.08369921 0.17230957 0.08386824
Это разумный подход? Существует ли лучший общий подход к выводу параметров нестабильных нелинейных моделей, подобных этому? (Я полагаю, что вместо этого я мог бы сделать второй слой передискретизации здесь, вместо того, чтобы в последний раз полагаться на теорию, но это может занять много времени в зависимости от модели. Даже если я не уверен, будут ли эти стандартные ошибки быть полезным для всего, так как они приблизятся к 0, если я просто увеличу количество загрузочных репликаций.)
Большое спасибо, и, кстати, я инженер, поэтому, пожалуйста, прости меня за то, что я новичок здесь.
источник
nls
подгонок может потерпеть неудачу, но из тех, которые сходятся, смещение будет огромным, а прогнозируемые стандартные ошибки / КИ крайне малы.nlsBoot
использует специальное требование 50% успешных совпадений, но я согласен с вами, что (не) сходство условных распределений также вызывает озабоченность.