模拟满足R中线性方程的数据吗？

Question

我想生成满足以下方程的数据（w1，w2，w3，w4）：

w1 + w2 + w3 + w4 = 1使得对于所有i，wi> = 0。

那么我会如何继续？

Answer 1

gend=function(n){#n is the number of such samples you would like to generate. 
    D=matrix(NA,n,4) 
    for(i in 1:n){ 
    a=runif(3) 
    b=a[2];d=a[3];a=a[1] 
    while((b+a)>1)b=runif(1) 
    while((b+a+d)>1)d=runif(1) 
    D[i,]=c(a,b,d,1-sum(a,b,d)) 
    } 
D 
} 
S=gend(10) 
S 
      [,1]  [,2]  [,3]   [,4] 
[1,] 0.8599033 0.04730250 0.008108873 0.0846853156 
[2,] 0.3181092 0.37239277 0.243917543 0.0655804751 
[3,] 0.4733250 0.07116504 0.145719661 0.3097903084 
[4,] 0.9482858 0.04967905 0.001082051 0.0009531255 
[5,] 0.7291447 0.16336750 0.106824697 0.0006631124 
[6,] 0.1896754 0.60669666 0.062790600 0.1408372964 
[7,] 0.8452493 0.12875891 0.022590556 0.0034
[8,] 0.4452403 0.52793687 0.011378945 0.0154438773 
[9,] 0.4618185 0.23065407 0.224276892 0.0832504972 
[10,] 0.2191794 0.09512862 0.110790973 0.5749009629 

rowSums(S) 
[1] 1 1 1 1 1 1 1 1 1 1 

Answer 2

没有循环，你可以尝试这样的事：

require(data.table) 
n <- 200 # how much data you need 
myFunc <- function(n) { 
    v <- replicate(3, rnorm(1e6)) # from normal distribution simulate 3 variables 
    v <- as.data.table(v) 
    v <- v[V1 >= 0 & V2 >= 0 & V3 >= 0] 
    v[, s := V1 + V2 + V3] 
    v <- v[s <= 1] 
    v[, V4 := 1 - s] 
    v[, s := NULL] 
    v <- v[1:n] # select needed count 
    v[] 
} 
set.seed(32) 
myFunc(n 
#    V1   V2   V3   V4 
# 1: 0.167742364 0.10237134 0.54388085 0.186005445 
# 2: 0.268182684 0.15659897 0.23427071 0.340947639 
# 3: 0.449311437 0.35508438 0.19126053 0.004343644 
# 4: 0.248128979 0.18794309 0.20178770 0.362140224 
# 5: 0.235724248 0.02147451 0.68609822 0.056703023 
# ---            
# 196: 0.748858750 0.10044298 0.14843887 0.002259405 
# 197: 0.527629846 0.15944044 0.11169268 0.2
# 198: 0.406024037 0.02289703 0.33207412 0.239004809 
# 199: 0.024428994 0.31282766 0.61436293 0.048380418 
# 200: 0.008766677 0.54692203 0.07853701 0.365774287 

all(sapply(myFunc(n), function(x) all(x <= 1 & x >=0))) 
# [1] TRUE 
all.equal(rowSums(myFunc(n)), rep(1, n)) 
# [1] TRUE 

速度对比：

microbenchmark::microbenchmark(gend(n), 
           myF(n)) 
# Unit: milliseconds 
# expr  min  lq  mean median  uq  max neval cld 
# gend(n) 3.996232 7.034548 21.23260 10.32951 16.37486 272.49600 100 b 
# myF(n) 7.723631 9.177185 11.00396 11.27347 13.05472 14.43866 100 a