[Analogue-commits] r286 - in pkg: R inst tests/Examples

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Sat Sep 8 00:48:00 CEST 2012


Author: gsimpson
Date: 2012-09-08 00:47:59 +0200 (Sat, 08 Sep 2012)
New Revision: 286

Modified:
   pkg/R/stdError.R
   pkg/inst/ChangeLog
   pkg/tests/Examples/analogue-Ex.Rout.save
Log:
weighted variance calc assumed weights summed to 1

Modified: pkg/R/stdError.R
===================================================================
--- pkg/R/stdError.R	2012-08-15 14:14:50 UTC (rev 285)
+++ pkg/R/stdError.R	2012-09-07 22:47:59 UTC (rev 286)
@@ -31,17 +31,21 @@
     ords <- apply(object$Dij, 2, getOrd)
     SEQ <- seq_len(ncol(ords))
     ## weights = 1/Dij
-    wts <- 1 / sapply(SEQ, getWts, object$Dij, ords, k.seq, USE.NAMES = FALSE)
+    wi <- 1 / sapply(SEQ, getWts, object$Dij, ords, k.seq, USE.NAMES = FALSE)
     ## produce matrix of Env data for each site
     env <- sapply(SEQ, getEnv, object$Dij, ords, k.seq,
                   object$orig.y, USE.NAMES = FALSE)
     ## mean of env of k closest analogues
     ybar <- colMeans(env)
-    wtdSD <- sqrt(colSums(wts * sweep(env, 2, ybar, "-")^2) /
-                  colSums(wts))
+    ## sum weights
+    sum.wi <- colSums(wi)
+    sum.wi2 <- colSums(wi^2)
+    sum2.wi <- sum.wi^2
+    frac <- sum.wi / (sum2.wi - sum.wi2)
+    wtdSD <- sqrt(frac * colSums(wi * sweep(env, 2, ybar, "-")^2))
     names(wtdSD) <- names(object$orig.y)
     class(wtdSD) <- "stdError"
-    return(wtdSD)
+    wtdSD
 }
 
 `stdError.predict.mat` <- function(object, k, ...) {
@@ -69,17 +73,21 @@
     ords <- apply(object$Dij, 2, getOrd)
     SEQ <- seq_len(ncol(ords))
     ## weights = 1/Dij
-    wts <- 1 / sapply(SEQ, getWts, object$Dij, ords, k.seq, USE.NAMES = FALSE)
+    wi <- 1 / sapply(SEQ, getWts, object$Dij, ords, k.seq, USE.NAMES = FALSE)
     ## produce matrix of Env data for each site
     env <- sapply(SEQ, getEnv, object$Dij, ords, k.seq,
                   object$observed, USE.NAMES = FALSE)
     ## mean of env of k closest analogues
     ybar <- colMeans(env)
-    wtdSD <- sqrt(colSums(wts * sweep(env, 2, ybar, "-")^2) /
-                  colSums(wts))
+    ## sum weights
+    sum.wi <- colSums(wi)
+    sum.wi2 <- colSums(wi^2)
+    sum2.wi <- sum.wi^2
+    frac <- sum.wi / (sum2.wi - sum.wi2)
+    wtdSD <- sqrt(frac * colSums(wi * sweep(env, 2, ybar, "-")^2))
     names(wtdSD) <- colnames(object$predictions$model$predicted)
     class(wtdSD) <- "stdError"
     attr(wtdSD, "k") <- k
     attr(wtdSD, "auto") <- object$auto
-    return(wtdSD)
+    wtdSD
 }

Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog	2012-08-15 14:14:50 UTC (rev 285)
+++ pkg/inst/ChangeLog	2012-09-07 22:47:59 UTC (rev 286)
@@ -2,6 +2,8 @@
 
 Version 0.9-10
 
+	* stdError: calculation assumed weights summed to 1.
+
 	* caterpillarPlot: can now be called by the shorter name
 	caterpillar().
 

Modified: pkg/tests/Examples/analogue-Ex.Rout.save
===================================================================
--- pkg/tests/Examples/analogue-Ex.Rout.save	2012-08-15 14:14:50 UTC (rev 285)
+++ pkg/tests/Examples/analogue-Ex.Rout.save	2012-09-07 22:47:59 UTC (rev 286)
@@ -1,5 +1,5 @@
 
