MATLAB Answers

Izzy
0

How to get this code running for multinomial

Asked by Izzy
on 14 Jun 2018 at 22:11
Latest activity Answered by Walter Roberson
on 15 Jun 2018 at 6:05
      function f = ll_ZXL(x)
      S1 = [7 2 3 4 7 7 2 6 2 3 8 9 4 2 1 2 4 4 6 3 9 5];
      S2 = [6 12 6 10 2 4 14 8 16 4 4 2 8 10 10 6 2 4 2 10 2 2];
      F1 = [0 3 12 3 6 3 3 0 0 15 0 0 3 9 9 15 12 9 9 3 0 12];
      F2 = [0 4 0 0 0 0 0 0 0 0 0 0 4 0 4 0 4 4 0 4 0 0];
      nj = [10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10];
      s1 = 0;
      for i = 1:length(S1)
          s1 = s1 + S1(i)*log(x1(i))
      end
      for i = 1:length(S2)
          s1 = s1 + S2(i)*log(x2(i))
      end    
      for i = 1:length(F1)
          s1 = s1 + F1(i)*log(x3(i)) 
      end
      for i = 1:length(F2)
          s1 = s1 + F2(i)*log(x4(i))
      end
      full_sys_output = [3.695788846	5.005193162	4.888786299	4.656239685	5.001524662	4.318879291	4.019492864	3.571165404	3.213860691	3.293049783	4.608615026	3.298121834	4.59052974	4.834971645	6.268156238;
      4.398149409	5.229691373	7.470998131	7.210096038	4.303732461	4.210739644	4.286758278	6.24129195	4.209645611	5.65015171	6.906871609	6.806636386	5.865834622	4.032708516	6.107925457;
      3.760878278	3.284322706	3.660353991	3.828069763	3.582677723	3.976187001	3.733231258	3.195025977	3.403973018	3.446696879	3.222038638	3.024519672	3.566692183	3.016443419	3.757123465;
      2.012861929	2.12873155	1.80779512	1.978172956	2.464945868	1.683127278	1.894733877	2.209921266	2.439077229	2.276070693	1.699735443	1.771360366	2.152983668	1.847740816	1.983080134;
      1.723570081	2.309313186	2.220445389	2.397980796	2.002651622	2.425701234	1.803045452	2.048072691	2.273258285	2.04644325	2.068571042	1.525734283	2.46946355	2.170368503	2.179242482;
      2.518491511	3.211078273	2.857197726	2.823173374	3.477468343	2.986600175	2.590407015	3.323372406	2.531154286	3.098936612	3.424046471	2.870120431	3.272445343	2.776962941	3.42768197;
      1.6047521	1.919790873	2.495299909	2.386878654	2.098953413	2.11427893	1.618709028	2.12222083	2.042327031	2.368536273	1.794653915	2.268644998	1.915704595	2.484637799	1.978874302];
      n = length(full_sys_output);
      Z = log(full_sys_output/60);
      % length
      L = [0.074621236 0.050946986 0.071401538	0.080871238	0.071401538	0.050946986	0.031060616	0.04356062	0.072727296	0.0378788	0.04166668	0.072727296	0.070833356	0.038446982	0.041098498	0.07386366	0.07386366	0.052272744	0.07386366	0.051515168	0.075189418	0.075189418];
      %% 1
      EY1 = 0;
      VARY1 = 0;
      for i = 1:length(S1)
          if i == 1 || i == 2 || i == 3 || i == 4 || i == 5
              EY1 = EY1 + L(i)/5 - 4/45*L(i) * x(i);
              VARY = VARY + (4/45*L(i))^2 * x(i) * (1-x(i));
          end
      end
      EY2 = 0;
      VARY2 = 0;
      for i = 1:length(S2)
          if i == 1 || i == 2 || i == 3 || i == 4 || i == 5
              EY2 = EY2 + L(i)/5 - 4/45*L(i) * x(i);
              VARY2 = VARY2 + (4/45*L(i))^2 * x(i) * (1-x(i));
          end
      end
      for i = 1:length(F1)
          if i == 1 || i == 2 || i == 3 || i == 4 || i == 5
              EY = EY + L(i)/5 - 4/45*L(i) * x(i);
              VARY = VARY + (4/45*L(i))^2 * x(i) * (1-x(i));
          end
      end
      for i = 1:length(F2)
          if i == 1 || i == 2 || i == 3 || i == 4 || i == 5
              EY = EY + L(i)/5 - 4/45*L(i) * x(i);
              VARY = VARY + (4/45*L(i))^2 * x(i) * (1-x(i));
          end
      end
      h1 = log(EY) - 1/2*log(VARY/EY^2 + 1);
      h2 = log(VARY/EY^2 + 1);
      % SUM
      SUM = sum((Z(1,:)-h1).^2);
      f1 = - n/2*log(2*pi) - n/2 * log(h2) - 1/(2*h2)*SUM;
      %% 2
      EY = 0;
      VARY = 0;
      for i = 1:length(S1)
          if i == 6 || i == 7 || i == 8 || i == 9 || i == 10 || i == 11 || i == 12
              EY = EY + L(i)/5 - 4/45*L(i) * x(i);
              VARY = VARY + (4/45*L(i))^2 * x(i) * (1-x(i));
          end
      end
      for i = 1:length(S2)
          if i == 6 || i == 7 || i == 8 || i == 9 || i == 10 || i == 11 || i == 12
              EY = EY + L(i)/5 - 4/45*L(i) * x(i);
              VARY = VARY + (4/45*L(i))^2 * x(i) * (1-x(i));
          end
      end
      for i = 1:length(F1)
          if i == 6 || i == 7 || i == 8 || i == 9 || i == 10 || i == 11 || i == 12
              EY = EY + L(i)/5 - 4/45*L(i) * x(i);
              VARY = VARY + (4/45*L(i))^2 * x(i) * (1-x(i));
          end
      end
      for i = 1:length(F2)
          if i == 6 || i == 7 || i == 8 || i == 9 || i == 10 || i == 11 || i == 12
              EY = EY + L(i)/5 - 4/45*L(i) * x(i);
              VARY = VARY + (4/45*L(i))^2 * x(i) * (1-x(i));
          end
      end
      h1 = log(EY) - 1/2*log(VARY/EY^2 + 1);
      h2 = log(VARY/EY^2 + 1);
      % SUM
      SUM = sum((Z(2,:)-h1).^2);
      f2 = - n/2*log(2*pi) - n/2 * log(h2) - 1/(2*h2)*SUM;
      %% 3
      EY = 0;
      VARY = 0;
      for i = 1:length(S1)
          if i == 13 || i == 14 || i == 15
              EY = EY + L(i)/5 - 4/45*L(i) * x(i);
              VARY = VARY + (4/45*L(i))^2 * x(i) * (1-x(i));
          end
      end
      for i = 1:length(S2)
          if i == 13 || i == 14 || i == 15
              EY = EY + L(i)/5 - 4/45*L(i) * x(i);
              VARY = VARY + (4/45*L(i))^2 * x(i) * (1-x(i));
          end
      end
      for i = 1:length(F1)
          if i == 13 || i == 14 || i == 15
              EY = EY + L(i)/5 - 4/45*L(i) * x(i);
              VARY = VARY + (4/45*L(i))^2 * x(i) * (1-x(i));
          end
      end
      for i = 1:length(F2)
          if i == 13 || i == 14 || i == 15
              EY = EY + L(i)/5 - 4/45*L(i) * x(i);
              VARY = VARY + (4/45*L(i))^2 * x(i) * (1-x(i));
          end
      end
      h1 = log(EY) - 1/2*log(VARY/EY^2 + 1);
      h2 = log(VARY/EY^2 + 1);
      % SUM
      SUM = sum((Z(3,:)-h1).^2);
      f3 = - n/2*log(2*pi) - n/2 * log(h2) - 1/(2*h2)*SUM;
      %% 4
      EY = 0;
      VARY = 0;
      for i = 1:length(S1)
          if i == 21 || i == 20
              EY = EY + L(i)/5 - 4/45*L(i) * x(i);
              VARY = VARY + (4/45*L(i))^2 * x(i) * (1-x(i));
          end
      end
      for i = 1:length(S2)
          if i == 21 || i == 20
              EY = EY + L(i)/5 - 4/45*L(i) * x(i);
              VARY = VARY + (4/45*L(i))^2 * x(i) * (1-x(i));
          end
      end
      for i = 1:length(F1)
          if i == 21 || i == 20
              EY = EY + L(i)/5 - 4/45*L(i) * x(i);
              VARY = VARY + (4/45*L(i))^2 * x(i) * (1-x(i));
          end
      end
      for i = 1:length(F2)
          if i == 21 || i == 20
              EY = EY + L(i)/5 - 4/45*L(i) * x(i);
              VARY = VARY + (4/45*L(i))^2 * x(i) * (1-x(i));
          end
      end
      h1 = log(EY) - 1/2*log(VARY/EY^2 + 1);
      h2 = log(VARY/EY^2 + 1);
      % SUM
      SUM = sum((Z(4,:)-h1).^2);
      f4 = - n/2*log(2*pi) - n/2 * log(h2) - 1/(2*h2)*SUM;
      %% 5
      EY = 0;
      VARY = 0;
      for i = 1:length(S1)
          if i == 18 || i == 16
              EY = EY + L(i)/5 - 4/45*L(i) * x(i);
              VARY = VARY + (4/45*L(i))^2 * x(i) * (1-x(i));
          end
      end
      for i = 1:length(S2)
          if i == 18 || i == 16
              EY = EY + L(i)/5 - 4/45*L(i) * x(i);
              VARY = VARY + (4/45*L(i))^2 * x(i) * (1-x(i));
          end
      end
      for i = 1:length(F1)
          if i == 18 || i == 16
              EY = EY + L(i)/5 - 4/45*L(i) * x(i);
              VARY = VARY + (4/45*L(i))^2 * x(i) * (1-x(i));
          end
      end
      for i = 1:length(F2)
          if i == 18 || i == 16
              EY = EY + L(i)/5 - 4/45*L(i) * x(i);
              VARY = VARY + (4/45*L(i))^2 * x(i) * (1-x(i));
          end
      end
      h1 = log(EY) - 1/2*log(VARY/EY^2 + 1);
      h2 = log(VARY/EY^2 + 1);
      % SUM
      SUM = sum((Z(5,:)-h1).^2);
      f5 = - n/2*log(2*pi) - n/2 * log(h2) - 1/(2*h2)*SUM;
      %% 6
      EY = 0;
      VARY = 0;
      for i = 1:length(S1)
          if i == 22 || i == 11 || i == 12
              EY = EY + L(i)/5 - 4/45*L(i) * x(i);
              VARY = VARY + (4/45*L(i))^2 * x(i) * (1-x(i));
          end
      end
      for i = 1:length(S2)
          if i == 22 || i == 11 || i == 12
              EY = EY + L(i)/5 - 4/45*L(i) * x(i);
              VARY = VARY + (4/45*L(i))^2 * x(i) * (1-x(i));
          end
      end
      for i = 1:length(F1)
          if i == 22 || i == 11 || i == 12
              EY = EY + L(i)/5 - 4/45*L(i) * x(i);
              VARY = VARY + (4/45*L(i))^2 * x(i) * (1-x(i));
          end
      end
      for i = 1:length(F2)
          if i == 22 || i == 11 || i == 12
              EY = EY + L(i)/5 - 4/45*L(i) * x(i);
              VARY = VARY + (4/45*L(i))^2 * x(i) * (1-x(i));
          end
      end
      h1 = log(EY) - 1/2*log(VARY/EY^2 + 1);
      h2 = log(VARY/EY^2 + 1);
      % SUM
      SUM = sum((Z(6,:)-h1).^2);
      f6 = - n/2*log(2*pi) - n/2 * log(h2) - 1/(2*h2)*SUM;
      disp(h2);
      %% 7
      numberofroutes = 7;
      EY = 0;
      EYarray = zeros(numberofroutes);
      VARYarray = zeros(numberofroutes);
      VARY = 0;
      Farray = zeros(1,7);
      for i = 1:length(S1)
          if  i == 1 || i == 2 || i == 3 || i == 4 || i == 5   
              EYarray(1) = EYarray(1) + L(i)/5 - 4/45*L(i) * x(i);
              VARYarray(1) = VARYarray(1) + (4/45*L(i))^2 * x(i) * (1-x(i));
          elseif i == 6 || i == 7 || i == 8 || i == 9 || i == 10 || i == 11 || i == 12
              EYarray(2) = EYarray(2) + L(i)/5 - 4/45*L(i) * x(i);
              VARYarray(2) = VARYarray(2) + (4/45*L(i))^2 * x(i) * (1-x(i)); 
          elseif i == 13 || i == 14 || i == 15
              EYarray(3) = EYarray(3) + L(i)/5 - 4/45*L(i) * x(i);
              VARYarray(3) = VARYarray(3) + (4/45*L(i))^2 * x(i) * (1-x(i)); 
          elseif i == 21 || i == 20
              EYarray(4) = EYarray(4) + L(i)/5 - 4/45*L(i) * x(i);
              VARYarray(4) = VARYarray(4) + (4/45*L(i))^2 * x(i) * (1-x(i)); 
          elseif i == 18 || i == 16
              EYarray(5) = EYarray(5) + L(i)/5 - 4/45*L(i) * x(i);
              VARYarray(5) = VARYarray(5) + (4/45*L(i))^2 * x(i) * (1-x(i)); 
          elseif i == 22 || i == 11 || i == 12
               EYarray(6) = EYarray(6) + L(i)/5 - 4/45*L(i) * x(i);
              VARYarray(6) = VARYarray(6) + (4/45*L(i))^2 * x(i) * (1-x(i)); 
          elseif i == 17 || i == 19
              EYarray(7) = EYarray(7) + L(i)/5 - 4/45*L(i) * x(i);
              VARYarray(7) = VARYarray(7) + (4/45*L(i))^2 * x(i) * (1-x(i));
          else 
          end
      end
      for i = 1:length(S2)
          if  i == 1 || i == 2 || i == 3 || i == 4 || i == 5   
              EYarray(1) = EYarray(1) + L(i)/5 - 4/45*L(i) * x(i);
              VARYarray(1) = VARYarray(1) + (4/45*L(i))^2 * x(i) * (1-x(i));
          elseif i == 6 || i == 7 || i == 8 || i == 9 || i == 10 || i == 11 || i == 12
              EYarray(2) = EYarray(2) + L(i)/5 - 4/45*L(i) * x(i);
              VARYarray(2) = VARYarray(2) + (4/45*L(i))^2 * x(i) * (1-x(i)); 
          elseif i == 13 || i == 14 || i == 15
              EYarray(3) = EYarray(3) + L(i)/5 - 4/45*L(i) * x(i);
              VARYarray(3) = VARYarray(3) + (4/45*L(i))^2 * x(i) * (1-x(i)); 
          elseif i == 21 || i == 20
              EYarray(4) = EYarray(4) + L(i)/5 - 4/45*L(i) * x(i);
              VARYarray(4) = VARYarray(4) + (4/45*L(i))^2 * x(i) * (1-x(i)); 
          elseif i == 18 || i == 16
              EYarray(5) = EYarray(5) + L(i)/5 - 4/45*L(i) * x(i);
              VARYarray(5) = VARYarray(5) + (4/45*L(i))^2 * x(i) * (1-x(i)); 
          elseif i == 22 || i == 11 || i == 12
               EYarray(6) = EYarray(6) + L(i)/5 - 4/45*L(i) * x(i);
              VARYarray(6) = VARYarray(6) + (4/45*L(i))^2 * x(i) * (1-x(i)); 
          elseif i == 17 || i == 19
              EYarray(7) = EYarray(7) + L(i)/5 - 4/45*L(i) * x(i);
              VARYarray(7) = VARYarray(7) + (4/45*L(i))^2 * x(i) * (1-x(i));
          else 
          end
      end
      for i = 1:length(F1)
          if  i == 1 || i == 2 || i == 3 || i == 4 || i == 5   
              EYarray(1) = EYarray(1) + L(i)/5 - 4/45*L(i) * x(i);
              VARYarray(1) = VARYarray(1) + (4/45*L(i))^2 * x(i) * (1-x(i));
          elseif i == 6 || i == 7 || i == 8 || i == 9 || i == 10 || i == 11 || i == 12
              EYarray(2) = EYarray(2) + L(i)/5 - 4/45*L(i) * x(i);
              VARYarray(2) = VARYarray(2) + (4/45*L(i))^2 * x(i) * (1-x(i)); 
          elseif i == 13 || i == 14 || i == 15
              EYarray(3) = EYarray(3) + L(i)/5 - 4/45*L(i) * x(i);
              VARYarray(3) = VARYarray(3) + (4/45*L(i))^2 * x(i) * (1-x(i)); 
          elseif i == 21 || i == 20
              EYarray(4) = EYarray(4) + L(i)/5 - 4/45*L(i) * x(i);
              VARYarray(4) = VARYarray(4) + (4/45*L(i))^2 * x(i) * (1-x(i)); 
          elseif i == 18 || i == 16
              EYarray(5) = EYarray(5) + L(i)/5 - 4/45*L(i) * x(i);
              VARYarray(5) = VARYarray(5) + (4/45*L(i))^2 * x(i) * (1-x(i)); 
          elseif i == 22 || i == 11 || i == 12
               EYarray(6) = EYarray(6) + L(i)/5 - 4/45*L(i) * x(i);
              VARYarray(6) = VARYarray(6) + (4/45*L(i))^2 * x(i) * (1-x(i)); 
          elseif i == 17 || i == 19
              EYarray(7) = EYarray(7) + L(i)/5 - 4/45*L(i) * x(i);
              VARYarray(7) = VARYarray(7) + (4/45*L(i))^2 * x(i) * (1-x(i));
          else 
          end
      end
      for i = 1:length(F2)
          if  i == 1 || i == 2 || i == 3 || i == 4 || i == 5   
              EYarray(1) = EYarray(1) + L(i)/5 - 4/45*L(i) * x(i);
              VARYarray(1) = VARYarray(1) + (4/45*L(i))^2 * x(i) * (1-x(i));
          elseif i == 6 || i == 7 || i == 8 || i == 9 || i == 10 || i == 11 || i == 12
              EYarray(2) = EYarray(2) + L(i)/5 - 4/45*L(i) * x(i);
              VARYarray(2) = VARYarray(2) + (4/45*L(i))^2 * x(i) * (1-x(i)); 
          elseif i == 13 || i == 14 || i == 15
              EYarray(3) = EYarray(3) + L(i)/5 - 4/45*L(i) * x(i);
              VARYarray(3) = VARYarray(3) + (4/45*L(i))^2 * x(i) * (1-x(i)); 
          elseif i == 21 || i == 20
              EYarray(4) = EYarray(4) + L(i)/5 - 4/45*L(i) * x(i);
              VARYarray(4) = VARYarray(4) + (4/45*L(i))^2 * x(i) * (1-x(i)); 
          elseif i == 18 || i == 16
              EYarray(5) = EYarray(5) + L(i)/5 - 4/45*L(i) * x(i);
              VARYarray(5) = VARYarray(5) + (4/45*L(i))^2 * x(i) * (1-x(i)); 
          elseif i == 22 || i == 11 || i == 12
               EYarray(6) = EYarray(6) + L(i)/5 - 4/45*L(i) * x(i);
              VARYarray(6) = VARYarray(6) + (4/45*L(i))^2 * x(i) * (1-x(i)); 
          elseif i == 17 || i == 19
              EYarray(7) = EYarray(7) + L(i)/5 - 4/45*L(i) * x(i);
              VARYarray(7) = VARYarray(7) + (4/45*L(i))^2 * x(i) * (1-x(i));
          else 
          end
      end
      for g = 1:numberofroutes
          h1 = log(EYarray(g)) - 1/2*log(VARYarray(g)/EYarray(g)^2 + 1);
          h2 = log(VARYarray/EYarray(g)^2 + 1);
          SUM = sum((Z(g,:)-h1).^2);
          Farray(1,g) = - n/2*log(2*pi) - n/2 * log(h2(g,1)) - 1/(2*h2(g,1))*SUM;
      end
      %% f
      f = sum(Farray)+ s1+s2+f1+f2;
      f = -f;

