$offsymxref offsymlist **The DUAL QC case (note: no weights) *Based on "qcmc100n.gms" * Note, creating the data based on Bill's method (e.g., "make" file in ARE) * (Nov, 1994) with 3 categories. set i index /1*337/ j index /1*3/ j0 index /1*3/ j1 index /1*3/ j2 index /1*3/ k index /1*4/ k2 index /1*4/ h index /1*3/ h0 index /1*3/ h1 index /1*3/ l index /1*3/ w index /1*3/ d index /1*3/ n index /1*5/ g index /1*1/ ; set i1(n) g1(g) ; alias (k1,k); alias (k3,k); i1(n)=yes; *i1("1")=no; parameter x(i,k) observed aggregate outcomes (inputs outputs etc.) xdata(i,k) keep data y and x ind(i) index alpha(j,k) parameter space for each RHS var and observation alphag(n,j,k) parameter space for each RHS var and observation psi(i,n) w0(i,n,h) ei0(i,n) alphaeg(h,n) parameter space for each errors alphab0(j0,k) parameter space for RHS vars correspond to y=0 alphab1(j1,k) parameter space for RHS var correspond to y=1 alphab2(j2,k) parameter space for RHS var correspond to y=2 qb(n,k,j) prior probabilities for betas qb0(k,j0) prior probabilities for beta0 qb1(k,j1) prior probabilities for beta1 qd1(n,d) prior prob. for delta1 qd2(n,d) prior prob. for delta2 qe(i,h) prior probabilities for errors qe0(i,h0) prior probabilities for errors of beta0 qee1(i,h1) prior probabilities for errors of beta1 qw(i,w) prior probabilities for x'beta qe1(i,l) prior probabilities for errors error(i) estimated errors error1(i) estimated errors for betas a no. of parameters (length of alpha = max of j) beta(g,n,k) estimated espected value of oefficients omega(k) partition function for k var delta Weight factor for objective functional sume(i) normalization for ME dist denom(i) denominator p1(i) temporary for data generation p2(i) temporary for data generation p3(i) temporary for data generation limit(i) temporary for data generation u(i) temporary for data generation p0(k,j) estimated probabilities recp(i,n) recovered p's recw(i,n,h) recovered w's omeg(i) partition fn for i psi(i,n) partition fn for w tablep(i,n,n) table for predictin results tablepred(n,n) predicted vs. exact cases param(k,j) parameter (beta) space paramw(i,w) parameter (x'beta) space entropy Entropy value (value of objective fn) MSE(g) mean sq. error entp(g) normalized ent for p's entpi(g,i) normalized ent for each i entpe(g) normalized ent for errors aventp average ent for p's varentp var of ent for p's aventpi(i) average ent per i varentpi(i) var of ent per i aventpe average ent per errors varentpe var ent per errors avmse average MSE varmse var MSE dist(g,n,k) distance between correct and estimated betas var(k) variance of coefficients ent(n,k) normalized entropy of coefficients ypred(i) predicted LHS var pred(i) predicting var ypred1(i) predicted minus correct y sumyp(g) total number of misses minmiss min number of misses maxmiss max number of misses predprop predicted proportion of ones exactprop correct proportion of ones xbeta(n,i) estimated x'beta exbeta(n,i) estimated exp(x'beta) xbeta1(n,i) estimated recovered probabilities xbetam(i) max xbeta xbetae(n,i) estimated x'beta+estimated errors SEbeta(i) square error for each observation SEbeta1 square error for all observations avbeta(n,k) average values of betas varbeta(n,k) var of betas totvar variance of betas avtablep(n,n) average value of prediction vs correct vartablep(n,n) var of value of prediction vs correct avmiss averages number of misses varmiss var of misses JJ number of parameters in param space (size of j) JJ0 number of parameters in param space (size of j0) JJ1 number of parameters in param space (size of j1) II number of observations (size of i) HH number of param space for errors (size of h) HH0 number of param space for errors (size of h) HH1 number of param space for errors (size of h) LL number of param space for consistency errors (size of l) WW number of points for w index DD number of points for d index NN total number of categories (i.e., index n) GG number of experiments KK **Below params for calculating Cov and Info matrices omega1(g,i) IME1(g,n) IME2(g,n) IME3(g,n) IME4(g,n) IML1(g,n) IML2(g,n) IEE(g,n) INFO(g,n,k,k) INFO1(g,n,k,k) TINFO(g,n,k,k) avIME1(n) avIME2(n) avIME3(n) avIME4(n) avIML1(n) avIML2(n) xsq(i) xsqk(k,k) ; KK=4; NN=5; JJ=3; JJ0=3; JJ1=3; II=337; GG=1; HH=3; HH0=3; HH1=3; LL=3; WW=3; DD=3; a=12; delta=0.5; *Generating the data *First step: generating the x's parameter rho1 determines correlation between x's t1(i) t2(i) t3(i) x(i,k) the multinomial data xx(i,k2) the multinomial data beta0(n,k) the correct betas e(i,n) the errors s(i,n) value of "continuous" LHS var b0(i) the correct LHS variable (category) y(i,n) the 0 1 specification b(i) transposing b0 to creat the y's ; rho1 = .85; table beta0(n,k) the correct betas * 1 2 3 4 *1 0 0 0 0 *2 -1 1 2 -1 *3 1 -1 -2 1 **Just for ommitting var exp 1 2 3 1 0 0 0 2 -1 1 2 3 1 -1 -2 ; options decimals=7; display beta0; parameter alphae(h) param space for errors / 1 0.0 2 0.0 3 0.0 *1 -.1 *2 0 *3 .1 *1 -.1 *2 -.0666 *3 -.0333 *4 0 *5 .0333 *6 .0666 *7 .1 /; parameter alphae1(h) param space for beta0 errors / *1 0.0 *2 0.0 *3 0.0 1 -1000.0 2 0.0 3 1000.0 /; *Uniform priors *qb(n,k,j)=1/JJ; *qe(i,h)=1/HH; variables e1(i,n) error coefficients lm(n,k) Lagrange multipliers ; positive variables p(i,n) probabilities of i choosing j pe(i,n,h) probabilities for errors ; variables obj objective function ; equations objective objective function objdual dual objective function normalize(n,k) Normalizing first category addp(i) adding up for p(in) adde(i,n) adding up for errors e1k0(i,n) e1k coefficients e1k1(i,n) e1k coefficients e1k2(i,n) e1k coefficients e1k3(i,n) e1k coefficients e1k4(i,n) e1k coefficients consist(n,k) consistency for n and k ; ********** *The ME objective for Soofgi with BANG objective.. obj =e= -sum(i,sum(n,p(i,n)*log((p(i,n)+1.e-8)))) - sum(i,sum(n,sum(h,pe(i,n,h)*log((pe(i,n,h)+1.e-8))))); ********** ********** *The ME objective for Soofgi with BANG and weight factor for objective *objective.. obj =e=(1-delta)*(-sum(i,sum(n,p(i,n)*log((p(i,n)+1.e-4))))) -delta*(sum(i,sum(n,sum(h,pe(i,n,h)*log((pe(i,n,h)+1.e-4)))))); ********** ********* *Dual Objective With Noise objdual.. obj =e= - sum(i,sum(n,sum(k,y(i,n)*x(i,k)*lm(n,k)))) + sum(i,log(sum(n,exp(sum(k,x(i,k)*lm(n,k)))))) + sum(i,sum(n,log(sum(h,exp(sum(k,x(i,k)*alphae(h)*lm(n,k))))))); ********* ********* *Dual Objective (PURE) *objdual.. obj =e= - sum(i,sum(n,sum(k,y(i,n)*x(i,k)*lm(n,k)))) * + sum(i,log(sum(n,exp(sum(k,x(i,k)*lm(n,k)))))); ********* normalize(n,k).. lm("1",k) =e=0; addp(i).. sum(n,p(i,n)) =e= 1; adde(i,n).. sum(h,pe(i,n,h)) =e= 1; e1k0(i,"1").. sum(h,alphae1(h)*pe(i,"1",h)) =e= e1(i,"1"); e1k1(i,"2").. sum(h,alphae(h)*pe(i,"2",h)) =e= e1(i,"2"); e1k2(i,"3").. sum(h,alphae(h)*pe(i,"3",h)) =e= e1(i,"3"); e1k3(i,"4").. sum(h,alphae(h)*pe(i,"4",h)) =e= e1(i,"4"); e1k4(i,"5").. sum(h,alphae(h)*pe(i,"5",h)) =e= e1(i,"5"); *The Soofgi with BANG case consist(n,k).. sum(i,p(i,n)*x(i,k)) + sum(i,e1(i,n)*x(i,k)) =e= sum(i,y(i,n)*x(i,k)); **Note "l" means less or equal **Note "g" means greater or equal *option nlp=conopt ; *Note: remove zero when alpha space is equal for all model base test /objective addp e1k0 e1k1 e1k2 e1k3 e1k4 adde consist /; model dual /objdual normalize /; *model ols leas squares /obj1/; * conopt and minos are two non-linear optimization packages * minos is used by default. options limrow=0,limcol=0; options domlim=30; options iterlim=14000; options reslim=900000; *options work=12000000; * conopt optimizer option is below options solprint=off, nlp=conopt; *options nlp=conopt; *rho1 = .85; *t1(i) = normal(0,1); *t2(i) = normal(0,1); *t3(i) = normal(0,1); table xx(i,k2) koop-poirier data 1 2 3 4 1 1 3.000000 11.00000 1.000000 2 1 14.00000 12.00000 1.000000 3 1 44.00000 12.00000 1.000000 4 1 18.00000 12.00000 1.000000 5 1 24.00000 14.00000 0. 6 1 38.00000 13.00000 1.000000 7 1 8.000000 14.00000 0. 8 1 19.00000 14.00000 1.000000 9 1 8.000000 12.00000 1.000000 10 1 3.000000 12.00000 1.000000 11 1 6.000000 12.00000 1.000000 12 1 40.00000 11.00000 1.000000 13 1 2.000000 12.00000 1.000000 14 1 22.00000 12.00000 1.000000 15 1 4.000000 12.00000 1.000000 16 1 22.00000 9.000000 1.000000 17 1 39.00000 6.000000 1.000000 18 1 3.000000 12.00000 1.000000 19 1 3.000000 10.00000 1.000000 20 1 7.000000 9.000000 1.000000 21 1 27.00000 15.00000 0. 22 1 52.00000 6.000000 0. 23 1 3.000000 12.00000 0. 24 1 35.00000 12.00000 1.000000 25 1 8.000000 14.00000 1.000000 26 1 34.00000 11.00000 1.000000 27 1 14.00000 14.00000 1.000000 28 1 9.000000 9.000000 1.000000 29 1 9.000000 14.00000 1.000000 30 1 5.000000 13.00000 1.000000 31 1 52.00000 14.00000 1.000000 32 1 10.00000 12.00000 1.000000 33 1 25.00000 12.00000 1.000000 34 1 32.00000 12.00000 1.000000 35 1 29.00000 12.00000 1.000000 36 1 4.000000 12.00000 1.000000 37 1 11.00000 15.00000 1.000000 38 1 11.00000 16.00000 1.000000 39 1 19.00000 16.00000 1.000000 40 1 18.00000 12.00000 1.000000 41 1 25.00000 10.00000 1.000000 42 1 8.000000 9.000000 1.000000 43 1 6.000000 10.00000 1.000000 44 1 2.000000 12.00000 1.000000 45 1 23.00000 8.000000 1.000000 46 1 8.000000 12.00000 1.000000 47 1 29.00000 3.000000 1.000000 48 1 30.00000 13.00000 1.000000 49 1 17.00000 12.00000 1.000000 50 1 9.000000 12.00000 1.000000 51 1 11.00000 15.00000 1.000000 52 1 9.000000 14.00000 1.000000 53 1 17.00000 10.00000 1.000000 54 1 4.000000 11.00000 1.000000 55 1 30.00000 9.000000 1.000000 56 1 22.00000 16.00000 0. 57 1 29.00000 10.00000 1.000000 58 1 6.000000 10.00000 1.000000 59 1 11.00000 16.00000 1.000000 60 1 5.000000 12.00000 1.000000 61 1 12.00000 12.00000 1.000000 62 1 26.00000 8.000000 1.000000 63 1 35.00000 8.000000 1.000000 64 1 17.00000 12.00000 1.000000 65 1 46.00000 6.000000 1.000000 66 1 6.000000 11.00000 1.000000 67 1 37.00000 6.000000 1.000000 68 1 32.00000 11.00000 1.000000 69 1 43.00000 8.000000 1.000000 70 1 4.000000 12.00000 1.000000 71 1 46.00000 6.000000 1.000000 72 1 51.00000 13.00000 1.000000 73 1 39.00000 10.00000 1.000000 74 1 37.00000 12.00000 1.000000 75 1 10.00000 12.00000 1.000000 76 1 4.000000 12.00000 1.000000 77 1 4.000000 12.00000 1.000000 78 1 49.00000 14.00000 0. 79 1 32.00000 12.00000 1.000000 80 1 9.000000 12.00000 1.000000 81 1 9.000000 12.00000 1.000000 82 1 8.000000 13.00000 1.000000 83 1 5.000000 12.00000 1.000000 84 1 34.00000 9.000000 1.000000 85 1 19.00000 8.000000 1.000000 86 1 41.00000 7.000000 0. 87 1 37.00000 14.00000 1.000000 88 1 4.000000 9.000000 1.000000 89 1 43.00000 11.00000 1.000000 90 1 14.00000 12.00000 1.000000 91 1 9.000000 12.00000 1.000000 92 1 33.00000 8.000000 1.000000 93 1 15.00000 13.00000 1.000000 94 1 12.00000 12.00000 1.000000 95 1 19.00000 13.00000 1.000000 96 1 23.00000 8.000000 0. 97 1 26.00000 13.00000 1.000000 98 1 13.00000 13.00000 1.000000 99 1 22.00000 12.00000 1.000000 100 1 4.000000 11.00000 1.000000 101 1 22.00000 12.00000 1.000000 102 1 10.00000 11.00000 1.000000 103 1 21.00000 9.000000 1.000000 104 1 38.00000 6.000000 0. 105 1 11.00000 12.00000 1.000000 106 1 47.00000 9.000000 1.000000 107 1 18.00000 13.00000 1.000000 108 1 8.000000 12.00000 1.000000 109 1 13.00000 12.00000 1.000000 110 1 10.00000 12.00000 1.000000 111 1 41.00000 11.00000 1.000000 112 1 49.00000 11.00000 1.000000 113 1 4.000000 13.00000 1.000000 114 1 9.000000 12.00000 1.000000 115 1 33.00000 12.00000 1.000000 116 1 2.000000 12.00000 1.000000 117 1 11.00000 15.00000 1.000000 118 1 56.00000 6.000000 1.000000 119 1 31.00000 13.00000 1.000000 120 1 13.00000 14.00000 1.000000 121 1 33.00000 11.00000 1.000000 122 1 41.00000 12.00000 1.000000 123 1 6.000000 12.00000 1.000000 124 1 21.00000 12.00000 1.000000 125 1 25.00000 13.00000 1.000000 126 1 13.00000 15.00000 1.000000 127 1 2.000000 12.00000 1.000000 128 1 23.00000 12.00000 1.000000 129 1 32.00000 12.00000 1.000000 130 1 46.00000 12.00000 1.000000 131 1 13.00000 12.00000 1.000000 132 1 29.00000 12.00000 1.000000 133 1 30.00000 12.00000 1.000000 134 1 50.00000 10.00000 1.000000 135 1 32.00000 10.00000 1.000000 136 1 29.00000 12.00000 1.000000 137 1 9.000000 16.00000 0. 138 1 49.00000 8.000000 1.000000 139 1 9.000000 14.00000 0. 140 1 41.00000 14.00000 1.000000 141 1 9.000000 12.00000 1.000000 142 1 5.000000 11.00000 1.000000 143 1 17.00000 12.00000 1.000000 144 1 9.000000 11.00000 1.000000 145 1 30.00000 12.00000 1.000000 146 1 29.00000 7.000000 0. 147 1 9.000000 14.00000 1.000000 148 1 37.00000 12.00000 1.000000 149 1 44.00000 7.000000 0. 150 1 22.00000 12.00000 1.000000 151 1 26.00000 12.00000 1.000000 152 1 10.00000 12.00000 1.000000 153 1 33.00000 13.00000 1.000000 154 1 41.00000 8.000000 1.000000 155 1 39.00000 12.00000 1.000000 156 1 29.00000 12.00000 0. 157 1 38.00000 14.00000 1.000000 158 1 12.00000 12.00000 0. 159 1 9.000000 12.00000 0. 160 1 10.00000 14.00000 1.000000 161 1 9.000000 16.00000 0. 162 1 20.00000 12.00000 1.000000 163 1 9.000000 11.00000 1.000000 164 1 41.00000 14.00000 1.000000 165 1 6.000000 14.00000 1.000000 166 1 10.00000 12.00000 1.000000 167 1 11.00000 14.00000 0. 168 1 21.00000 12.00000 1.000000 169 1 20.00000 13.00000 1.000000 170 1 31.00000 14.00000 1.000000 171 1 4.000000 16.00000 1.000000 172 1 12.00000 13.00000 1.000000 173 1 17.00000 14.00000 1.000000 174 1 40.00000 6.000000 1.000000 175 1 53.00000 12.00000 1.000000 176 1 35.00000 14.00000 1.000000 177 1 12.00000 14.00000 1.000000 178 1 13.00000 15.00000 1.000000 179 1 48.00000 8.000000 1.000000 180 1 23.00000 12.00000 1.000000 181 1 11.00000 12.00000 1.000000 182 1 9.000000 12.00000 1.000000 183 1 9.000000 12.00000 1.000000 184 1 4.000000 12.00000 1.000000 185 1 34.00000 16.00000 1.000000 186 1 12.00000 12.00000 1.000000 187 1 21.00000 13.00000 0. 188 1 12.00000 15.00000 1.000000 189 1 17.00000 12.00000 1.000000 190 1 21.00000 12.00000 1.000000 191 1 20.00000 12.00000 1.000000 192 1 35.00000 12.00000 0. 193 1 44.00000 15.00000 1.000000 194 1 6.000000 16.00000 1.000000 195 1 5.000000 14.00000 1.000000 196 1 42.00000 11.00000 1.000000 197 1 34.00000 12.00000 1.000000 198 1 37.00000 16.00000 1.000000 199 1 19.00000 13.00000 1.000000 200 1 32.00000 12.00000 1.000000 201 1 25.00000 12.00000 1.000000 202 1 19.00000 12.00000 1.000000 203 1 50.00000 12.00000 1.000000 204 1 6.000000 12.00000 1.000000 205 1 49.00000 12.00000 1.000000 206 1 3.000000 11.00000 1.000000 207 1 49.00000 18.00000 1.000000 208 1 39.00000 15.00000 1.000000 209 1 20.00000 15.00000 1.000000 210 1 10.00000 12.00000 1.000000 211 1 5.000000 12.00000 1.000000 212 1 10.00000 13.00000 1.000000 213 1 30.00000 16.00000 1.000000 214 1 31.00000 15.00000 1.000000 215 1 9.000000 12.00000 1.000000 216 1 8.000000 12.00000 1.000000 217 1 49.00000 13.00000 1.000000 218 1 11.00000 16.00000 1.000000 219 1 2.000000 12.00000 1.000000 220 1 6.000000 12.00000 1.000000 221 1 12.00000 10.00000 1.000000 222 1 5.000000 17.00000 1.000000 223 1 3.000000 12.00000 1.000000 224 1 6.000000 16.00000 1.000000 225 1 38.00000 8.000000 1.000000 226 1 9.000000 16.00000 1.000000 227 1 10.00000 16.00000 1.000000 228 1 37.00000 14.00000 1.000000 229 1 13.00000 14.00000 1.000000 230 1 11.00000 13.00000 1.000000 231 1 21.00000 16.00000 0. 232 1 4.000000 16.00000 1.000000 233 1 2.000000 16.00000 1.000000 234 1 6.000000 16.00000 1.000000 235 1 13.00000 17.00000 1.000000 236 1 13.00000 16.00000 1.000000 237 1 18.00000 16.00000 1.000000 238 1 44.00000 14.00000 1.000000 239 1 9.000000 16.00000 1.000000 240 1 38.00000 16.00000 1.000000 241 1 25.00000 14.00000 1.000000 242 1 32.00000 13.00000 1.000000 243 1 18.00000 17.00000 1.000000 244 1 22.00000 19.00000 1.000000 245 1 20.00000 12.00000 1.000000 246 1 11.00000 18.00000 1.000000 247 1 9.000000 16.00000 1.000000 248 1 14.00000 14.00000 1.000000 249 1 2.000000 15.00000 1.000000 250 1 8.000000 14.00000 1.000000 251 1 26.00000 12.00000 0. 252 1 10.00000 15.00000 1.000000 253 1 8.000000 12.00000 1.000000 254 1 25.00000 16.00000 1.000000 255 1 8.000000 16.00000 1.000000 256 1 20.00000 14.00000 1.000000 257 1 11.00000 12.00000 1.000000 258 1 23.00000 20.00000 1.000000 259 1 25.00000 19.00000 1.000000 260 1 31.00000 13.00000 1.000000 261 1 2.000000 16.00000 1.000000 262 1 25.00000 14.00000 1.000000 263 1 3.000000 16.00000 1.000000 264 1 16.00000 19.00000 0. 265 1 6.000000 14.00000 1.000000 266 1 17.00000 19.00000 1.000000 267 1 23.00000 15.00000 1.000000 268 1 16.00000 14.00000 1.000000 269 1 19.00000 13.00000 1.000000 270 1 9.000000 14.00000 1.000000 271 1 14.00000 19.00000 1.000000 272 1 45.00000 14.00000 1.000000 273 1 17.00000 16.00000 1.000000 274 1 42.00000 13.00000 1.000000 275 1 15.00000 12.00000 1.000000 276 1 17.00000 12.00000 1.000000 277 1 8.000000 16.00000 1.000000 278 1 15.00000 17.00000 1.000000 279 1 37.00000 18.00000 1.000000 280 1 33.00000 16.00000 1.000000 281 1 16.00000 11.00000 1.000000 282 1 13.00000 18.00000 1.000000 283 1 44.00000 15.00000 1.000000 284 1 28.00000 16.00000 1.000000 285 1 41.00000 11.00000 1.000000 286 1 25.00000 16.00000 1.000000 287 1 66.00000 12.00000 1.000000 288 1 31.00000 12.00000 1.000000 289 1 25.00000 20.00000 1.000000 290 1 17.00000 17.00000 1.000000 291 1 23.00000 16.00000 1.000000 292 1 34.00000 17.00000 1.000000 293 1 21.00000 16.00000 1.000000 294 1 4.000000 15.00000 1.000000 295 1 41.00000 8.000000 1.000000 296 1 8.000000 12.00000 1.000000 297 1 2.000000 19.00000 1.000000 298 1 20.00000 16.00000 1.000000 299 1 38.00000 15.00000 1.000000 300 1 23.00000 16.00000 0. 301 1 31.00000 18.00000 1.000000 302 1 19.00000 20.00000 1.000000 303 1 4.000000 20.00000 1.000000 304 1 29.00000 19.00000 0. 305 1 47.00000 12.00000 0. 306 1 27.00000 20.00000 1.000000 307 1 11.00000 16.00000 1.000000 308 1 24.00000 13.00000 1.000000 309 1 12.00000 14.00000 1.000000 310 1 14.00000 19.00000 1.000000 311 1 22.00000 16.00000 1.000000 312 1 11.00000 18.00000 1.000000 313 1 16.00000 16.00000 1.000000 314 1 18.00000 19.00000 1.000000 315 1 17.00000 13.00000 1.000000 316 1 41.00000 13.00000 1.000000 317 1 6.000000 16.00000 1.000000 318 1 3.000000 16.00000 1.000000 319 1 6.000000 19.00000 1.000000 320 1 36.00000 18.00000 1.000000 321 1 3.000000 12.00000 1.000000 322 1 6.000000 18.00000 1.000000 323 1 6.000000 13.00000 1.000000 324 1 9.000000 17.00000 1.000000 325 1 6.000000 16.00000 1.000000 326 1 14.00000 19.00000 1.000000 327 1 7.000000 19.00000 1.000000 328 1 11.00000 19.00000 0. 329 1 14.00000 14.00000 1.000000 330 1 13.00000 16.00000 1.000000 331 1 24.00000 15.00000 1.000000 332 1 7.000000 18.00000 1.000000 333 1 19.00000 13.00000 1.000000 334 1 43.00000 12.00000 1.000000 335 1 31.00000 12.00000 1.000000 336 1 39.00000 7.000000 1.000000 337 1 12.00000 16.00000 1.000000 ; x(i,"1") = xx(i,"1"); x(i,"2") = xx(i,"2"); x(i,"3") = xx(i,"3"); x(i,"4") = xx(i,"4"); **Note adjust the index k before experimenting with number of vars **Omitting var 1 *x(i,"1") = xx(i,"2"); *x(i,"2") = xx(i,"3"); *x(i,"3") = xx(i,"4"); **Omitting var 2 *x(i,"1") = xx(i,"1"); *x(i,"2") = xx(i,"3"); *x(i,"3") = xx(i,"4"); **Omitting var 3 *x(i,"1") = xx(i,"1"); *x(i,"2") = xx(i,"2"); *x(i,"3") = xx(i,"4"); **Omitting var 4 *x(i,"1") = xx(i,"1"); *x(i,"2") = xx(i,"2"); *x(i,"3") = xx(i,"3"); *Calculate Probabilities *First, calculate x'beta ***s(i,n) = sum(k,x(i,k)*beta0(n,k)); *Denominator is: ***denom(i) = 1 + exp(s(i,"2")) + exp(s(i,"3")); *probability of the categories is: ***p1(i) = 1/denom(i); ***p2(i) = exp(s(i,"2"))/denom(i); ***p3(i) = exp(s(i,"3"))/denom(i); ***limit(i) = p1(i) + p2(i); ***u(i) = uniform(0,1); **Creating the LHS category variable b0 ***b0(i)$(u(i) le p1(i)) = 0; ***b0(i)$((u(i) gt p1(i)) and (u(i) le limit(i))) = 1; ***b0(i)$(u(i) gt limit(i)) = 2; parameter b0(i) / 1 0. 2 0. 3 0. 4 0. 5 0. 6 0. 7 0. 8 0. 9 0. 10 0. 11 0. 12 0. 13 0. 14 0. 15 0. 16 0. 17 0. 18 0. 19 0. 20 0. 21 0. 22 0. 23 0. 24 0. 25 0. 26 0. 27 0. 28 0. 29 0. 30 0. 31 0. 32 1.000000 33 1.000000 34 1.000000 35 1.000000 36 1.000000 37 1.000000 38 1.000000 39 1.000000 40 1.000000 41 1.000000 42 1.000000 43 1.000000 44 1.000000 45 1.000000 46 1.000000 47 1.000000 48 1.000000 49 1.000000 50 1.000000 51 1.000000 52 1.000000 53 1.000000 54 1.000000 55 1.000000 56 1.000000 57 1.000000 58 1.000000 59 1.000000 60 1.000000 61 1.000000 62 1.000000 63 1.000000 64 1.000000 65 1.000000 66 1.000000 67 1.000000 68 1.000000 69 1.000000 70 1.000000 71 1.000000 72 1.000000 73 1.000000 74 1.000000 75 1.000000 76 1.000000 77 1.000000 78 1.000000 79 1.000000 80 1.000000 81 1.000000 82 1.000000 83 1.000000 84 1.000000 85 1.000000 86 1.000000 87 1.000000 88 1.000000 89 1.000000 90 1.000000 91 1.000000 92 1.000000 93 1.000000 94 1.000000 95 1.000000 96 1.000000 97 1.000000 98 1.000000 99 1.000000 100 1.000000 101 2.000000 102 2.000000 103 2.000000 104 2.000000 105 2.000000 106 2.000000 107 2.000000 108 2.000000 109 2.000000 110 2.000000 111 2.000000 112 2.000000 113 2.000000 114 2.000000 115 2.000000 116 2.000000 117 2.000000 118 2.000000 119 2.000000 120 2.000000 121 2.000000 122 2.000000 123 2.000000 124 2.000000 125 2.000000 126 2.000000 127 2.000000 128 2.000000 129 2.000000 130 2.000000 131 2.000000 132 2.000000 133 2.000000 134 2.000000 135 2.000000 136 2.000000 137 2.000000 138 2.000000 139 2.000000 140 2.000000 141 2.000000 142 2.000000 143 2.000000 144 2.000000 145 2.000000 146 2.000000 147 2.000000 148 2.000000 149 2.000000 150 2.000000 151 2.000000 152 2.000000 153 2.000000 154 2.000000 155 2.000000 156 2.000000 157 2.000000 158 2.000000 159 2.000000 160 2.000000 161 2.000000 162 2.000000 163 2.000000 164 2.000000 165 2.000000 166 2.000000 167 2.000000 168 2.000000 169 2.000000 170 2.000000 171 2.000000 172 2.000000 173 2.000000 174 2.000000 175 2.000000 176 2.000000 177 2.000000 178 2.000000 179 2.000000 180 2.000000 181 2.000000 182 2.000000 183 2.000000 184 2.000000 185 3.000000 186 3.000000 187 3.000000 188 3.000000 189 3.000000 190 3.000000 191 3.000000 192 3.000000 193 3.000000 194 3.000000 195 3.000000 196 3.000000 197 3.000000 198 3.000000 199 3.000000 200 3.000000 201 3.000000 202 3.000000 203 3.000000 204 3.000000 205 3.000000 206 3.000000 207 3.000000 208 3.000000 209 3.000000 210 3.000000 211 3.000000 212 3.000000 213 3.000000 214 3.000000 215 3.000000 216 3.000000 217 3.000000 218 3.000000 219 3.000000 220 3.000000 221 3.000000 222 3.000000 223 3.000000 224 3.000000 225 3.000000 226 4.000000 227 4.000000 228 4.000000 229 4.000000 230 4.000000 231 4.000000 232 4.000000 233 4.000000 234 4.000000 235 4.000000 236 4.000000 237 4.000000 238 4.000000 239 4.000000 240 4.000000 241 4.000000 242 4.000000 243 4.000000 244 4.000000 245 4.000000 246 4.000000 247 4.000000 248 4.000000 249 4.000000 250 4.000000 251 4.000000 252 4.000000 253 4.000000 254 4.000000 255 4.000000 256 4.000000 257 4.000000 258 4.000000 259 4.000000 260 4.000000 261 4.000000 262 4.000000 263 4.000000 264 4.000000 265 4.000000 266 4.000000 267 4.000000 268 4.000000 269 4.000000 270 4.000000 271 4.000000 272 4.000000 273 4.000000 274 4.000000 275 4.000000 276 4.000000 277 4.000000 278 4.000000 279 4.000000 280 4.000000 281 4.000000 282 4.000000 283 4.000000 284 4.000000 285 4.000000 286 4.000000 287 4.000000 288 4.000000 289 4.000000 290 4.000000 291 4.000000 292 4.000000 293 4.000000 294 4.000000 295 4.000000 296 4.000000 297 4.000000 298 4.000000 299 4.000000 300 4.000000 301 4.000000 302 4.000000 303 4.000000 304 4.000000 305 4.000000 306 4.000000 307 4.000000 308 4.000000 309 4.000000 310 4.000000 311 4.000000 312 4.000000 313 4.000000 314 4.000000 315 4.000000 316 4.000000 317 4.000000 318 4.000000 319 4.000000 320 4.000000 321 4.000000 322 4.000000 323 4.000000 324 4.000000 325 4.000000 326 4.000000 327 4.000000 328 4.000000 329 4.000000 330 4.000000 331 4.000000 332 4.000000 333 4.000000 334 4.000000 335 4.000000 336 4.000000 337 4.000000 /; **Starting MC loops Loop(g, **End of data generation *display s; *display b0; *Converting the categiries b(i) into y(i,n) with 0 1 values Loop(i, b(i) = b0(i); ); b(i)$(b0(i)=0) = 8; b(i)$(b0(i)=1) = 9; y(i,n) = 0; Loop(i, y(i,n)=b(i); ); Loop(n, y(i,"1")$(y(i,n) EQ 8)=1; y(i,"2")$(y(i,n) EQ 9)=1; y(i,"3")$(y(i,n) EQ 2)=1; y(i,"4")$(y(i,n) EQ 3)=1; y(i,"5")$(y(i,n) EQ 4)=1; y(i,n)$(y(i,n) NE 1)=0; ); * display b0; * display y; **just stores data *xdata(i,"1")=b0(i); *xdata(i,"2")=x(i,"2"); *xdata(i,"3")=x(i,"3"); *xdata(i,"4")=x(i,"4"); options decimals=7; *display xdata; *initialization p.l(i,n)=1/NN; pe.l(i,n,h)=1/HH; tablepred(n,n) =0; lm.l(n,k)=0; *ME case solve base maximizing obj using nlp; *Dual GME Problem *solve dual minimizing obj using nlp; *CE case *solve base minimizing obj using nlp; *solve ols maximizing obj using nlp; options decimals=7; display p.l; ***Primal beta(g,n,k) = consist.m(n,k); **Dual *beta(g,n,k) = lm.l(n,k); *** *beta(g,n,k) = (1/(1-delta))*consist.m(n,k); *When normalizing Beta0=0 use the following beta(g,"1",k) = 0; options decimals=7; *display beta; options decimals=7; psi(i,n)=sum(h,exp(sum(k,x(i,k)*alphae(h)*lm.l(n,k)))); w0(i,n,h)=exp(sum(k,x(i,k)*alphae(h)*lm.l(n,k)))/psi(i,n); ei0(i,n)=sum(h,alphae(h)*w0(i,n,h)); display w0; display ei0; ***Calculating INF Matrix ME and ML IME1(g,n) = II + sum(i,exp(sum(k,beta(g,n,k)*x(i,k)))); omega1(g,i) = sum(n,exp(sum(k,beta(g,n,k)*x(i,k)))); IME3(g,n) = sum(i,(exp(sum(k,beta(g,n,k)*x(i,k)))) / (omega1(g,i)*omega1(g,i))); *IME3(g,n) = sum(i,1/p.l(i,n)); IME2(g,n) = sum(i,(1+exp(-sum(k,beta(g,n,k)*x(i,k))))/(exp(-sum(k,beta(g,n,k)*x( i,k))))); IEE(g,n) = sum(i,sum(h,1/pe.l(i,n,h))); **Note 1 and 2 are similar xsq(i) = sum(k,x(i,k)*x(i,k)); xsqk(k,k1) = sum(i,x(i,k)*x(i,k1)); *INFO(g,n,k,k1) = sum(i,(exp(sum(k3,beta(g,n,k3)*x(i,k3)))) / * (omega1(g,i)*omega1(g,i)))*xsqk(k1,k); INFO(g,n,k,k1) = sum(i,((exp(sum(k3,beta(g,n,k3)*x(i,k3))) / (omega1(g,i)*omega1(g,i))*(x(i,k)*x(i,k1))))); INFO1(g,n,k,k1) = sum(i,sum(h,((exp(sum(k3,beta(g,n,k3)*x(i,k3)*alphae(h))) / (psi(i,n)*psi(i,n))*((x(i,k)*alphae(h))*(x(i,k1)*alphae(h))))))); TINFO(g,n,k,k1) = INFO(g,n,k,k1) + INFO1(g,n,k,k1); IML1(g,n) = sum(i,((exp(sum(k,beta(g,n,k)*x(i,k)))) / (omega1(g,i)*omega1(g,i)))*xsq(i)); **original ML1 below *IML1(g,n) = sum(i, ((exp(sum(k,beta(g,n,k)*x(i,k))))/ * ((1+exp(sum(k,beta(g,n,k)* * x(i,k))))* (1+exp(sum(k,beta(g,n,k)*x(i,k)))))) * xsq(i)); IML2(g,n) = sum(i, ((exp(-sum(k,beta(g,n,k)*x(i,k))))/ ((1+exp(-sum(k,beta(g,n,k )*x(i,k))))* (1+exp(-sum(k,beta(g,n,k)*x(i,k)))))) ); ***End **PREDICTING with the recovered Betas xbeta(n,i) = sum(k,x(i,k)*beta(g,n,k)); ***Predicting according to Xbeta pred(i) = max(xbeta("1",i),xbeta("2",i),xbeta("3",i),xbeta("4",i),xbeta("5",i)); ypred(i)$(pred(i) eq xbeta("1",i)) = 0; ypred(i)$(pred(i) eq xbeta("2",i)) = 1; ypred(i)$(pred(i) eq xbeta("3",i)) = 2; ypred(i)$(pred(i) eq xbeta("4",i)) = 3; ypred(i)$(pred(i) eq xbeta("5",i)) = 4; ypred1(i)=ypred(i)-b0(i); **options decimals=7; tablep(i,"1","1")$(ypred(i) eq 0 and b0(i) eq 0) = 1; tablep(i,"1","2")$(ypred(i) eq 0 and b0(i) eq 1) = 1; tablep(i,"1","3")$(ypred(i) eq 0 and b0(i) eq 2) = 1; tablep(i,"1","4")$(ypred(i) eq 0 and b0(i) eq 3) = 1; tablep(i,"1","5")$(ypred(i) eq 0 and b0(i) eq 4) = 1; tablep(i,"2","1")$(ypred(i) eq 1 and b0(i) eq 0) = 1; tablep(i,"2","2")$(ypred(i) eq 1 and b0(i) eq 1) = 1; tablep(i,"2","3")$(ypred(i) eq 1 and b0(i) eq 2) = 1; tablep(i,"2","4")$(ypred(i) eq 1 and b0(i) eq 3) = 1; tablep(i,"2","5")$(ypred(i) eq 1 and b0(i) eq 4) = 1; tablep(i,"3","1")$(ypred(i) eq 2 and b0(i) eq 0) = 1; tablep(i,"3","2")$(ypred(i) eq 2 and b0(i) eq 1) = 1; tablep(i,"3","3")$(ypred(i) eq 2 and b0(i) eq 2) = 1; tablep(i,"3","4")$(ypred(i) eq 2 and b0(i) eq 3) = 1; tablep(i,"3","5")$(ypred(i) eq 2 and b0(i) eq 4) = 1; tablep(i,"4","1")$(ypred(i) eq 3 and b0(i) eq 0) = 1; tablep(i,"4","2")$(ypred(i) eq 3 and b0(i) eq 1) = 1; tablep(i,"4","3")$(ypred(i) eq 3 and b0(i) eq 2) = 1; tablep(i,"4","4")$(ypred(i) eq 3 and b0(i) eq 3) = 1; tablep(i,"4","5")$(ypred(i) eq 3 and b0(i) eq 4) = 1; tablep(i,"5","1")$(ypred(i) eq 4 and b0(i) eq 0) = 1; tablep(i,"5","2")$(ypred(i) eq 4 and b0(i) eq 1) = 1; tablep(i,"5","3")$(ypred(i) eq 4 and b0(i) eq 2) = 1; tablep(i,"5","4")$(ypred(i) eq 4 and b0(i) eq 3) = 1; tablep(i,"5","5")$(ypred(i) eq 4 and b0(i) eq 4) = 1; options decimals=7; display tablep; *loop(i1, *tablepred(g,n,i1) = sum(i,tablep(i,n,i1)); *); options decimals=7; display tablepred; **ypred1(i)$(ypred1(i) ne 0) = 1; **sumyp(g)=sum(i,ypred1(i)); **display sumyp; **Predicting directly with the recovered p's **pred(i) = max(p.l(i,"1"),p.l(i,"2"),p.l(i,"3")); **ypred(i)$(pred(i) eq p.l(i,"1")) = 0; **ypred(i)$(pred(i) eq p.l(i,"2")) = 1; **ypred(i)$(pred(i) eq p.l(i,"3")) = 2; *display ypred; *display b0; **tablep(i,"1","1")$(ypred(i) eq 0 and b0(i) eq 0) = 1; **tablep(i,"1","2")$(ypred(i) eq 0 and b0(i) eq 1) = 1; **tablep(i,"1","3")$(ypred(i) eq 0 and b0(i) eq 2) = 1; **tablep(i,"2","1")$(ypred(i) eq 1 and b0(i) eq 0) = 1; **tablep(i,"2","2")$(ypred(i) eq 1 and b0(i) eq 1) = 1; **tablep(i,"2","3")$(ypred(i) eq 1 and b0(i) eq 2) = 1; **tablep(i,"3","1")$(ypred(i) eq 2 and b0(i) eq 0) = 1; **tablep(i,"3","2")$(ypred(i) eq 2 and b0(i) eq 1) = 1; **tablep(i,"3","3")$(ypred(i) eq 2 and b0(i) eq 2) = 1; options decimals=7; *display tablep; loop(i1, tablepred(n,i1) = sum(i,tablep(i,n,i1)); ); options decimals=7; display tablepred; ypred1(i)=ypred(i)-b0(i); options decimals=7; *display ypred1; ypred1(i)$(ypred1(i) ne 0) = 1; sumyp(g)=sum(i,ypred1(i)); *display sumyp; dist(g,n,k) = beta(g,n,k) - (1*beta0(n,k)); MSE(g) = sum(n,sum(k,dist(g,n,k)*dist(g,n,k))) / ((NN-1)*KK); *display dist; *display MSE; ***Calculating P's and w's for Norm entropies omeg(i) = sum(n,exp(sum(k,x(i,k)*lm.l(n,k)))); recp(i,n) = exp(sum(k,x(i,k)*lm.l(n,k))) / omeg(i); psi(i,n) = sum(h,exp(sum(k,x(i,k)*lm.l(n,k)*alphae(h)))); recw(i,n,h) = exp(sum(k,x(i,k)*lm.l(n,k)*alphae(h))) / psi(i,n); entp(g) = -sum(i,sum(n,recp(i,n)*log(1.e-4+recp(i,n))))/(log(NN)*II); entpi(g,i) = -sum(n,p.l(i,n)*log(1.e-4+p.l(i,n)))/log(NN); options decimals=7; *display recp; *display recw; *display entpi; entpe(g) = -sum(i,sum(n,sum(h,recw(i,n,h)*log(1.e-4+recw(i,n,h)))))/(log(HH)*II*NN); options decimals=7; *display entpe; *entropy = sum(k,sum(j,pb.l(k,j)*log((1.e-7+pb.l(k,j))/1.e-7+qb(k,j)))); ); **End Monte Carlo Loop *Loop(g, *maxmiss$(sumyp(g) gt sumyp(g-1))=sumyp(g); *minmiss$(sumyp(g) lt sumyp(g-1))=sumyp(g); *); *display sumyp; *display maxmiss; *display minmiss; *minmiss = min(sumyp(g)); *maxmiss = max(sumyp(g)); options decimals=7; display sumyp; *display entp; *display entpi; *display entpe; display MSE; avbeta(n,k) = sum(g,beta(g,n,k))/GG; avIME1(n) = sum(g,IME1(g,n))/GG; avIME2(n) = sum(g,IME2(g,n))/GG; avIME3(n) = sum(g,IME3(g,n))/GG; avIML1(n) = sum(g,IML1(g,n))/GG; avIML2(n) = sum(g,IML2(g,n))/GG; varbeta(n,k) = sum(g,(beta(g,n,k) - avbeta(n,k))*(beta(g,n,k) - avbeta(n,k)))/GG; totvar = sum(n,sum(k,varbeta(n,k))); *loop(i1, *tablepred(g,n,i1) = sum(i,tablep(i,n,i1)); *); *avtablep(n,n) = sum(g,tablepred(g,n,n))/GG; *vartablep(n,n) = sum(g,(tablepred(g,n,n) - avtablep(n,n))*(tablepred(g,n,n) - avtablep(n,n))); *minmiss = min(sumyp(g)); *maxmiss = max(sumyp(g)); avmiss = sum(g,sumyp(g))/GG; varmiss = sum(g,(sumyp(g)-avmiss)*(sumyp(g)-avmiss)) / GG; avmse = sum(g,mse(g))/GG; varmse = sum(g,(mse(g)-avmse)*(mse(g)-avmse)) / GG; aventp = sum(g,entp(g))/GG; varentp = sum(g,(entp(g)-aventp)*(entp(g)-aventp)) / GG; aventpi(i) = sum(g,entpi(g,i))/GG; varentpi(i) = sum(g,(entpi(g,i)-aventpi(i))*(entpi(g,i)-aventpi(i))) / GG; aventpe = sum(g,entpe(g))/GG; varentpe = sum(g,(entpe(g)-aventpe)*(entpe(g)-aventpe)) / GG; options decimals=7; display avbeta; display varbeta; display totvar; display avmiss; display IME1; display IME2; display IME3; display IEE; display IML1; display IML2; display xsq; display xsqk; display INFO; display INFO1; display TINFO; display avIME1; display avIME2; display avIME3; display avIML1; display avIML2; *display minmiss; *display maxmiss; display varmiss; display alphae; display avmse; display varmse; display aventp; display varentp; display aventpi; display varentpi; display aventpe; display varentpe; options decimals=7; *display pe.l; options decimals=7; **display e1.l; *Note: below is the entropy for each beta *ent(n,k) = -sum(j,pb.l(n,k,j)*log(1.e-7+pb.l(n,k,j)))/log(JJ); *options decimals=7; *display ent; *entropy = sum(k,sum(j,pb.l(k,j)*log((1.e-7+pb.l(k,j))/1.e-7+qb(k,j)))); *entropy = -sum(n,sum(k,sum(j,pb.l(n,k,j)*log(1.e-7+pb.l(n,k,j))))); *options decimals=7; *display entropy;