-R version 2.15.1 Patched (2012-07-27 r60018) -- "Roasted Marshmallows"
+R version 2.15.1 Patched (2012-08-29 r60485) -- "Roasted Marshmallows"
 Copyright (C) 2012 The R Foundation for Statistical Computing
 ISBN 3-900051-07-0
 Platform: x86_64-unknown-linux-gnu (64-bit)
@@ -6528,21 +6528,21 @@
 > ## standard errors
 > stdError(ik.mat)
    V14.61   V17.196   V18.110   V16.227    V14.47    V23.22     V2.12    V23.29 
-2.9570077 2.2342181 2.0330039 2.8513418 2.2138447 2.4593532 2.3161163 0.8285100 
+3.7017647 2.7940379 2.4904961 3.5144020 2.7300757 3.0418167 2.9199864 1.0952777 
    V12.43      R9.7    A157.3    V23.81    V23.82    V12.53    V23.83    V12.56 
-1.8521500 2.9870952 1.3093628 1.4874007 1.5322853 2.1145062 0.4004632 2.4418377 
+2.4168224 3.8347560 1.6422409 1.8922602 1.9662188 2.6825814 0.5043179 2.9932954 
   A152.84    V16.50   V22.122    V16.41     V4.32    V12.66   V19.245      V4.8 
-2.7213261 1.8144450 1.6398700 2.2182611 0.0000000 0.8432195 1.9932477 1.7840542 
+3.3706748 2.2382896 2.0969686 3.0440301 0.0000000 1.0495554 2.6358967 2.1929384 
   A180.15    V18.34   V20.213   V19.222   A180.39   V16.189    V12.18     V7.67 
-2.3981335 1.3210885 2.0899961 1.8534830 0.8091183 0.5628718 0.1068089 0.8383123 
+2.9604059 1.6382546 2.5610904 2.2776091 0.9980639 0.6977540 0.1364087 1.0693862 
   V17.165   V19.310   V16.190  A153.154   V19.308   V22.172    V10.98   V22.219 
-0.9361208 0.8748343 1.6496071 0.4790917 0.4814322 0.4293269 1.2054665 0.4599242 
+1.1469297 1.0946184 2.0219976 0.6101154 0.6061351 0.5269860 1.4875516 0.5643260 
    V16.33   V22.204   V20.167    V10.89    V12.79   V19.216    V14.90   A180.72 
-1.2390022 1.1486092 1.1037839 0.9112438 0.4250058 0.8771831 0.6509596 0.4911457 
+1.5194344 1.4172128 1.3626793 1.1202298 0.5381647 1.0771118 0.8011406 0.6104648 
    V16.21   A180.76   V15.164   A180.78     V14.5    V3.128   A179.13     V9.31 
-0.4907379 0.2565399 0.3519882 0.6308207 0.4650968 1.1928807 0.4847921 0.2452820 
+0.6381386 0.3198533 0.4701962 0.7742820 0.5714590 1.4633510 0.5972338 0.3021622 
   V20.230     V20.7   V20.234    V18.21   V12.122 
-0.2428485 0.4741689 0.4771624 0.4564864 0.0000000 
+0.2980597 0.5809639 0.5854821 0.5723716 0.0000000 
 attr(,"class")
 [1] "stdError"
 > 
@@ -6550,34 +6550,38 @@
 > coreV12.mat <- predict(ik.mat, V12.122, k = 3)
 > ## standard errors
 > stdError(coreV12.mat)
