1 The SAS System 12:35 Thursday, October 27, 2005 NOTE: Copyright (c) 2002-2003 by SAS Institute Inc., Cary, NC, USA. NOTE: SAS (r) 9.1 (TS1M3) Licensed to UNIVERSITY OF FLORIDA-T/R, Site 0007029013. NOTE: This session is executing on the Linux 2.6.12-gentoo-r10 platform. NOTE: SAS 9.1.3 Service Pack 3 You are running SAS 9. Some SAS 8 files will be automatically converted by the V9 engine; others are incompatible. Please see http://support.sas.com/rnd/migration/planning/platform/64bit.html PROC MIGRATE will preserve current SAS file attributes and is recommended for converting all your SAS libraries from any SAS 8 release to SAS 9. For details and examples, please see http://support.sas.com/rnd/migration/index.html This message is contained in the SAS news file, and is presented upon initialization. Edit the file "news" in the "misc/base" directory to display site-specific news and information in the program log. The command line option "-nonews" will prevent this display. NOTE: SAS initialization used: real time 4.74 seconds cpu time 0.00 seconds 1 options nodate nonumber ps=500 ls=76; 2 3 data one; 4 retain brandnum 0; 5 * infile 'h:\public_html\data\winwine.dat' lrecl=218; 6 infile 'winwine.dat' lrecl=218; 7 * infile 'a:\winwine.dat' lrecl=218; 8 input brand $ 1-25 price score01-score23; 9 brandnum=brandnum+1; 10 if brandnum<=5 then do; 11 color=1; 12 white=1; 13 type=1; 14 end; 15 else if brandnum <=10 then do; 16 color=2; 17 white=0; 18 type=1; 19 end; 20 else if brandnum <=15 then do; 21 color=1; 22 white=1; 23 type=2; 24 end; 25 else if brandnum <=20 then do; 26 color=2; 27 white=0; 28 type=2; 29 end; 30 cprice=price-9.464; 31 score=score01; output; 32 score=score02; output; 33 score=score03; output; 34 score=score04; output; 35 score=score05; output; 36 score=score06; output; 37 score=score07; output; 38 score=score08; output; 39 score=score09; output; 40 score=score10; output; 41 score=score11; output; 42 score=score12; output; 43 score=score13; output; 44 score=score14; output; 45 score=score15; output; 46 score=score16; output; 47 score=score17; output; 48 score=score18; output; 49 score=score19; output; 50 score=score20; output; 51 score=score21; output; 52 score=score22; output; 53 score=score23; output; 54 drop score01-score23; 55 56 run; NOTE: The infile 'winwine.dat' is: File Name=/home/winner/public_html/data/winwine.dat, Owner Name=winner,Group Name=winner, Access Permission=rwxrw-r--, File Size (bytes)=4400 NOTE: 20 records were read from the infile 'winwine.dat'. The minimum record length was 218. The maximum record length was 218. One or more lines were truncated. NOTE: The data set WORK.ONE has 460 observations and 8 variables. NOTE: DATA statement used (Total process time): real time 0.00 seconds cpu time 0.02 seconds 56 ! 57 58 proc sort; 59 by brandnum; 60 NOTE: There were 460 observations read from the data set WORK.ONE. NOTE: The data set WORK.ONE has 460 observations and 8 variables. NOTE: PROCEDURE SORT used (Total process time): real time 0.00 seconds cpu time 0.00 seconds 61 data two; 62 set one; 63 by brandnum; 64 retain rater; 65 if first.brandnum then rater=0; 66 rater=rater+1; 67 run; NOTE: There were 460 observations read from the data set WORK.ONE. NOTE: The data set WORK.TWO has 460 observations and 9 variables. NOTE: DATA statement used (Total process time): real time 0.00 seconds cpu time 0.00 seconds 67 ! 68 69 proc sort data=two; 70 by rater; 71 run; NOTE: There were 460 observations read from the data set WORK.TWO. NOTE: The data set WORK.TWO has 460 observations and 9 variables. NOTE: PROCEDURE SORT used (Total process time): real time 0.00 seconds cpu time 0.00 seconds 71 ! 72 73 data twoa; 74 set two; 75 by rater; 76 retain numrate numwhite numred; 77 if first.rater then do; numrate=0; numwhite=0; numred=0; end; 78 if score ne . then do; 79 numrate=numrate+1; 80 if white=1 then numwhite=numwhite+1; 81 if white=0 then numred=numred+1; 82 end; 83 if last.rater then output; 84 keep rater numrate numwhite numred; 85 run; NOTE: There were 460 observations read from the data set WORK.TWO. NOTE: The data set WORK.TWOA has 23 observations and 4 variables. NOTE: DATA statement used (Total process time): real time 0.00 seconds cpu time 0.00 seconds 85 ! 86 87 proc print; 88 NOTE: There were 23 observations read from the data set WORK.TWOA. NOTE: The PROCEDURE PRINT printed page 1. NOTE: PROCEDURE PRINT used (Total process time): real time 0.01 seconds cpu time 0.01 seconds 89 proc sort data=two; 90 by rater brandnum; 91 run; NOTE: There were 460 observations read from the data set WORK.TWO. NOTE: The data set WORK.TWO has 460 observations and 9 variables. NOTE: PROCEDURE SORT used (Total process time): real time 0.00 seconds cpu time 0.01 seconds 91 ! 92 93 data three; 94 set two; 95 if score ne .; 96 if brandnum=1 then brandnum01=1; else brandnum01=0; 97 if brandnum=2 then brandnum02=1; else brandnum02=0; 98 if brandnum=3 then brandnum03=1; else brandnum03=0; 99 if brandnum=4 then brandnum04=1; else brandnum04=0; 100 if brandnum=5 then brandnum05=1; else brandnum05=0; 101 if brandnum=6 then brandnum06=1; else brandnum06=0; 102 if brandnum=7 then brandnum07=1; else brandnum07=0; 103 if brandnum=8 then brandnum08=1; else brandnum08=0; 104 if brandnum=9 then brandnum09=1; else brandnum09=0; 105 if brandnum=10 then brandnum10=1; else brandnum10=0; 106 if brandnum=11 then brandnum11=1; else brandnum11=0; 107 if brandnum=12 then brandnum12=1; else brandnum12=0; 108 if brandnum=13 then brandnum13=1; else brandnum13=0; 109 if brandnum=14 then brandnum14=1; else brandnum14=0; 110 if brandnum=15 then brandnum15=1; else brandnum15=0; 111 if brandnum=16 then brandnum16=1; else brandnum16=0; 112 if brandnum=17 then brandnum17=1; else brandnum17=0; 113 if brandnum=18 then brandnum18=1; else brandnum18=0; 114 if brandnum=19 then brandnum19=1; else brandnum19=0; 115 if brandnum=20 then brandnum20=1; else brandnum20=0; 116 117 if rater=1 then do; 118 if white=1 then do; z01=1; z02=0; end; 119 else do; z01=0; z02=1; end; end; 120 121 else if rater=2 then do; 122 if white=1 then do; z03=1; z04=0; end; 123 else do; z03=0; z04=1; end; end; 124 125 else if rater=3 then do; 126 if white=1 then do; z05=1; z06=0; end; 127 else do; z05=0; z06=1; end; end; 128 129 else if rater=4 then do; 130 if white=1 then do; z07=1; z08=0; end; 131 else do; z07=0; z08=1; end; end; 132 133 else if rater=5 then do; 134 if white=1 then do; z09=1; z10=0; end; 135 else do; z09=0; z10=1; end; end; 136 137 else if rater=6 then do; 138 if white=1 then do; z11=1; z12=0; end; 139 else do; z11=0; z12=1; end; end; 140 141 else if rater=7 then do; 142 if white=1 then do; z13=1; z14=0; end; 143 else do; z13=0; z14=1; end; end; 144 145 else if rater=8 then do; 146 if white=1 then do; z15=1; z16=0; end; 147 else do; z15=0; z16=1; end; end; 148 149 else if rater=9 then do; 150 if white=1 then do; z17=1; z18=0; end; 151 else do; z17=0; z18=1; end; end; 152 153 else if rater=10 then do; 154 if white=1 then do; z19=1; z20=0; end; 155 else do; z19=0; z20=1; end; end; 156 157 else if rater=11 then do; 158 if white=1 then do; z21=1; z22=0; end; 159 else do; z21=0; z22=1; end; end; 160 161 else if rater=12 then do; 162 if white=1 then do; z23=1; z24=0; end; 163 else do; z23=0; z24=1; end; end; 164 165 else if rater=13 then do; 166 if white=1 then do; z25=1; z26=0; end; 167 else do; z25=0; z26=1; end; end; 168 169 else if rater=14 then do; 170 if white=1 then do; z27=1; z28=0; end; 171 else do; z27=0; z28=1; end; end; 172 173 else if rater=15 then do; 174 if white=1 then do; z29=1; z30=0; end; 175 else do; z29=0; z30=1; end; end; 176 177 else if rater=16 then do; 178 if white=1 then do; z31=1; z32=0; end; 179 else do; z31=0; z32=1; end; end; 180 181 else if rater=17 then do; 182 if white=1 then do; z33=1; z34=0; end; 183 else do; z33=0; z34=1; end; end; 184 185 else if rater=18 then do; 186 if white=1 then do; z35=1; z36=0; end; 187 else do; z35=0; z36=1; end; end; 188 189 else if rater=19 then do; 190 if white=1 then do; z37=1; z38=0; end; 191 else do; z37=0; z38=1; end; end; 192 193 else if rater=20 then do; 194 if white=1 then do; z39=1; z40=0; end; 195 else do; z39=0; z40=1; end; end; 196 197 else if rater=21 then do; 198 if white=1 then do; z41=1; z42=0; end; 199 else do; z41=0; z42=1; end; end; 200 201 else if rater=22 then do; 202 if white=1 then do; z43=1; z44=0; end; 203 else do; z43=0; z44=1; end; end; 204 205 else if rater=23 then do; 206 if white=1 then do; z45=1; z46=0; end; 207 else do; z45=0; z46=1; end; end; 208 209 210 keep white cprice score rater brandnum01-brandnum20 z01-z46; 211 run; NOTE: There were 460 observations read from the data set WORK.TWO. NOTE: The data set WORK.THREE has 400 observations and 70 variables. NOTE: DATA statement used (Total process time): real time 0.00 seconds cpu time 0.00 seconds 211 ! 212 213 proc sort data=three; 214 by rater descending white; 215 216 quit; NOTE: There were 400 observations read from the data set WORK.THREE. NOTE: The data set WORK.THREE has 400 observations and 70 variables. NOTE: PROCEDURE SORT used (Total process time): real time 0.00 seconds cpu time 0.00 seconds 216 ! 217 218 219 proc iml; NOTE: IML Ready 220 use three; 220 ! 221 read all into xyzmat; 221 ! 222 223 use twoa; 223 ! 224 read all into numrate; 224 ! 225 226 numwhite=numrate(|,3|); 226 ! 227 numred=numrate(|,4|); 227 ! 228 229 print numwhite numred; 229 ! 230 231 x0=j(nrow(xyzmat),1,1); 231 ! 232 233 y=xyzmat(|,3|); 233 ! 234 235 xbrand=xyzmat(|,5:24|); 235 ! 236 237 x=xbrand; 237 ! 238 239 240 z=xyzmat(|,25:70|); 240 ! 241 242 do i1=1 to 400; 242 ! 243 do i2=1 to 46; 243 ! 244 if z(|i1,i2|)=. then z(|i1,i2|)=0; 244 ! 245 end; 245 ! 246 end; 246 ! 247 248 bols=inv(x`*x)*x`*y; 248 ! 249 250 e=y-x*bols; 250 ! 251 white_e=xyzmat(|,1|)#e; 251 ! 252 red_e=(j(nrow(xyzmat),1,1)-xyzmat(|,1|))#e; 252 ! 253 254 n_white=xyzmat(|,1|)`*xyzmat(|,1|); 254 ! 255 n_red=nrow(xyzmat)-n_white; 255 ! 256 257 white_var=white_e`*white_e/n_white; 257 ! 258 red_var=red_e`*red_e/n_red; 258 ! 259 260 ybarw=j(23,1,0); 260 ! 261 ybarr=j(23,1,0); 261 ! 262 ybars=j(46,1,0); 262 ! 263 264 do i1=1 to 23; 264 ! 265 if numwhite(|i1,|)>0 then do; 265 ! 266 ybarw(|i1,|)=z(|,2*i1-1|)`*y/numwhite(|i1,|); 266 ! 267 ybars(|2*i1-1,|)=ybarw(|i1,|); 267 ! 268 end; 268 ! 269 if numred(|i1,|)>0 then do; 269 ! 270 ybarr(|i1,|)=z(|,2*i1|)`*y/numred(|i1,|); 270 ! 271 ybars(|2*i1,|)=ybarr(|i1,|); 271 ! 272 end; 272 ! 273 end; 273 ! 274 275 eran=y-(x*bols)-(z*ybars)+(j(nrow(xyzmat),1,1/nrow(xyzmat))`*y); 275 ! 276 277 ran_var=eran`*eran/nrow(xyzmat); 277 ! 278 279 print white_var red_var n_white n_red ran_var; 279 ! 280 281 * print e white_e red_e eran; 282 283 print ybarw numwhite ybarr numred ybars; 283 ! 284 285 286 287 288 * print y; 289 * print x; 290 * print z; 291 292 thold=j(3,1,0); 292 ! 293 294 thold(|1,|)=(white_var-ran_var)/2; 294 ! 295 thold(|2,|)=red_var-ran_var; 295 ! 296 thold(|3,|)=ran_var; 296 ! 297 298 * thold(|3,|)=1; 299 300 thdiff=1000; 300 ! 301 302 iter=0; 302 ! 303 start newraph; 303 ! 304 305 do until(thdiff < .0001); 305 ! 306 307 iter=iter+1; 307 ! 308 gmat=j(46,46,0); 308 ! 309 310 do i1=1 to 46; 310 ! 311 do i2=1 to 46; 311 ! 312 313 if i1=i2 then do; 313 ! 314 if mod(i1,2)=1 then do; 314 ! gmat(|i1,i2|)=thold(|1,1|); 314 ! end; 314 ! 315 if mod(i1,2)=0 then do; 315 ! gmat(|i1,i2|)=thold(|2,1|); 315 ! end; 315 ! 316 end; 316 ! The SAS System 317 318 319 end; 319 ! 320 end; 320 ! 321 322 * print gmat; 323 324 rmat=thold(|3,1|)*i(nrow(xyzmat)); 324 ! 325 326 vmat=z*gmat*z`+rmat; 326 ! 327 328 * print 'Start det(vmat)'; 329 330 * detvmat=det(vmat); 331 332 * print 'end det(vmat)'; 333 /* 334 print gmat; 335 print rmat; 336 print vmat; 337 print detvmat; 338 */ 339 340 lorowi=0; 340 ! 341 hirowi=0; 341 ! 342 ivmat=j(400,400,0); 342 ! 343 detvmat=1; 343 ! 344 345 do i1=1 to 23; 345 ! 346 * print i1 lorowi hirowi; 347 hirowi=hirowi+numrate(|i1,2|); 347 ! 348 vmati=j(numrate(|i1,2|),numrate(|i1,2|),0); 348 ! 349 vmati=vmat(|lorowi+1:hirowi,lorowi+1:hirowi|); 349 ! 350 ivmat(|lorowi+1:hirowi,lorowi+1:hirowi|)=inv(vmati); 350 ! 351 detvmat=detvmat*det(vmati); 351 ! 352 lorowi=hirowi; 352 ! 353 end; 353 ! 354 355 356 357 358 b=inv(x`*ivmat*x)*x`*ivmat*y; 358 ! 359 360 rvec=y-x*b; 360 ! 361 362 p1=log(detvmat); 362 ! 363 p2=rvec`*ivmat*rvec; 363 ! 364 p3a=x`*ivmat*x; 364 ! 365 p3b=det(p3a); 365 ! 366 p3=log(p3b); 366 ! 367 * print p3a p3b; 368 * p3=log(det(x`*ivmat*x)); 369 370 m2lml=p1+p2; 370 ! 371 m2lrml=p1+p2+p3; 371 ! 372 373 374 375 qmat=root(inv(x`*ivmat*x)); 375 ! 376 xs=x*qmat`; 376 ! 377 378 dv1=j(46,46,0); 378 ! 379 dv2=j(46,46,0); 379 ! 380 381 382 do i1=1 to 46; 382 ! 383 do i2=1 to 46; 383 ! 384 385 if i1=i2 then do; 385 ! 386 if mod(i1,2)=1 then do; 386 ! dv1(|i1,i2|)=1; 386 ! end; 386 ! 387 if mod(i1,2)=0 then do; 387 ! dv2(|i1,i2|)=1; 387 ! end; 387 ! 388 end; 388 ! 389 390 391 end; 391 ! 392 end; 392 ! 393 394 * print dv1; 395 * print dv2; 396 397 398 vdot1=z*dv1*z`; 398 ! 399 vdot2=z*dv2*z`; 399 ! 400 401 402 403 404 lg11=trace(ivmat*vdot1); 404 ! 405 lg12=trace(ivmat*vdot2); 405 ! 406 lg13=trace(ivmat); 406 ! 407 408 lg21=-rvec`*ivmat*vdot1*ivmat*rvec; 408 ! 409 lg22=-rvec`*ivmat*vdot2*ivmat*rvec; 409 ! 410 lg23=-rvec`*ivmat*ivmat*rvec; 410 ! 411 412 lg31=-trace(xs`*ivmat*vdot1*ivmat*xs); 412 ! 413 lg32=-trace(xs`*ivmat*vdot2*ivmat*xs); 413 ! 414 lg33=-trace(xs`*ivmat*ivmat*xs); 414 ! 415 416 lg1=lg11+lg21+lg31; 416 ! 417 lg2=lg12+lg22+lg32; 417 ! 418 lg3=lg13+lg23+lg33; 418 ! 419 lg=lg1//lg2//lg3; 419 ! 420 421 bg111=-trace(ivmat*vdot1*ivmat*vdot1); 421 ! 422 bg112=-trace(ivmat*vdot1*ivmat*vdot2); 422 ! 423 bg113=-trace(ivmat*vdot1*ivmat); 423 ! 424 425 bg121=-trace(ivmat*vdot2*ivmat*vdot1); 425 ! 426 bg122=-trace(ivmat*vdot2*ivmat*vdot2); 426 ! 427 bg123=-trace(ivmat*vdot2*ivmat); 427 ! 428 429 430 bg131=-trace(ivmat*ivmat*vdot1); 430 ! 431 bg132=-trace(ivmat*ivmat*vdot2); 431 ! 432 bg133=-trace(ivmat*ivmat); 432 ! 433 434 bg11=bg111||bg112||bg113; 434 ! 435 bg12=bg121||bg122||bg123; 435 ! 436 bg13=bg131||bg132||bg133; 436 ! 437 bg1=bg11//bg12//bg13; 437 ! 438 439 bg211=2*(rvec`*ivmat*vdot1*ivmat*vdot1*ivmat*rvec)- 440 2*(rvec`*ivmat*vdot1*ivmat*xs*xs`*ivmat*vdot1*ivmat*rvec); 440 ! 441 bg212=2*(rvec`*ivmat*vdot1*ivmat*vdot2*ivmat*rvec)- 442 2*(rvec`*ivmat*vdot1*ivmat*xs*xs`*ivmat*vdot2*ivmat*rvec); 442 ! 443 bg213=2*(rvec`*ivmat*vdot1*ivmat*ivmat*rvec)- 444 2*(rvec`*ivmat*vdot1*ivmat*xs*xs`*ivmat*ivmat*rvec); 444 ! 445 446 bg221=2*(rvec`*ivmat*vdot2*ivmat*vdot1*ivmat*rvec)- 447 2*(rvec`*ivmat*vdot2*ivmat*xs*xs`*ivmat*vdot1*ivmat*rvec); 447 ! 448 bg222=2*(rvec`*ivmat*vdot2*ivmat*vdot2*ivmat*rvec)- 449 2*(rvec`*ivmat*vdot2*ivmat*xs*xs`*ivmat*vdot2*ivmat*rvec); 449 ! 450 bg223=2*(rvec`*ivmat*vdot2*ivmat*ivmat*rvec)- 451 2*(rvec`*ivmat*vdot2*ivmat*xs*xs`*ivmat*ivmat*rvec); 451 ! 452 453 bg231=2*(rvec`*ivmat*ivmat*vdot1*ivmat*rvec)- 454 2*(rvec`*ivmat*ivmat*xs*xs`*ivmat*vdot1*ivmat*rvec); 454 ! 455 bg232=2*(rvec`*ivmat*ivmat*vdot2*ivmat*rvec)- 456 2*(rvec`*ivmat*ivmat*xs*xs`*ivmat*vdot2*ivmat*rvec); 456 ! 457 bg233=2*(rvec`*ivmat*ivmat*ivmat*rvec)-2*(rvec`*ivmat*ivmat*xs*xs 457 ! `*ivmat*ivmat*rvec); 458 459 bg21=bg211||bg212||bg213; 459 ! 460 bg22=bg221||bg222||bg223; 460 ! 461 bg23=bg231||bg232||bg233; 461 ! 462 bg2=bg21//bg22//bg23; 462 ! 463 464 bg311=(2*trace(xs`*ivmat*vdot1*ivmat*vdot1*ivmat*xs)) 465 -trace(xs`*ivmat*vdot1*ivmat*xs*xs`*ivmat*vdot1*ivmat*xs) 465 ! ; 466 bg312=(2*trace(xs`*ivmat*vdot1*ivmat*vdot2*ivmat*xs)) 467 -trace(xs`*ivmat*vdot1*ivmat*xs*xs`*ivmat*vdot2*ivmat*xs) 467 ! ; 468 bg313=(2*trace(xs`*ivmat*vdot1*ivmat*ivmat*xs)) 469 -trace(xs`*ivmat*vdot1*ivmat*xs*xs`*ivmat*ivmat*xs); 469 ! 470 471 bg321=(2*trace(xs`*ivmat*vdot2*ivmat*vdot1*ivmat*xs)) 472 -trace(xs`*ivmat*vdot2*ivmat*xs*xs`*ivmat*vdot1*ivmat*xs) 472 ! ; 473 bg322=(2*trace(xs`*ivmat*vdot2*ivmat*vdot2*ivmat*xs)) 474 -trace(xs`*ivmat*vdot2*ivmat*xs*xs`*ivmat*vdot2*ivmat*xs) 474 ! ; 475 bg323=(2*trace(xs`*ivmat*vdot2*ivmat*ivmat*xs)) 476 -trace(xs`*ivmat*vdot2*ivmat*xs*xs`*ivmat*ivmat*xs); 476 ! 477 478 bg331=(2*trace(xs`*ivmat*ivmat*vdot1*ivmat*xs)) 479 -trace(xs`*ivmat*ivmat*xs*xs`*ivmat*vdot1*ivmat*xs); 479 ! 480 bg332=(2*trace(xs`*ivmat*ivmat*vdot2*ivmat*xs)) 481 -trace(xs`*ivmat*ivmat*xs*xs`*ivmat*vdot2*ivmat*xs); 481 ! 482 bg333=(2*trace(xs`*ivmat*ivmat*ivmat*xs))-trace(xs`*ivmat*ivmat*x 482 ! s*xs`*ivmat*ivmat*xs); 483 484 bg31=bg311||bg312||bg313; 484 ! 485 bg32=bg321||bg322||bg323; 485 ! 486 bg33=bg331||bg332||bg333; 486 ! 487 bg3=bg31//bg32//bg33; 487 ! 488 489 bg=bg1+bg2+bg3; 489 ! 490 491 thnew=thold-inv(bg)*lg; 491 ! 492 493 thdiff=(thnew-thold)`*(thnew-thold); 493 ! 494 495 print thnew thold thdiff; 495 ! 496 497 498 499 thold=thnew; 499 ! 500 501 end; 501 ! 502 finish; NOTE: Module NEWRAPH defined. 502 ! 503 run newraph; 503 ! 504 505 invbg=inv(bg); 505 ! 506 dinvbg=diag(invbg); 506 ! 507 sethnew=sqrt(dinvbg*j(3,1,1)); 507 ! 508 509 print iter thnew sethnew b; 509 ! 510 511 512 513 quit; NOTE: Exiting IML. NOTE: 255 workspace compresses. NOTE: The PROCEDURE IML printed page 2. NOTE: PROCEDURE IML used (Total process time): real time 22.44 seconds cpu time 20.48 seconds 513 ! NOTE: SAS Institute Inc., SAS Campus Drive, Cary, NC USA 27513-2414 NOTE: The SAS System used: real time 27.23 seconds cpu time 20.53 seconds