关于蒙特卡洛模拟中，有变化的是方程组系数，怎么解决 模拟一次的程序见内容

新疆弹棉花的

for i = 1:100%生成1到100的X的数;
​X(i)=i;
end
X=X'; %转置;
Y = [];
X0 = ones (100,1);
Y0 = ones (100,1);
B1 = ones(100,1); %生成3个矩阵，分别为xy均值矩阵，b1的矩阵做准备;
B2 = ones(100,1);
%for i=1:1000，蒙特卡洛模拟要做1000次，所以会有1000个方程组要解;
u = normrnd(0,1,100,1);
Y = 3*B1 + 4*X ;
Ys=Y + u;
下面是指模拟一次的程序，是可以解出来的
syms s b1 b2;
% p =log(((2*pi*s)^(-50))*exp(-(1/(2*s))*sum((Ys-b1*B1-b2*X).*(Ys-b1*B1-b2*X))));
%m=diff(p,b1)*(2*s)
%n=diff(p,b2)*(2*s)
m
>>m=.......
n
>>n=.....
[b1,b2]=solve('.......=0','。。。。=0')
%联立方程组m,n，求解m=0,n=0时的b1,b2;
%只有一次的话，用solve可以，但是要把公式输进去才能算，符号不行，但是不能输1000次哇;
%求问怎么办;

hu8105

我也想知道蒙特卡洛如何解！