结果如图所示
式中x11-x1212是前文不同概率分布生成的随机数
因为要生成的样本数蛮多的，所以想用并行计算尝试一下会不会快一点
代码如下：
xx=Table[{x11[[i]],x22[[j]],x33,x44,x55,x66,x77,x88,x99,x1010,x1111,x1212},{i,1,n+1},{j,1,n+1}];//AbsoluteTiming
x=ParallelTable[{x11[[i]],x22[[j]],x33,x44,x55,x66,x77,x88,x99,x1010,x1111,x1212},{i,1,n+1},{j,1,n+1}];//AbsoluteTiming
前文代码举例：
n=100;
s=500;
u1=RandomVariate[NormalDistribution[],{s}];
x11=Table[m1[[i]]+s1*(a31[[3]]*(u1[[m]])^2+a31[[2]]*u1[[m]]+a31[[1]]),{i,1,n+1},{m,1,s}];
u3=RandomVariate[NormalDistribution[],{s}];
x33=Table[m3+s3*(a31[[3]]*(u3[[m]])^2+a31[[2]]*u3[[m]]+a31[[1]]),{m,1,s}];

并行不是提速的万金油，如果被并行的运算比较简单，核间通信的成本可能会远高于并行的收益。
你的例子只要善用列表整体操作应该就可以大幅提速，不过你的示例代码不完整，我就不折腾了。