How to get this code running for multinomial

1 view (last 30 days)
Izzy
Izzy on 14 Jun 2018
Commented: Rena Berman on 27 Jul 2018
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)

Sign in to answer this question.

Answers (1)

Walter Roberson
Walter Roberson on 15 Jun 2018

Categories

Find more on Pulsed Waveforms in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!