Here is the code to run

clc
clear all
options = optimoptions('fmincon','Algorithm','sqp');  
% # x 1
S1 = [7 2 3 4 7 7 2 6 2 3 8 9 4 2 1 2 4 4 6 3 9 5];
% # x 2
S2 = [6 12 6 10 2 4 14 8 16 4 4 2 8 10 10 6 2 4 2 10 2 2];
% # x 3
F1 = [0 3 12 3 6 3 3 0 0 15 0 0 3 9 9 15 12 9 9 3 0 12];
% # x 4
F2 = [0 4 0 0 0 0 0 0 0 0 0 0 4 0 4 0 4 4 0 4 0 0];
% # 
nj = [9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9];
L = [0.074621236	0.050946986	0.071401538	0.080871238	0.071401538	0.050946986	0.031060616	0.04356062	0.072727296	0.0378788	0.04166668	0.072727296	0.070833356	0.038446982	0.041098498	0.07386366	0.07386366	0.052272744	0.07386366	0.051515168	0.075189418	0.075189418];
x1 = S1./nj;
x2 = S2./nj;
x3 = F1./nj;
x4 = F2./nj;
x11 = x1/1.0001;
x21 = x2/1.0001;
x31 = x3/1.0001;
x41 = x4/1.0001;
lb1 = zeros(1,length(S1));
ub1 = ones(1,length(S1));
lb2 = zeros(1,length(S2));
ub2 = ones(1,length(S2));
lb3 = zeros(1,length(F1));
ub3 = ones(1,length(F1));
lb4 = zeros(1,length(F2));
ub4 = ones(1,length(F2));
b1 = ones(length(x1),1);
b2 = ones(length(x2),1);
b3 = ones(length(x3),1);
b4 = ones(length(x4),1);
y1 = fmincon(@ll_ZXL,x11,[],[],[],[],lb1,ub1,[],options)
y2 = fmincon(@ll_ZXL,x21,[],[],[],[],lb2,ub2,[],options)
y3 = fmincon(@ll_ZXL,x31,[],[],[],[],lb3,ub3,[],options)
y4 = fmincon(@ll_ZXL,x41,[],[],[],[],lb4,ub4,[],options)
loglike_y1 = -ll_ZXL(y1)
loglike_x1 = -ll_ZXL(x11)
loglike_y2 = -ll_ZXL(y2)
loglike_x2 = -ll_ZXL(x21)
loglike_y3 = -ll_ZXL(y3)
loglike_x3 = -ll_ZXL(x31)
loglike_y4 = -ll_ZXL(y4)
loglike_x4 = -ll_ZXL(x41)
error1 = (y1-x1)./x1*100;
error2 = (y2-x2)./x2*100;
error3 = (y3-x3)./x3*100;
error4 = (y4-x4)./x4*100;
T1 = table([1:length(S1)]', L', x1',y1',error1','VariableNames',{'Link_Num';'Length';'Sample_Mean';'MLE';'Percentage_Error'})
T2 = table([1:length(S2)]', L', x2',y2',error2','VariableNames',{'Link_Num';'Length';'Sample_Mean';'MLE';'Percentage_Error'})
T3 = table([1:length(F1)]', L', x3',y3',error3','VariableNames',{'Link_Num';'Length';'Sample_Mean';'MLE';'Percentage_Error'})
T4 = table([1:length(F2)]', L', x4',y4',error4','VariableNames',{'Link_Num';'Length';'Sample_Mean';'MLE';'Percentage_Error'})

Here is the error message:

Undefined function or variable 'x1'.
Error in ll_ZXL (line 24)
    s1 = s1 + S1(i)*log(x1(i)) %+ S2(i)*log(x2(i)) + F1(i)*log(x3(i)) + F2(i)*log(x4(i))%;%(nj(i)-S1(i))*log(1-x(i));
Error in fmincon (line 534)
      initVals.f = feval(funfcn{3},X,varargin{:});
Error in MLE_Test (line 45)
y1 = fmincon(@ll_ZXL,x11,[],[],[],[],lb1,ub1,[],options)

  0 Comments

Sign in to comment.

1 Answer

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today