function [A,B] = Gibbs_Regress(N)

  %%%%% construct data
          A = [];
          B = []; 
          x = 1:100;
          x_1 = ones(1,100)';
          x_2 = sin(3.14159*x/4)';
          x_3 = cos(3.14159*x/8)';
          sd = 4;
          
           u = normrnd(0,sd, [1,100])';
           
           b1 = 1;
           b2 = 1;
           b3 = 3;
           
           y = b1*x_1+b2*x_2+b3*x_3+u;        
  %%%%% prior information
  
       beta_0 = [1 1 1]';
       cov_0 = eye(3,3);
       invcov_0 = inv(cov_0);
       nu_0  = 2;
       delta_0 = 2;     
  %%%%%  data manipulation
   X = [x_1,x_2,x_3];  %% covariate matrix
  
   nu_1 = nu_0+100;
  %%%%% Gibbs Sampler
  var = 10; %%%%%% starting sample of sigma^2
  for i = 1:N
      cov_1 = inv(X'*X/var+invcov_0);
      beta_hat = cov_1*(X'*y/var+invcov_0*beta_0);
      beta = mvnrnd(beta_hat',cov_1)';               %%% sample conditional of beta
       error  = (y-X*beta);
      delta_1 = delta_0+ error'*error;
      
      temp = gamrnd(nu_1/2, (delta_1/2)^(-1));     %%% sample conditional of sigma^2
       var = temp^(-1);
       A = [A;beta'];
       B = [B, var];    
  end