Математическое моделирование в среде MATLAB компьютерный эксперимент на основе бутстреп – подхода (метода), страница 4

    4.9254  -17.7279    8.4965  166.0618

    5.0167  -17.4568    8.4011  164.8657

    6.4715  -17.0078    7.7622  165.8663

    6.1602  -15.8387    9.3116  162.9867

    8.7723  -16.4936    9.9522  167.6259

    5.5873  -20.0521    8.8876  164.5282

    9.3188  -16.3801    8.6333  164.9673

    8.2628  -16.2089    9.3292  168.0092

    3.0824  -16.4716    9.1210  164.6361

    4.8471  -20.0319   10.9658  164.4028

    6.3940  -18.3714    9.4913  167.7776

    5.9826  -17.7898    7.9912  166.5001

    4.4461  -17.4977    8.7767  167.7601

    7.9668  -16.2557    5.3272  165.0921

    6.1681  -18.5945    9.7106  164.1607

    6.9955  -18.5933    9.1761  167.0502

    5.1166  -17.8129    8.1140  165.9005

    5.2079  -17.5418    8.0186  164.7044

    6.6627  -17.0928    7.3797  165.7050

    6.3514  -15.9236    8.9291  162.8255

    8.9635  -16.5786    9.5697  167.4646

    5.7785  -20.1370    8.5052  164.3669

    9.5100  -16.4651    8.2508  164.8060

    8.4540  -16.2939    8.9467  167.8480

    3.2736  -16.5566    8.7386  164.4749

3. Рассчитать оценки параметров модели методом МНК (взвешенным) в соответствии со снятыми данными, провести сравнительный анализ оценок параметров модели.

Для расчета оценок воспользоваться алгоритмами предложенными в лабораторной работе №3, п.4 Для сравнительного анализа оценок параметров модели воспользоваться разностью (), которая показывает насколько полученные оценки параметров модели после проведения компьютерных экспериментов с бутстреп выборкой стали ближе к истинным параметрам.

Просчитать 100  раз, и определить сколько раз .

3.1. По гистограмме распределения ошибок

sigma = 1.5;

Q=[6 -4 -3 0.5];

N = 16;

rnd=[];

Ysred=[];

Otkl=[];

k_intrv=5;

KB=100;

RNX=[];

RNXT=[];

YB=[];

Yall=[];

MeanYB=[];

MeanYall=[];

p=[];

maxp=0;

sq_m=0;

PLAN=[

            0      16;

            2.765  16;

            7.235  16;

            10     16;

];

for l= 1: KB

random = randn(16,4)*sigma;

Ymodel=zeros(N,4);

for j=1:N

            for i=1:4

              Ymodel(j,i) = Q(1)+Q(2)*PLAN(i,1)+Q(3)*PLAN(i,1)^2+Q(4)*PLAN(i,1)^3+random(j,i);

            end

end

Ysred = mean(Ymodel);

%Ysred

for i=1:N

   for j=1:4

      Otkl(i,j)= Ymodel(i,j) - Ysred(j);

  end

end

%Otkl

MinOtkl = min(min(Otkl));

MaxOtkl = max(max(Otkl));

Delta = (MaxOtkl - MinOtkl)/k_intrv;

Chast = zeros (1,k_intrv);

LeftVal = MinOtkl;

for l = 1:k_intrv     

   for i = 1:N

      for j = 1:4

            if (Otkl(i,j) >= LeftVal) & (Otkl(i,j) < LeftVal + Delta)

            Chast(l) = Chast(l) + 1;

         end

         if Otkl(i,j) == MaxOtkl & l == k_intrv

            Chast(k_intrv) = Chast(k_intrv) + 1;

         end

      end  

   end   

   LeftVal = LeftVal + Delta;

end

%Chast;

%point1

for i= 1: k_intrv

   p(i)=Chast(i)/(N*4);

end

maxP=max(p);

sum=0;

N1=N;

   for i= 1: 4

      sum_b=zeros(1,k_intrv);

      k=0;

     while N~=0

       r_x = rand;

       r_y = rand;

       ran_x = MinOtkl+(MaxOtkl-MinOtkl)*r_x;

       ran_y = maxP*r_y;

       for j= 1: k_intrv

          a = MinOtkl + (j-1)*Delta;

          b = MinOtkl + j*Delta;

          if (a <= ran_x) & (ran_x < b)

            if  ran_y <= p(j)

             if  sum_b(j) < Chast(j)

                        sum_b(j) = sum_b(j) + 1;

                           k = k + 1;

               RNX(i,k) = ran_x;

               N = N-1;

               break;

             end

            else break;

           end

         end

       end

     end

   N = N1;  

   end

   RNXT = RNX';

   for i=1:N

     for j=1:4

       YB(i,j) = Ysred(j) + RNXT(i,j);

     end

   end

            Yall = [Ymodel;YB];

            MeanYB = mean(YB);

            MeanYall = mean(Yall);

      n = 4;

            k = 4;

            m = 1;

            X = zeros(1,m);

            f=[1 X(1) X(1)^2 X(1)^3];

            F = ones(n,k);

            for i= 1: n

               for j= 1: m

         X(j) = PLAN(i,j);

               end

               f=[1 X(1) X(1)^2 X(1)^3];

               for j= 1: k

         F(i,j) = f(j);

      end

   end

   Q_M=inv(F'*F)*F'*Ysred';

   D_MB=[];

   for i= 1: k

            sum = 0;

            for j= 1: 2*N

            sum = sum + (Yall(j,i)-MeanYall(i))^2;

            end

            D_MB(i) = (1/((2*N)*(2*N-1)))*sum;

   end

   D = diag(D_MB);

   M = F'*inv(D)*F;

   Q_MB=inv(M)*F'*inv(D)*MeanYall';

   sqm=0;

   sqmb=0;

   for i= 1: k

      sqm = sqm + (Q(i)-Q_M(i))^2;

      sqmb = sqmb + (Q(i)-Q_MB(i))^2;

   end

   if sqm > sqmb

     sq_m=sq_m+1;

  end

end

YB

Yall

MeanYB

MeanYall

sq_m

Q_M

Q_MB

Результаты

YB =

    6.3445  -18.5623    6.7798  165.3493

    7.1925  -18.2450    7.1858  163.9161

    6.1192  -17.4951    9.3537  164.3617

    6.8087  -18.9182   11.1488  167.2663

    5.0421  -16.0366    9.3526  165.7689

    8.5740  -17.6436    9.8665  165.6685

    8.7000  -15.7616    9.0134  166.5463

    6.3938  -16.6296    9.6259  167.1241

    9.3566  -19.5029    9.8847  165.7145

    7.1967  -17.3653    9.1331  166.3693

    6.8748  -15.5068    6.9022  165.3785

    6.3774  -16.0153    9.3000  165.8999

    6.5430  -15.8328    9.3214  166.5810

    7.2523  -17.0393    9.2903  163.8583

    8.1870  -14.3908   10.3854  165.5108

    6.2888  -19.7145    8.7638  167.4628

Yall =

    8.5444  -16.7598    9.8265  166.9943

    7.0834  -18.9879    9.7658  166.1934

    3.6379  -18.0983    9.1599  164.3775

    4.3880  -16.7007   10.3927  166.5148

    4.9292  -16.7954    8.8020  165.7917