Кодирование интервала предсказания из обобщенной модели присадок с очень большим набором данных

Question

У меня есть небольшой набор данных людей, их местонахождения и знают ли они друг друга. Это подмножество набора данных с 1000 человек. Учитывая, что каждый человек может знать любого другого человека, число потенциальных связей растет только под n^2. Я хочу подгонять модель с небольшим подмножеством, чтобы получить вероятность привязки в зависимости от расстояния, чтобы я мог выполнять моделирование с более широким набором данных.

У меня есть две проблемы:

1. Я не знаю, как создать интервал предсказания из подогнанной объекта GAM.
2. Генерация интервала прогнозирования с использованием posterior simulation или с использованием this technique из R-sig-mixed является вычислительно запрещающей.

Ниже приведен пример моей проблемы, генерирующий интервал, используя технику из R-sig-mixed. Будьте предупреждены, что на последнем шаге появится ошибка в том, что вы не сможете выделить огромный вектор, если только вы не находитесь на действительно впечатляющей машине.

#Some fake location data 
set.seed(13) 
x = runif(50)*2 
y = runif(50)*2 
d = cbind(ID = 1:50,as.matrix(dist(data.frame(x,y)))) 

Я хочу моделировать ссылки как функцию расстояния. Более поддельные данные:

library(reshape) 
mdata <- melt(as.data.frame(d), id=c("ID"),measure.vars = colnames(d)[2:ncol(d)],variable.name="distance") 
mdata$popularity = rnorm(25,sd=.3) 
colnames(mdata)[colnames(mdata)=="variable"] = "knows" 
colnames(mdata)[colnames(mdata)=="value"] = "distance" 
mdata = subset(mdata,ID!=knows) 
a = exp(1/(mdata$distance/runif(nrow(mdata))^mdata$distance)+mdata$popularity+rnorm(nrow(mdata),sd=.001)) 
mdata$prlink = a/(1+a) 
with(mdata,plot(distance,prlink)) 
mdata$link = runif(nrow(mdata))<mdata$prlink 
mdata$ID = as.factor(mdata$ID) 
mdata$knows = as.factor(mdata$knows) 
mdata$dum=1 #this facilitates predicting from the population of the model, later 

Теперь я модель данных:

library(mgcv) 
mod = gam(link~s(distance)+s(ID,bs="re",by=dum)+s(knows,bs="re",by=dum),data=mdata,family=binomial(link="logit")) 
plot(mod,pages=1) 
summary(mod) 

Теперь я хочу, чтобы применить подогнанную модель для моего основного набора данных:

x = runif(1000)*2 
y = runif(1000)*2 
d = cbind(ID = 1:1000,as.matrix(dist(data.frame(x,y)))) 
mdata <- melt(as.data.frame(d),id.vars = "ID") 
colnames(mdata)[colnames(mdata)=="variable"] = "knows" 
colnames(mdata)[colnames(mdata)=="value"] = "distance" 
mdata = subset(mdata,ID!=knows) 
mdata$dum=0; mdata$ID=1; mdata$knows=2 #These are needed for prediction, even though I am predicting from the population of the model, not one of the levels. 

Некоторые инструменты хронометража. ..

tic <- function(gcFirst = TRUE, type=c("elapsed", "user.self", "sys.self")) 
{ 
    type <- match.arg(type) 
    assign(".type", type, envir=baseenv()) 
    if(gcFirst) gc(FALSE) 
    tic <- proc.time()[type]   
    assign(".tic", tic, envir=baseenv()) 
    invisible(tic) 
} 

toc <- function() 
{ 
    type <- get(".type", envir=baseenv()) 
    toc <- proc.time()[type] 
    tic <- get(".tic", envir=baseenv()) 
    print(toc - tic) 
    invisible(toc) 
} 
tic() 
p = predict(mod,newdata=mdata,type="response") 
toc() 

Просто предскажите оценка очков занимает 31 секунду на моей машине. Теперь, чтобы попытаться получить интервалы предсказания, сначала нужно получить проектную матрицу ...

tic() 
Designmat = predict(mod,newdata=mdata,type="lpmatrix") 
toc() 

Это заняло у меня 47 секунд и замер мой компьютер, пока он работал.

Теперь вот способ для получения интервала предсказания, что я нашел на R-сиг-смешанный ...

ВНИМАНИЕ: Следующий код ПЫТАЙТЕСЬ ALLOCATE массивного количества памяти и может привести к сбою УСТРОЙСТВА.

tic() 
predvar <- diag(Designmat %*% vcov(mod) %*% t(Designmat)) 
SE <- sqrt(predvar) 
SE2 <- sqrt(predvar+mod$sig2) 
tfrac <- qt(0.975, mod$df.residual) 
interval = tfrac*SE2 
toc() 

>Error: cannot allocate vector of size 7435.7 Gb 

Есть ли другой способ ???

Answer 1

Вам не нужно рассчитывать Designmat %*% vcov(mod) %*% t(Designmat). В конце концов, вам нужна только диагональ. Попробуйте это:

tmp <- Designmat %*% vcov(mod) 

library(compiler) 
diagMult <- cmpfun(function(m1, m2) sapply(seq_len(nrow(m1)), 
              function(i) m1[i,] %*% m2[,i])) 
predvar <- diagMult(tmp, t(Designmat)) 

(не полностью протестирована функция должна быть реализована с Rcpp для повышения скорости, если не уже скомпилированные версии в какой-то пакет.).