Friday, September 01, 2006

 

econgeek has left a new comment on your post "9/01/2006 06:28:00 AM":

and why exactly are you posting this here?

The post below was meant as a joke. More exactly it was the punchline to the joke at the end of the post immediately below it (posted just before it). Unfortunately, that post has been sent to the archive so this post is not just a semi debugged gauss program but also a very feeble punchline without a joke.


A Gauss program

n = 1000;
tb=20;

tg=1|1;

sn = 1;
se = 1;


z = rndn(n,2);
x = z*tg+sn*(rndn(n,1));
y = z*tg*tb+se*(rndn(n,1));




b1 = 0;
b2 =200;

c=8;


n=rows(y);
onen2= ones(n+2,1);

z = z|zeros(2,2);
x = x|zeros(2,1);
y = y|zeros(2,1);







r2s=(y-x*b1).*(y-x*b1);
r2s = sortc(r2s,1);
v1 = (0.5*(r2s[ floor((n+5)/2)]+r2s[floor((n+6)/2)]));



w1=(onen2 - (y-x*b1).^2/(c^2*v1)).* (((y-x*b1).^2).<(c*c*v1));

yw1 = y.*w1;
xw1 = x.*w1;
zw1 = z.*w1;



if b2==b1;
goto ddone;
endif;

gam1 = xw1/zw1;

/* here I orthogonalise zw1 so that it is easy to add an new observation that doesn't change gamma*/

mat1 = gam1~(0|1);
zw1a= zw1*mat1;
rhoz = zw1a[.,2]/zw1a[.,1];

mat2 = (1|0)~(-rhoz|1);

zw1a=zw1a*mat2;


czx1 = (zw1a[.,1])'xw1;
czy1 = (zw1a[.,1])'yw1;

zw1a[n+1,1]= meanc(abs(zw1a[.,1]));

x[n+1,1]=(czy1-czx1*b2)/((b2-b1)*zw1a[n+1,1]);
y[n+1,1]=x[n+1,1]*b1;

xw1[n+1,1]=x[n+1,1];
yw1[n+1,1]=y[n+1,1];

/*back to original z*/

zw1=zw1a*inv(mat2)*inv(mat1);

xhw1 = zw1*(xw1/zw1);


z[n+1,.]=zw1[n+1,.];

r2s=(y-x*b2).*(y-x*b2);
r2s = sortc(r2s,1);
v2 = (0.5*(r2s[ floor((n+3)/2)]+r2s[floor((n+4)/2)]));

w2=(onen2 - (y-x*b2).^2/(c^2*v2)).* (((y-x*b2).^2).<(c*c*v2));


/*do or redo obs n+2*/


zw2 = z.*w2;
xw2 = x.*w2;
yw2 = y.*w2;

rhoz2 = zw2[.,2]/zw1[.,1];

mat3 = (1|0)~(-rhoz2|1);

zw2a = zw2*mat3;



zw2a1 = zw2a[.,1];
zw2a2 = zw2a[.,2];


if ((zw2a1'yw2/zw2a1'xw2)==b2)*((zw2a2'yw2/zw2a2'xw2)==b2)*(w2[n+2]==0)==1;
goto ddone;
endif;

if ((zw2a1'yw2/zw2a1'xw2)==(zw2a2'y/zw2a1'x))*((zw2a1'y/zw2a1'x)==b2)*(w2[n+2]>0)==1;
goto redw1;
endif;

if ((zw2a1'yw2/zw2a1'xw2)==(zw2a2'y/zw2a1'x))*((zw2a1'y/zw2a1'x==b2)==0);

xw2a[n+2]=(xw2a/zw2a)*(meanc(abs(zw2a)));
zw2a[n+2,1]=meanc(abs(zw2a[.,1]));
yw2a[n+2]=b2*xw2a[n+2];
endif;

lam = (b2*zw2a2'xw2-zw2a2'yw2)/(zw2a1'yw2-b2*zw2a1'xw2);

mat4 = (lam|1)~(0|1);

zw2b=zw2a*mat4;

rhoz3 = (zw2b[.,2]/zw2b[.,1]);

mat5 = (1|0)~(-rhoz3|1);


zw2c = zw2b*mat5;


zw2c[n+2,2]= meanc(zw2c[.,2]);

xw2[n+2]=-(zw2c[1:n+1,2])'xw2[1:n+1]/zw2c[n+2,2];
yw2[n+2] = xw2[n+2]*b2;


/*back to original z*/

zw2 = zw2c*inv(mat5)*inv(mat4)*inv(mat3);

z[n+2,.]=zw2[n+2,.];
x[n+2]=xw2[n+2];
y[n+2]=yw2[n+2];


xhw2 = zw2*(xw2/zw2);

"check same final";

(yw2/xhw2)~b2;



redw1:

r2s=(y-x*b1).*(y-x*b1);
r2s = sortc(r2s,1);
v1 = (0.5*(r2s[ floor((n+3)/2)]+r2s[floor((n+4)/2)]));

w1=(onen2 - (y-x*b1).^2/(c^2*v1)).* (((y-x*b1).^2).<(c*c*v1));

yw1 = y.*w1;
xw1 = x.*w1;
zw1 = z.*w1;

xhw1 = zw1*(xw1/zw1);


ddone:;

"final check";
"estimate with w(" b1 "), estimate with w(" b2 "), and b2";



(yw1/xhw1)~(yw2/xhw2)~b2;
"estimate with w(" b1 ") -b2 and estimate with w(" b2 ") - b2";
((yw2/xhw2)-b2)~((yw1/xhw1)-b2);

"z[n+1,.], x[n+1], and y[n+1]"
z[n+1,.] x[n+1] y[n+1];


"z[n+2,.], x[n+2], and y[n+2]"
z[n+2,.] x[n+2] y[n+2];

posted by Robert permalink and comments6:28 AM