-        0        10        20        30        40        50        60        70 
-0.6666560 0.5127652 0.4940833 1.3084467 1.0824814 1.0650554 0.7312066 1.0129625 
-       80        90       100       110       120       130       140       150 
-1.1916917 1.1149534 0.4033714 0.2334125 1.3081399 0.2394602 1.4197788 0.7572307 
-      160       170       180       190       200       210       220       230 
-0.6499565 0.6538152 0.9742193 0.2401338 1.3854266 1.8150936 0.4194549 0.5012481 
-      240       250       260       270       280       290       300       310 
-1.1805085 1.0541290 1.0460351 1.2970070 1.5004813 0.4862513 1.0462475 0.9377745 
-      320       330       340       350       360       370       380       390 
-1.2997606 1.1837605 1.0368135 0.6092353 1.0319833 1.0469552 0.4029888 1.8028169 
-      400       410       420       430       440       450       460       470 
-0.3996499 1.1499003 0.2427915 0.4139073 0.2583339 0.0000000 1.2141587 1.0439356 
-      480       490       500       510       520       530       540       550 
-0.2338637 0.9428443 0.4827475 0.4148127 0.3879866 1.0328590 0.0937982 0.8497036 
-      560       570       580       590       600       610       620       630 
-0.2416557 0.2278804 0.4058629 0.9718368 0.2303837 0.7352933 0.3277687 0.2261525 
-      640       650       660       670       680       690       700       710 
-0.4791844 0.6390568 0.2347800 0.4065762 1.8441614 1.0495782 0.7924059 1.1846794 
-      720       730       740       750       760       770       780       790 
-1.2390412 0.8534425 0.4012654 0.2426828 0.2340115 1.1977802 0.9582384 0.4754369 
-      800       810       820       830       840       850       860       870 
-0.4855805 1.4611535 0.2325244 1.1660600 0.8679123 0.4122362 1.0629969 1.0521058 
-      880       890       900       910       920       930       940       950 
-1.2560759 1.6556741 0.5943474 0.8451956 0.9962813 0.2244707 0.2450447 0.3782048 
-      960       970       980       990      1000      1010      1020      1030 
-0.4063936 0.2338417 0.2274927 0.2340763 1.1603260 0.2382300 0.8493407 1.1706837 
-     1040      1050      1060      1070      1080      1090 
-0.4089408 1.0715114 0.2426835 1.0801634 1.4499069 1.0887329 
+         0         10         20         30         40         50         60 
+72.3307103  0.6391695  0.6113086  1.6025626  1.3348327  1.3084087  0.8978862 
+        70         80         90        100        110        120        130 
+ 1.2413641  1.4600559  1.3768654  0.4941674  0.2879595  1.6022391  0.2934360 
+       140        150        160        170        180        190        200 
+ 1.7397088  0.9275737  0.7977268  0.8014461  1.1946764  0.2943270  1.7023691 
+       210        220        230        240        250        260        270 
+ 2.2240272  0.5164777  0.6192205  1.4696838  1.3014299  1.2840133  1.5909909 
+       280        290        300        310        320        330        340 
+ 1.8396982  0.5974174  1.2814680  1.1491139  1.5922597  1.4498985  1.2714704 
+       350        360        370        380        390        400        410 
+ 0.7462475  1.2643008  1.2826104  0.4937233  2.2117480  0.4924661  1.4086468 
+       420        430        440        450        460        470        480 
+ 0.2979207  0.5081370  0.3230163  0.0000000  1.4923076  1.2789171  0.2864583 
+       490        500        510        520        530        540        550 
+ 1.1551482  0.5925977  0.5086959  0.4774772  1.2652373  0.1148866  1.0407605 
+       560        570        580        590        600        610        620 
+ 0.2967980  0.2796920  0.4972385  1.2027792  0.2835480  0.9006515  0.4014728 
+       630        640        650        660        670        680        690 
+ 0.2813064  0.5873118  0.7853038  0.2875564  0.4981700  2.2594850  1.2891113 
+       700        710        720        730        740        750        760 
+ 0.9713892  1.4526717  1.5179814  1.0454000  0.4918791  0.2977649  0.2869260 
+       770        780        790        800        810        820        830 
+ 1.4679702  1.1750746  0.5846293  0.5960886  1.7904245  0.2851860  1.4506051 
+       840        850        860        870        880        890        900 
+ 1.0694963  0.5053518  1.3027428  1.2888246  1.5458275  2.0294194  0.7310171 
+       910        920        930        940        950        960        970 
+ 1.0429069  1.2370137  0.2763181  0.3011330  0.4632124  0.4977645  0.2866663 
+       980        990       1000       1010       1020       1030       1040 
+ 0.2804741  0.2868464  1.4217870  0.2919910  1.0469965  1.4342437  0.5022240 
+      1050       1060       1070       1080       1090 
+ 1.3229458  0.2985295  1.3235215  1.7766052  1.3336282 
 attr(,"class")
 [1] "stdError"
 attr(,"k")
@@ -7250,7 +7254,7 @@
 > ### * <FOOTER>
 > ###
 > cat("Time elapsed: ", proc.time() - get("ptime", pos = 'CheckExEnv'),"\n")
-Time elapsed:  19.233 0.237 20.124 0 0 
+Time elapsed:  18.418 0.198 19.118 0 0 
 > grDevices::dev.off()
 null device 
           1 



More information about the Analogue-commits mailing list