MATH 483/583 Linear Statistical Modeling Morrow, Summer 2006

with(linalg):

Warning, the protected names norm and trace have been redefined and unprotected

X:=matrix(25,3,[1,7,560, 1,3,220, 1,3,340, 1,4,80, 1,6,150, 1,7,330, 1,2,110, 1,7,210, 1,30,1460, 1,5,605, 1,16,688, 1,10,215, 1,4,255, 1,6,462, 1,9,448, 1,10,776, 1,6,200, 1,7,132, 1,3,36, 1,17,770, 1,10,140, 1,26,810, 1,9,450, 1,8,635, 1,4,150]):

y:=matrix(25,1,[16.68,11.50,12.03,14.88,13.75,18.11,8.00,17.83,79.24,21.50,40.33,21.00, 13.50,19.75,24.00,29.00,15.35,19.00,9.50,35.10,17.90,52.32,18.75,19.83,10.75]):

XpX:=multiply(transpose(X),X); ypy:=multiply(transpose(y),y); Xpy:=multiply(transpose(X),y);

NiM+SSRYcFhHNiItSSdtYXRyaXhHNiRJKnByb3RlY3RlZEdGKUkoX3N5c2xpYkdGJTYjNyU3JSIjRCIkPiMiJkstIjclRi8iJWJJIicqKlE4NyVGMEYzIigpb0Ru

NiM+SSR5cHlHNiItSSdtYXRyaXhHNiRJKnByb3RlY3RlZEdGKUkoX3N5c2xpYkdGJTYjNyM3IyQiKiFIMUo9ISIl

NiM+SSRYcHlHNiItSSdtYXRyaXhHNiRJKnByb3RlY3RlZEdGKUkoX3N5c2xpYkdGJTYjNyU3IyQiJmdmJiEiIzcjJCInV3Z0RjA3IyQiKXBycUxGMA==

C:=evalf(inverse(XpX)); beta:=multiply(C,Xpy); betapXpy:=multiply(transpose(beta),Xpy); SSRes:=matadd(ypy,betapXpy,1,-1); one:=col(X,1): sumy:=multiply(transpose(y),one)[1];sigmahat2:=SSRes[1,1]/23;

NiM+SSJDRzYiLUknbWF0cml4RzYkSSpwcm90ZWN0ZWRHRilJKF9zeXNsaWJHRiU2IzclNyUkIitoPTpLNiEjNSQhK05LZltXISM3JCErIilwRG4kKSEjOTclRjEkIisiSCR5VkZGMyQhK2QnM2R5JUY2NyVGNEY6JCIrVVp1RzchIzo=

NiM+SSViZXRhRzYiLUknbWF0cml4RzYkSSpwcm90ZWN0ZWRHRilJKF9zeXNsaWJHRiU2IzclNyMkIio5SjdNIyEiKTcjJCIqQTJmaCJGMDcjJCIqaCNbUTkhIzU=

NiM+SSliZXRhcFhweUc2Ii1JJ21hdHJpeEc2JEkqcHJvdGVjdGVkR0YpSShfc3lzbGliR0YlNiM3IzcjJCIrTSgqbzI9ISIm

NiM+SSZTU1Jlc0c2Ii1JJ21hdHJpeEc2JEkqcHJvdGVjdGVkR0YpSShfc3lzbGliR0YlNiM3IzcjJCIpbUpQQiEiJg==

NiM+SSVzdW15RzYiJCImZ2YmISIj

NiM+SSpzaWdtYWhhdDJHNiIkIis0WUE7NSEiKQ==

SSR:= betapXpy[1,1]-sumy^2/25;

NiM+SSRTU1JHNiIkIiolNCIzYiYhIiY=

Example 3.8

sebeta1hat:=sqrt(sigmahat2*C[2,2]);

NiM+SStzZWJldGExaGF0RzYiJCIrIVE/KXA7ISM1

Example 3.9

x0:=vector([1,8,275]); h00:=multiply(multiply(transpose(x0),C),x0);

