answer=input("Generate random data ? (1/0)");
if answer==1 then

  // do only once
  lx=input("\nHow many data points ?");
  x=(-0.5*lx:0.5*lx)/lx;
  
  N=length(x);
  x2=ones(1:N);
  I=eye(2,2);

  // this is the input for the ideal line 1,0
  Di=[x',x2'];

  // generate input data by adding random noise to straight line

  sigma_noise=input("Sigma for data generation ?");
  x3=zeros(1:N);

  for k=1:N,
    x3(k)=grand(1,1,'nor',0,sigma_noise);
  end;
  t=x+x3;

  xinit('Data');
  plot2d(x,x,strf='055',rect=[-0.6,-0.6,0.6,0.6]);
  plot2d(x,t,style=-1,strf='055',rect=[-0.6,-0.6,0.6,0.6]);
  
end,

//
//  Define weight ranges
//

minw1=0.5;
maxw1=1.5;
minw2=-0.5;
maxw2=0.5;

lw=50;
w1=(minw1*lw:maxw1*lw)/lw;
w2=(minw2*lw:maxw2*lw)/lw;



//
// Define prior function and normalising constant
//
s_w=input("Sigma for prior ? (1/sqrt(alpha))");
alpha=1/(s_w^2);
Z_w=1/(2*%pi*s_w^2);
deff('[y0]=prior(w3,w4)','y0=Z_w*exp(-(w3^2+w4^2)/(2*s_w^2))');
xinit('Prior');
fplot3d(w1,w2,prior,-120,89.8,"W1@W2@P(w|H)",[2,1,4],[minw1,maxw1,minw2,maxw2,0,Z_w]);


// Get sigma_nu
sigma_nu=input("Sigma for noise model ? (1/sqrt(beta))");
beta=1/(sigma_nu^2);
Z_nu=1/(2*%pi*(sigma_nu^2))^(N/2);

//Define functions
deff('[y1]=lin(w3,w4)','y1=w3*x+w4');
deff('[y2]=delta(t,x)','y2=(t-x)^2');
deff('[y3]=ED(w3,w4)','y3=sum(delta(t,lin(w3,w4)))');
deff('[y4]=likelihood(w3,w4)','y4=Z_nu*exp(-ED(w3,w4)/(sigma_nu^2))');
deff('[y5]=posterior(w3,w4)','y5=likelihood(w3,w4)*prior(w3,w4)');

z1=feval(w1,w2,ED);
[g1,g2]=min(z1);
w1_mp=w1(g2(1));
w2_mp=w2(g2(2));
mprintf('\nValues maximising likelihood: a=%4.2f b=%4.2f \n',w1_mp,w2_mp);
A=likelihood(w1_mp,w2_mp);
deff('[y6]=norm_likelihood(w3,w4)','y6=likelihood(w3,w4)/A');

z2=feval(w1,w2,posterior);
[g1,g2]=max(z2);
w1_MP=w1(g2(1));
w2_MP=w2(g2(2));
mprintf('Values maximising posterior distribution: a=%4.2f b=%4.2f\n',w1_MP,w2_MP);
B=posterior(w1_MP,w2_MP);
deff('[y7]=norm_posterior(w3,w4)','y7=posterior(w3,w4)/B');

wait=input("\n\t\treturn to see P(D|w,h)\n","string");
xinit('likelihood');
fplot3d(w1,w2,norm_likelihood,-70,55,"W1@W2@P(D|w,H)",[2,1,4],[minw1,maxw1,minw2,maxw2,0,1]);

wait=input("\n\t\treturn to see P(w|D,H)\n","string");
xinit('posterior');
fplot3d(w1,w2,norm_posterior,-70,55,"W1@W2@P(w|D,H)",[2,1,4],[minw1,maxw1,minw2,maxw2,0,1]);


//Calculate nabla nabla M and invert
x_square=x.^2;
temp=[sum(x_square),sum(x);sum(x),sum(x2)];
inv_V=alpha*I+beta*temp;
V=inv(inv_V);

//Calculate best estimate sigma_nu

Chi_square_target=N-(2-alpha*trace(V));
sigma_nu_opt=1/sqrt(Chi_square_target/(2*ED(w1_mp,w2_mp)));
mprintf("\n\t\tOptimised value for sigma_nu: %5.4f\n\n",sigma_nu_opt);

wait=input("\n\t\tPress Return to see predictions\n\n","string");

// new x for predictions
x_save=x;
x=-2:0.1:2;
N_predict=length(x);
Di_save=Di;
x2_predict=ones(1:N_predict);
Di=[x',x2_predict'];
SN=(sigma_nu^2)*ones(1:N_predict);

Y=lin(w1_MP,w2_MP);
Y_UP_1=lin(w1_MP,w2_MP)+sqrt(diag(Di*V*Di')');
Y_LOW_1=lin(w1_MP,w2_MP)-sqrt(diag(Di*V*Di')');
Y_UP_2=lin(w1_MP,w2_MP)+sqrt(diag(Di*V*Di')'+SN);
Y_LOW_2=lin(w1_MP,w2_MP)-sqrt(diag(Di*V*Di')'+SN);

xinit('Predictions');
plot2d(x_save,t,style=-1,strf='055',rect=[min(x),min(Y),max(x),max(Y)]);
plot2d(x,Y,strf='055',rect=[min(x),min(Y),max(x),max(Y)]);
plot2d(x,Y_UP_1,style=2,strf='055',rect=[min(x),min(Y),max(x),max(Y)]);
plot2d(x,Y_LOW_1,style=2,strf='055',rect=[min(x),min(Y),max(x),max(Y)]);
plot2d(x,Y_UP_2,style=3,strf='055',rect=[min(x),min(Y),max(x),max(Y)]);
plot2d(x,Y_LOW_2,style=3,strf='055',rect=[min(x),min(Y),max(x),max(Y)]);

//

x=x_save;
Di=Di_save;

answer=input("\n\t\tReturn to quit\n\n");
xdel(winsid());