Solving a set of Linear Equations using R programm

I need to solve a system of linear equations:

aQ + bP = c
dQ + eP = f

Where:

a ~ N(100;10)
b ~ N(-1;0.1)
c ~ N(10;1)
d ~ N(10;0.1)
e ~ N(100;10)
f ~ N(10;0.1)

So far I have written:

a <- rnorm(100, mean=100, sd=10)
b <- rnorm(100, mean=-1, sd=0.1)
c <- rnorm(100, mean=10, sd=1)
d <- rnorm(100, mean=10, sd=0.1)
e <- rnorm(100, mean=100, sd=10)
f <- rnorm(100, mean=10, sd=0.1)

P <- vector()
Q <- vector()

for (i in 1:100) {
    coefs <- matrix(c(a[i],d[i],b[i],e[i]),2,2)
    ys <- array(c(c,f),2)
    solve(coefs[i], ys[i])
}

The problem is that the for loop is only giving me one solution for P and Q and I would like that the for loop to calculate the 100 set of values do a,b,c,d,e and f and store the data on the vectors Q and P.

回答1:

You could try using apply()

# data  
df = data.frame(a,b,c,d,e,f)

# apply function
out = t(apply(df, 1, 
      function(x){ 
      coefs = matrix( c(x['a'], x['d'], x['b'], x['e']), 2, 2); 
      ys = array(c(x['c'],x['f']),2); 
      out = solve(coefs, ys); 
      names(out) = c('P','Q');
      out}))

Or else using sapply()

out = t(sapply(seq(100), function(i){
       coefs = matrix(c(a[i], d[i], b[i] ,e[i]),2,2); ys = array(c(c[i],f[i]),2);
       out = solve(coefs, ys); names(out) = c('P','Q'); 
       out}))

and if you want to use for loop, here is what you could do

# declare matrix to store the output
out = matrix(ncol=2, nrow=100) 

# populate declared matrix using for loop
for(i in 1:100){  
    coefs = matrix(c(a[i],d[i],b[i],e[i]),2,2)
    ys = array(c(c[i],f[i]),2)
    out[i,] = as.vector(solve(coefs, ys))
    colnames(out) = c('P','Q')
    out}

回答2:

Just to offer another possible approach, you could use mapply() to effectively iterate over each coefficient and RHS vector in parallel:

