Chi-squared contours (xfig example)

Chi-squared contours

require("isisscripts");

variable x = [0:2.2*PI:.01];

variable alpha,x0,y0,range,col;
variable a,b;
variable i;

define ellipse(x,i,alpha,x0,y0,a,b){
  return {x0+cos(alpha)*(a-i*1.5)*cos(x)-sin(alpha)*(b-i*1.5)*sin(x),
	  y0+sin(alpha)*(a-i*1.5)*cos(x)+cos(alpha)*(b-i*1.5)*sin(x)};
}

x0=0.;y0=0.;
alpha=PI/4.;
range=10.;
a=4.,b=10.;

col=["red","blue","green4"];

vmessage("sin(alpha)=%f, cos(alpha=%f)",sin(alpha),cos(alpha));


variable X=xfig_plot_new(15,15);
X.world([-range:range],[-range:range]);

variable x_par_0 = ellipse(x,0,alpha,x0,y0,a,b)[0];
variable y_par_0 = ellipse(x,0,alpha,x0,y0,a,b)[1];
variable x_par_1 = ellipse(x,1,alpha,x0,y0,a,b)[0];
variable y_par_1 = ellipse(x,1,alpha,x0,y0,a,b)[1];
variable x_par_2 = ellipse(x,2,alpha,x0,y0,a,b)[0];
variable y_par_2 = ellipse(x,2,alpha,x0,y0,a,b)[1];

X.yaxis(; major=[min(y_par_0),min(y_par_1),min(y_par_2),0,max(y_par_0),max(y_par_1),max(y_par_2)],
	ticlabels=[`$3\sigma$`,
		   `$-2\sigma$`,
		   `$-1\sigma$`,
		   `$\hat{b}=0$`,
		   `$+3\sigma$`,
		   `$+2\sigma$`,
		   `$+1\sigma$`],
	ticlabels2=0);
X.y2axis(;off);
X.x2axis(;off);
X.x1axis(; major=[min(x_par_0),min(x_par_1),min(x_par_2),0,max(x_par_0),max(x_par_1),max(x_par_2)],
	ticlabels=[`$-3\sigma$`,
		   `$-2\sigma$`,
		   `$-1\sigma$`,
		   `$\hat{a}=0$`,
		   `$+3\sigma$`,
		   `$+2\sigma$`,
		   `$+1\sigma$`],
	ticlabels2=0);

X.plot([-10,0],[0,0];line=1);
X.plot([0,0],[-10,0];line=1);

for (i=0;i<3;i++){
  variable x_par = ellipse(x,i,alpha,x0,y0,a,b)[0];
  variable y_par = ellipse(x,i,alpha,x0,y0,a,b)[1];

  X.plot(x_par,y_par;color=col[i],width=3);
  X.plot([min(x_par),min(x_par)],[-10,y_par[where(x_par==min(x_par))]];line=1,color=col[i]);
  X.plot([max(x_par),max(x_par)],[-10,y_par[where(x_par==max(x_par))]];line=1,color=col[i]);
  X.plot([-10,x_par[where(y_par==min(y_par))]],[min(y_par),min(y_par)];line=1,color=col[i]);
  X.plot([-10,x_par[where(y_par==max(y_par))]],[max(y_par),max(y_par)];line=1,color=col[i]);  
  X.plot([0,0],[0,0];sym="x");
}

X.render("chisqr_ellipse.eps");