← Back to branch summary

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

« back to all changes in this revision

Viewing changes to R/family.robust.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.robust.R
10
10
 
11
11
 
12
12
 
 
13
 
 
14
 
13
15
edhuber <- function(x, k = 0.862, mu = 0, sigma = 1, log = FALSE) {
14
16
  if (!is.logical(log.arg <- log)) stop("bad input for argument 'log'")
15
17
  rm(log)
130
132
 
131
133
  if (!is.Numeric(imethod, allow = 1, integ = TRUE, posit = TRUE) ||
132
134
      imethod > 4)
133
 
       stop("'imethod' must be 1 or 2 or 3 or 4")
 
135
       stop("argument 'imethod' must be 1 or 2 or 3 or 4")
134
136
 
135
137
  if (!is.Numeric(k, allow = 1, posit = TRUE))
136
138
      stop("bad input for argument 'k'")
179
181
      }
180
182
    }), list( .llocat = llocation, .lscale = lscale,
181
183
              .elocat = elocation, .escale = escale,
182
 
              .imethod=imethod ))),
183
 
    inverse = eval(substitute(function(eta, extra = NULL) {
 
184
              .imethod = imethod ))),
 
185
    linkinv = eval(substitute(function(eta, extra = NULL) {
184
186
      eta2theta(eta[,1], .llocat, earg = .elocat)
185
187
    }, list( .llocat = llocation,
186
188
             .elocat = elocation, .escale = escale ))),
257
259
 
258
260
 
259
261
 
 
262
 
 
263
 huber1 <- function(llocation = "identity",
 
264
                    elocation = list(),
 
265
                    k = 0.862,
 
266
                    imethod = 1) {
 
267
 
 
268
 
 
269
 print("hi 20110802")
 
270
  A1 <- (2 * dnorm(k) / k - 2 * pnorm(-k))
 
271
  eps <- A1 / (1 + A1)
 
272
 
 
273
  if (!is.Numeric(imethod, allow = 1, integ = TRUE, posit = TRUE) ||
 
274
      imethod > 4)
 
275
       stop("argument 'imethod' must be 1 or 2 or 3 or 4")
 
276
 
 
277
  if (!is.Numeric(k, allow = 1, posit = TRUE))
 
278
      stop("bad input for argument 'k'")
 
279
 
 
280
  if (mode(llocation)  !=  "character" && mode(llocation) != "name")
 
281
       llocation = as.character(substitute(llocation))
 
282
  if (!is.list(elocation)) elocation = list()
 
283
 
 
284
  new("vglmff",
 
285
    blurb = c("Huber least favorable distribution\n\n",
 
286
              "Links: ",
 
287
              namesof("location",  llocation,  earg = elocation), "\n\n",
 
288
              "Mean: location"),
 
289
    initialize = eval(substitute(expression({
 
290
      predictors.names <-
 
291
         c(namesof("location", .llocat, earg = .elocat, tag = FALSE))
 
292
 
 
293
      if (ncol(y <- cbind(y)) != 1)
 
294
           stop("response must be a vector or a one-column matrix")
 
295
 
 
296
      if (!length(etastart)) {
 
297
          junk = lm.wfit(x = x, y = y, w = w)
 
298
          location.init <- if ( .llocat == "loge") pmax(1/1024, y) else {
 
299
            if ( .imethod == 3) {
 
300
              rep(weighted.mean(y, w), len = n)
 
301
            } else if ( .imethod == 2) {
 
302
              rep(median(rep(y, w)), len = n)
 
303
            } else if ( .imethod == 1) {
 
304
              junk$fitted
 
305
            } else {
 
306
              y
 
307
            }
 
308
          }
 
309
          etastart <- cbind(
 
310
               theta2eta(location.init,  .llocat, earg = .elocat))
 
311
      }
 
312
    }), list( .llocat = llocation,
 
313
              .elocat = elocation,
 
314
              .imethod = imethod ))),
 
315
    linkinv = eval(substitute(function(eta, extra = NULL) {
 
316
      eta2theta(eta, .llocat, earg = .elocat)
 
317
    }, list( .llocat = llocation,
 
318
             .elocat = elocation ))),
 
319
    last = eval(substitute(expression({
 
320
      misc$link <-    c("location" = .llocat )
 
321
      misc$earg <- list("location" = .elocat )
 
322
      misc$expected <- TRUE
 
323
      misc$k.huber <- .k
 
324
      misc$imethod <- .imethod
 
325
    }), list( .llocat = llocation,
 
326
              .elocat = elocation,
 
327
              .k      = k,         .imethod = imethod ))),
 
328
   loglikelihood = eval(substitute(
 
329
     function(mu, y, w, residuals = FALSE, eta, extra = NULL) {
 
330
     location <- eta2theta(eta, .llocat, earg = .elocat)
 
331
     kay      <- .k
 
332
     if (residuals) stop("loglikelihood residuals not ",
 
333
                         "implemented yet") else {
 
334
       sum(w * dhuber(y, k = kay, mu = location,  sigma = 1,
 
335
                      log = TRUE))
 
336
     }
 
337
   }, list( .llocat = llocation,
 
338
            .elocat = elocation,
 
339
            .k      = k ))),
 
340
    vfamily = c("huber1"),
 
341
    deriv = eval(substitute(expression({
 
342
      mylocat <- eta2theta(eta, .llocat,  earg = .elocat)
 
343
      myk     <- .k
 
344
 
 
345
      zedd <- (y - mylocat) # / myscale
 
346
      cond2 <- (abs(zedd) <=  myk)
 
347
      cond3 <-     (zedd  >   myk)
 
348
 
 
349
      dl.dlocat        <- -myk + 0 * zedd # cond1
 
350
      dl.dlocat[cond2] <- zedd[cond2]
 
351
      dl.dlocat[cond3] <-  myk  # myk is a scalar
 
352
      dl.dlocat <- dl.dlocat # / myscale
 
353
 
 
354
 
 
355
    if (FALSE) {
 
356
      dl.dscale        <- (-myk * zedd)
 
357
      dl.dscale[cond2] <-      (zedd^2)[cond2]
 
358
      dl.dscale[cond3] <- ( myk * zedd)[cond3]
 
359
      dl.dscale <- (-1 + dl.dscale) / myscale
 
360
    }
 
361
 
 
362
      dlocat.deta <- dtheta.deta(mylocat, .llocat, earg = .elocat)
 
363
      ans <-
 
364
      c(w) * cbind(dl.dlocat * dlocat.deta)
 
365
      ans
 
366
    }), list( .llocat = llocation,
 
367
              .elocat = elocation,
 
368
              .eps    = eps,       .k      = k ))),
 
369
    weight = eval(substitute(expression({
 
370
      wz   <- matrix(as.numeric(NA), n, 1) # diag matrix; y is one-col too
 
371
 
 
372
 
 
373
 
 
374
 
 
375
      temp4 <- erf(myk / sqrt(2))
 
376
      ed2l.dlocat2 <- temp4 * (1 - .eps) # / myscale^2
 
377
 
 
378
 
 
379
      wz[, iam(1,1,M)] <- ed2l.dlocat2 * dlocat.deta^2
 
380
      ans
 
381
      c(w) * wz
 
382
    }), list( .eps = eps ))))
 
383
}
 
384
 
 
385
 
 
386