Search code examples
rfunctiondataframemontecarlopropagation

Error propagation dataframe using Montecarlo simulation


I am doing error analysis of predictive models and I need to calculate global error, this is, I need to calculate the resultant error from propagation of indirect measurements errors. My data "df" looks like something similar to this

x         y    x_se    y_se
1.43    0.64    2.00    0.09
2.25    0.47    2.00    0.07
2.36    0.72    2.00    0.07
2.05    0.63    2.00    0.07
1.54    0.90    2.00    0.09
1.31    0.94    2.00    0.10
1.48    0.68    2.00    0.09

Where 'x' and 'y' are the measured variables, and 'x_se' and 'y_se' the standard errors of these measurements

I have used function 'propagate' from package 'qpcR' for the first row

EXPR <-  expression((0.1*x)*(y/(y+0.1)))
x <- c(1.43, 2)
y <- c(0.64, 0.09)
DF <- cbind(x, y)
res <- propagate(expr = EXPR, data = DF, type = "stat",do.sim = TRUE, verbose = TRUE)
res$summary

or this one from the package 'metRology'

expr <- expression((0.1*a)*(b/(b+0.1)))
x <- list(a=1.43, b=0.64)
u <- list(2,0.1)
u.expr<-uncert(expr, x, u, method="MC",B=200)
u.expr

both of them work, but I do not know how to apply any of these functions to the entire df. Thanks for your help


Solution

  • This is one possible solution to obtain the res$summary for each row of your dataframe.

    You first create a custom function my_fun that does what you were trying to do for a single row of the dataframe. Then, you apply this function to each row of your dataframe. The end result would be a list with as many elements as your dataframe rows.

    my_fun <- function(r){
      EXPR <- expression((0.1 * x) * (y / (y + 0.1)))
      x <- c(r[["x"]], r[["x_se"]])
      y <- c(r[["y"]], r[["y_se"]])
      DF <- cbind(x, y)
      res <- propagate(expr = EXPR, data = DF, type = "stat", do.sim = TRUE, verbose = TRUE)
      return(res$summary)
    }
    
    set.seed(1)
    apply(df, 1, my_fun)
    

    Output

    [[1]]
                      Sim Perm       Prop
    Mean        0.1233798  NaN  0.1233898
    Std.dev.    0.1726114   NA  0.1730206
    Median      0.1234929   NA         NA
    MAD         0.1708902   NA         NA
    Conf.lower -0.2153721   NA -0.2157244
    Conf.upper  0.4686238   NA  0.4625041
    
    [[2]]
                      Sim Perm       Prop
    Mean        0.1849104  NaN  0.1849310
    Std.dev.    0.1645677   NA  0.1650419
    Median      0.1851380   NA         NA
    MAD         0.1652076   NA         NA
    Conf.lower -0.1401244   NA -0.1385452
    Conf.upper  0.5155102   NA  0.5084072
    
    [[3]]
                      Sim Perm       Prop
    Mean        0.2070037  NaN  0.2070098
    Std.dev.    0.1754712   NA  0.1756395
    Median      0.2057387   NA         NA
    MAD         0.1768489   NA         NA
    Conf.lower -0.1308923   NA -0.1372374
    Conf.upper  0.5528197   NA  0.5512569
    
    [[4]]
                      Sim Perm       Prop
    Mean        0.1766511  NaN  0.1766596
    Std.dev.    0.1723788   NA  0.1726441
    Median      0.1744966   NA         NA
    MAD         0.1736575   NA         NA
    Conf.lower -0.1552201   NA -0.1617167
    Conf.upper  0.5203006   NA  0.5150359
    
    [[5]]
                      Sim Perm       Prop
    Mean        0.1384704  NaN  0.1384753
    Std.dev.    0.1798246   NA  0.1800144
    Median      0.1376255   NA         NA
    MAD         0.1782177   NA         NA
    Conf.lower -0.2151249   NA -0.2143465
    Conf.upper  0.4946514   NA  0.4912970
    
    [[6]]
                      Sim Perm       Prop
    Mean        0.1182867  NaN  0.1182874
    Std.dev.    0.1805740   NA  0.1807828
    Median      0.1184392   NA         NA
    MAD         0.1805303   NA         NA
    Conf.lower -0.2378475   NA -0.2360404
    Conf.upper  0.4731363   NA  0.4726152
    
    [[7]]
                      Sim Perm       Prop
    Mean        0.1287621  NaN  0.1287730
    Std.dev.    0.1740355   NA  0.1743982
    Median      0.1256472   NA         NA
    MAD         0.1740543   NA         NA
    Conf.lower -0.2106699   NA -0.2130411
    Conf.upper  0.4720430   NA  0.4705872