NiM+SSN4MEc2Ii1JJ3ZlY3Rvckc2JEkqcHJvdGVjdGVkR0YpSShfc3lzbGliR0YlNiM3JSIiIiIiKSIkdiM=

NiM+SSRoMDBHNiIkIit5XkQoUiYhIzY=

y0hat:=multiply(transpose(beta),x0);

NiM+SSZ5MGhhdEc2Ii1JJ3ZlY3Rvckc2JEkqcHJvdGVjdGVkR0YpSShfc3lzbGliR0YlNiM3IyQiKzM7VkE+ISIp

vary0hat:=sigmahat2*h00;

NiM+SSl2YXJ5MGhhdEc2IiQiK0xOI1tbJiEjNQ==

Example 3.13

xd:=vector([1,8,1200]); hdd:=multiply(multiply(transpose(xd),C),xd);

NiM+SSN4ZEc2Ii1JJ3ZlY3Rvckc2JEkqcHJvdGVjdGVkR0YpSShfc3lzbGliR0YlNiM3JSIiIiIiKSIlKzc=

NiM+SSRoZGRHNiIkIitVLWl0JykhIzU=

H:=multiply(multiply(X,C),transpose(X)): diagH:=vector(25):

for i from 1 to 25 do diagH[i]:=H[i,i] od:

print(diagH);

NiMtSSd2ZWN0b3JHNiRJKnByb3RlY3RlZEdGJkkoX3N5c2xpYkc2IjYjNzskIitpeCw9NSEjNSQiKyJvaiwyKCEjNiQiK1drWnQpKkYwJCIrayN6dWApRjAkIitzJlw1XShGMCQiK1ZHcCdHJUYwJCIrQXUnKXoiKUYwJCIrTSJmRFAnRjAkIiskZUBIKVxGLSQiK3IlZkgnPkYtJCIrSS1FOCcpRjAkIishKnBiTzZGLSQiK1tLWTdoRjAkIisxO0xDeUYwJCIrXHEyNlRGMCQiK0NWU2Y7Ri0kIisiKT4/VmZGMCQiKzJjL0UnKkYwJCIrI1FkW2sqRjAkIis0J1tvLCJGLSQiK2Mpb0ZsIkYtJCIrJT1fZCJSRi0kIitnXStFVEYwJCIrL0UzMTdGLSQiK3RhTWttRjA=

C:=evalf(inverse(XpX)): R:=multiply(C,transpose(X)):

for i to 25 do print(R[1,i]); od;

NiMkIitQJVI9XyQhIzY=

NiMkIisvVDlZIikhIzY=

NiMkIittSzJVciEjNg==

NiMkIitlMnFzKSkhIzY=

NiMkIitBVEYoUighIzY=

NiMkIitVJjNqVyYhIzY=

NiMkIisjcCw5XiohIzY=

NiMkIishUXouWCchIzY=

NiMkIStHYy9DOSEjNQ==

NiMkIis9Oi5OUyEjNg==

NiMkIStwTCFIYiIhIzY=

NiMkIitDXidSMiYhIzY=

NiMkIisneUklM3UhIzY=

NiMkIitXKipvJ3klISM2

NiMkIityTkRwTiEjNg==

NiMkIiplUiQqeiQhIzY=

NiMkIit0NyIqeXAhIzY=

NiMkIitEYS0uciEjNg==

NiMkIisqUT5kbyohIzY=

NiMkIStseChRbyMhIzY=

NiMkIispUjQ6cSYhIzY=

NiMkIStiPklBcSEjNg==

NiMkIitkIT5EYiQhIzY=

NiMkIitSZVZcQyEjNg==

NiMkIitwRipwRykhIzY=

sum(R[1,x]^2, x=1..25);

NiMkIitoPTpLNiEjNQ==

sqrt(%);

NiMkIitUa3VrTCEjNQ==

R[1,1]/%;

NiMkIis6em9ZNSEjNQ==