← Back to branch summary

~ubuntu-branches/ubuntu/trusty/r-cran-vgam/trusty

« back to all changes in this revision

Viewing changes to R/family.ts.R

  • Committer: Package Import Robot
  • Author(s): Chris Lawrence
  • Date: 2011-11-04 13:13:06 UTC
  • mfrom: (1.2.9)
  • mto: This revision was merged to the branch mainline in revision 14.
  • Revision ID: package-import@ubuntu.com-20111104131306-w9fd83i51rw60gxf
Tags: upstream-0.8-4
Import upstream version 0.8-4

Show diffs side-by-side

added added

removed removed

Lines of Context:
R/family.ts.R
14
14
        rrar.Ak1 <- function(MM, coeffs, Ranks., aa) {
15
15
            ptr <- 0
16
16
            Ak1 <- diag(MM)
17
 
            for(j in 1:MM) {
 
17
            for(jay in 1:MM) {
18
18
                for(i in 1:MM) {
19
 
                    if (i > j && (MM+1)-(Ranks.[j]-1) <= i) {
 
19
                    if (i > jay && (MM+1)-(Ranks.[jay]-1) <= i) {
20
20
                        ptr <- ptr + 1
21
 
                        Ak1[i,j] <- coeffs[ptr]
 
21
                        Ak1[i,jay] <- coeffs[ptr]
22
22
                    }
23
23
                }
24
24
            }
26
26
            Ak1
27
27
        }
28
28
        rrar.Di <- function(i, Ranks.) {
29
 
            if (Ranks.[1]==Ranks.[i]) diag(Ranks.[i]) else 
 
29
            if (Ranks.[1] == Ranks.[i]) diag(Ranks.[i]) else 
30
30
            rbind(diag(Ranks.[i]), matrix(0, Ranks.[1]-Ranks.[i], Ranks.[i]))
31
31
        }
32
32
        rrar.Mi <- function(i, MM, Ranks., ki) {
33
 
            if (Ranks.[ki[i]]==MM)
 
33
            if (Ranks.[ki[i]] == MM)
34
34
                return(NULL)
35
35
            hi <- Ranks.[ki[i]] - Ranks.[ki[i+1]]
36
36
            Ji <- matrix(0, hi, Ranks.[1])
120
120
rrar <- function(Ranks=1, coefstart=NULL)
121
121
{
122
122
    lag.p <- length(Ranks)
 
123
 
123
124
    new("vglmff",
124
 
    blurb=c("Nested reduced-rank vector autoregressive model AR(", lag.p,
 
125
    blurb = c("Nested reduced-rank vector autoregressive model AR(", lag.p,
125
126
           ")\n\n",
126
127
           "Link:     ",
127
128
           namesof("mu_t", "identity"),
128
129
           ", t = ", paste(paste(1:lag.p, coll=",", sep="")) ,
129
130
           ""),
130
 
    initialize=eval(substitute(expression({
 
131
    initialize = eval(substitute(expression({
131
132
        Ranks. <- .Ranks
132
133
        plag <- length(Ranks.)
133
134
        nn <- nrow(x)   # original n
146
147
            ki <- udsrank <- unique(dsrank)
147
148
            uu <- length(udsrank)
148
149
            for(i in 1:uu)
149
 
               ki[i] <- max((1:plag)[dsrank==udsrank[i]])
 
150
               ki[i] <- max((1:plag)[dsrank == udsrank[i]])
150
151
            ki <- c(ki, plag+1)  # For computing a
151
152
            Ranks. <- c(Ranks., 0) # For computing a
152
153
            aa <- sum( (MM-Ranks.[ki[1:uu]]) * (Ranks.[ki[1:uu]]-Ranks.[ki[-1]]) )
162
163
 
163
164
        new.coeffs <- .coefstart  # Needed for iter=1 of $weight
164
165
        new.coeffs <- if (length(new.coeffs))
165
 
                          rep(new.coeffs, len=aa+sum(Ranks.)*MM) else
 
166
                          rep(new.coeffs, len = aa+sum(Ranks.)*MM) else
166
167
                          runif(aa+sum(Ranks.)*MM)
167
168
        temp8 <- rrar.Wmat(y.save,Ranks.,MM,ki,plag,aa,uu,nn,new.coeffs)
168
169
        X_vlm_save <- temp8$UU %*% temp8$Ht 
185
186
        w <- w[-indices]
186
187
        n.save <- n <- nn - plag
187
188
    }), list( .Ranks=Ranks, .coefstart=coefstart ))), 
188
 
    inverse=function(eta, extra=NULL) {
 
189
    linkinv=function(eta, extra=NULL) {
189
190
        aa <- extra$aa
190
191
        coeffs <- extra$coeffs
191
192
        MM <- extra$MM
262
263
 
263
264
 
264
265
garma <- function(link="identity",
265
 
                  earg=list(),
 
266
                  earg =list(),
266
267
                  p.ar.lag=1, q.lag.ma=0,
267
268
                  coefstart=NULL,
268
269
                  step=1.0)
278
279
    if (!is.list(earg)) earg = list()
279
280
 
280
281
    new("vglmff",
281
 
    blurb=c("GARMA(", p.ar.lag, ",", q.lag.ma, ")\n\n",
 
282
    blurb = c("GARMA(", p.ar.lag, ",", q.lag.ma, ")\n\n",
282
283
           "Link:     ",
283
 
           namesof("mu_t", link, earg= earg),
 
284
           namesof("mu_t", link, earg = earg),
284
285
           ", t = ", paste(paste(1:p.ar.lag, coll=",", sep=""))),
285
 
    initialize=eval(substitute(expression({
 
286
    initialize = eval(substitute(expression({
286
287
        plag <- .p.ar.lag
287
 
        predictors.names = namesof("mu", .link, earg= .earg, tag=FALSE)
 
288
        predictors.names = namesof("mu", .link, earg = .earg, tag=FALSE)
288
289
        indices <- 1:plag
289
290
        tt <- (1+plag):nrow(x) 
290
291
        pp <- ncol(x)
304
305
        w.save <- w  # Save the original
305
306
 
306
307
        new.coeffs <- .coefstart  # Needed for iter=1 of @weight
307
 
        new.coeffs <- if (length(new.coeffs)) rep(new.coeffs, len=pp+plag) else 
308
 
                      c(runif(pp), rep(0, plag)) 
 
308
        new.coeffs <- if (length(new.coeffs))
 
309
                        rep(new.coeffs, len = pp+plag) else
 
310
                        c(runif(pp), rep(0, plag)) 
309
311
        if (!length(etastart)) {
310
312
            etastart <- x[-indices,,drop=FALSE] %*% new.coeffs[1:pp]
311
313
        }
324
326
        for(i in 1:plag)
325
327
            more[[i]] <- i + max(unlist(attr(x.save, "assign")))
326
328
        attr(x, "assign") <- c(attr(x.save, "assign"), more)
327
 
    }), list( .link=link, .p.ar.lag=p.ar.lag, .coefstart=coefstart, .earg=earg ))), 
328
 
    inverse=eval(substitute(function(eta, extra=NULL) {
329
 
        eta2theta(eta, link= .link, earg= .earg)
330
 
    }, list( .link=link, .earg=earg ))),
331
 
    last=eval(substitute(expression({
 
329
    }), list( .link=link, .p.ar.lag=p.ar.lag,
 
330
              .coefstart=coefstart, .earg = earg ))), 
 
331
    linkinv = eval(substitute(function(eta, extra=NULL) {
 
332
        eta2theta(eta, link= .link, earg = .earg)
 
333
    }, list( .link=link, .earg = earg ))),
 
334
    last = eval(substitute(expression({
332
335
        misc$link <- c(mu = .link)
333
336
        misc$earg <- list(mu = .earg)
334
337
        misc$plag <- plag
335
 
    }), list( .link=link, .earg=earg ))),
336
 
    loglikelihood=eval(substitute(
 
338
    }), list( .link=link, .earg = earg ))),
 
339
    loglikelihood = eval(substitute(
337
340
        function(mu, y, w, residuals = FALSE, eta, extra=NULL) {
338
341
        if (residuals) switch( .link,
339
342
            identity=y-mu,
345
348
            loge=sum(w*(-mu + y*log(mu))),
346
349
            inverse=sum(w*(-mu + y*log(mu))),
347
350
            sum(w*(y*log(mu) + (1-y)*log1p(-mu))))
348
 
    }, list( .link=link, .earg=earg ))),
349
 
    middle2=eval(substitute(expression({
 
351
    }, list( .link=link, .earg = earg ))),
 
352
    middle2 = eval(substitute(expression({
350
353
        realfv <- fv
351
354
        for(i in 1:plag) {
352
355
            realfv <- realfv + old.coeffs[i+pp] *
354
357
        }
355
358
 
356
359
        true.eta <- realfv + offset  
357
 
        mu <- family@inverse(true.eta, extra)  # overwrite mu with correct one
358
 
    }), list( .link=link, .earg=earg ))),
359
 
    vfamily=c("garma", "vglmgam"),
360
 
    deriv=eval(substitute(expression({
 
360
        mu <- family@linkinv(true.eta, extra)  # overwrite mu with correct one
 
361
    }), list( .link=link, .earg = earg ))),
 
362
    vfamily = c("garma", "vglmgam"),
 
363
    deriv = eval(substitute(expression({
361
364
        dl.dmu <- switch( .link,
362
365
                      identity=y-mu,
363
366
                      loge=(y-mu)/mu,
364
367
                      inverse=(y-mu)/mu,
365
368
                      (y-mu) / (mu*(1-mu)))
366
 
        dmu.deta <- dtheta.deta(mu, .link, earg= .earg)
 
369
        dmu.deta <- dtheta.deta(mu, .link, earg = .earg)
367
370
        step <- .step      # This is another method of adjusting step lengths
368
371
        step * w * dl.dmu * dmu.deta
369
 
    }), list( .link=link, .step=step, .earg=earg ))),
370
 
    weight=eval(substitute(expression({
 
372
    }), list( .link=link, .step=step, .earg = earg ))),
 
373
    weight = eval(substitute(expression({
371
374
        x[,1:pp] <- x.save[tt,1:pp] # Reinstate 
372
375
 
373
376
        for(i in 1:plag) {
374
 
            temp = theta2eta(y.save[tt-i], .link, earg= .earg)
 
377
            temp = theta2eta(y.save[tt-i], .link, earg = .earg)
375
378
            x[,1:pp] <- x[,1:pp] - x.save[tt-i,1:pp] * new.coeffs[i+pp] 
376
379
            x[,pp+i] <- temp - x.save[tt-i,1:pp,drop=FALSE] %*% new.coeffs[1:pp]
377
380
        }
378
381
        class(x)=if(is.R()) "matrix" else "model.matrix" # Added 27/2/02; 26/2/04
379
382
 
380
 
        if (iter==1)
 
383
        if (iter == 1)
381
384
            old.coeffs <- new.coeffs 
382
385
 
383
386
        X_vlm_save <- lm2vlm.model.matrix(x, Blist, xij=control$xij)
387
390
                       loge=mu,
388
391
                       inverse=mu^2,
389
392
                       mu*(1-mu))
390
 
        w * dtheta.deta(mu, link= .link, earg= .earg)^2 / vary
 
393
        w * dtheta.deta(mu, link= .link, earg = .earg)^2 / vary
391
394
    }), list( .link=link,
392
 
              .earg=earg ))))
 
395
              .earg = earg ))))
393
396
}
394
397
 
395
398