set.seed(1);
a <- rnorm(100,100,10);
b <- rnorm(100,-1,0.1);
c <- rnorm(100,10,1);
d <- rnorm(100,10,0.1);
e <- rnorm(100,100,10);
f <- rnorm(100,10,0.1);
mapply(function(a,b,c,d,e,f) solve(matrix(c(a,d,b,e),2),matrix(c(c,f),2)),a,b,c,d,e,f);
##            [,1]       [,2]       [,3]      [,4]       [,5]       [,6]       [,7]       [,8]       [,9]      [,10]      [,11]      [,12]      [,13]      [,14]      [,15]      [,16]      [,17]      [,18]      [,19]      [,20]      [,21]      [,22]      [,23]     [,24]      [,25]      [,26]      [,27]      [,28]      [,29]      [,30]     [,31]      [,32]      [,33]      [,34]      [,35]      [,36]      [,37]      [,38]      [,39]      [,40]      [,41]      [,42]      [,43]     [,44]      [,45]      [,46]      [,47]      [,48]      [,49]     [,50]      [,51]      [,52]      [,53]     [,54]      [,55]      [,56]      [,57]     [,58]      [,59]      [,60]      [,61]      [,62]      [,63]      [,64]      [,65]     [,66]      [,67]      [,68]      [,69]      [,70]      [,71]     [,72]      [,73]      [,74]      [,75]     [,76]      [,77]      [,78]      [,79]      [,80]      [,81]      [,82]      [,83]      [,84]      [,85]      [,86]      [,87]      [,88]     [,89]     [,90]     [,91]      [,92]      [,93]      [,94]      [,95]      [,96]      [,97]      [,98]      [,99]     [,100]
## [1,] 0.11195915 0.11550647 0.12751804 0.0841975 0.07569079 0.13696420 0.10252485 0.09894877 0.09515799 0.10922284 0.08629273 0.10113763 0.10308528 0.11216529 0.09963416 0.11659483 0.09791256 0.08081292 0.09926686 0.09472402 0.07648024 0.09355051 0.09391840 0.1219382 0.08409615 0.11939194 0.09931966 0.09958410 0.10817085 0.09947289 0.0803567 0.07280706 0.09108175 0.10730635 0.11623282 0.10439468 0.11082056 0.08964156 0.10089985 0.09365293 0.10998488 0.11389442 0.09656601 0.0872372 0.12093390 0.08759978 0.09202209 0.09137391 0.10043375 0.1025373 0.09823991 0.11186173 0.09698966 0.1112308 0.09463053 0.09386283 0.07966344 0.1193879 0.09922995 0.09776010 0.08903184 0.09724774 0.09156914 0.10903049 0.12760021 0.1014739 0.11805636 0.07770935 0.09601201 0.07510938 0.09370570 0.1130247 0.08705415 0.14043486 0.11704624 0.1091617 0.08154497 0.10810454 0.08690775 0.11685334 0.11133037 0.09817079 0.10215157 0.11097263 0.08957527 0.08781713 0.08506901 0.11398943 0.1017468 0.1080845 0.1026900 0.09321191 0.09276108 0.08089681 0.10226763 0.09705664 0.12423443 0.11736748 0.11432446 0.10282015
## [2,] 0.08016778 0.07420647 0.09132015 0.0953508 0.09734921 0.09111316 0.09172938 0.09253545 0.08018053 0.09452145 0.09233637 0.08349215 0.07979583 0.09634328 0.09443451 0.08404651 0.08153545 0.09468065 0.10482579 0.07829014 0.08002979 0.09773864 0.09023393 0.1070831 0.08657002 0.07433625 0.10757136 0.09891601 0.09612082 0.09671603 0.1139721 0.08951823 0.10746614 0.08889339 0.08501627 0.09050828 0.08181816 0.09255519 0.09681833 0.08489033 0.09251537 0.07567294 0.09099927 0.0841442 0.08769969 0.13057528 0.10425570 0.09546615 0.08564604 0.1174630 0.08258592 0.09420732 0.09640516 0.1048073 0.10757819 0.08931113 0.08544322 0.1124914 0.09978848 0.08516907 0.09438840 0.07408315 0.08145526 0.08357922 0.08661255 0.0884171 0.09454044 0.08658720 0.07824592 0.08435285 0.09760972 0.1068059 0.09319386 0.08188822 0.07865432 0.1037524 0.08916170 0.08914854 0.07982980 0.08718701 0.09064425 0.10193109 0.09645453 0.09396509 0.07402690 0.07399806 0.09308106 0.09940533 0.1126132 0.0841854 0.1002329 0.07365820 0.10095744 0.08985023 0.06348183 0.10252098 0.08397957 0.09875275 0.08654524 0.09736881

Or if you prefer getting the results as a list, you could pass SIMPLIFY=F to mapply(), or just use Map(), which doesn't simplify by default:

Map(function(a,b,c,d,e,f) solve(matrix(c(a,d,b,e),2),matrix(c(c,f),2)),a,b,c,d,e,f);
## [[1]]
##            [,1]
## [1,] 0.11195915
## [2,] 0.08016778
##
## [[2]]
##            [,1]
## [1,] 0.11550647
## [2,] 0.07420647
##
## ... (snip) ...
##
## [[99]]
##            [,1]
## [1,] 0.11432446
## [2,] 0.08654524
##
## [[100]]
##            [,1]
## [1,] 0.10282015
## [2,] 0.09736881