function mcmcExample(varargin)

%;----EXAMPLE CODE----

%;A simple pseudo-inversion

 optargin = size(varargin,2);
 switch optargin
     case 0
         noise=0.25;
     case 1
         noise=varargin{1};
     otherwise
         error('Too many input arguments');
 end

%;To run, start IDL
%;type .compile mcmc.pro
%;then type mcmc_example
%;type retall once done to break out of debug mode

%;let's make a simple PAR, temperature time series
  par1 =sin(pi*((1:8760)-12)/12)*2000;
  par1(par1<0)=0;

  par2 = sin(pi*(1:8760)/8760);
  par = par1.*par2;
  temp = sin(pi*(1:8760)/8760)*30;

%;Set up the X values
  x = [temp; par];
  
%;The "true" parameters and output
  truth = [1.5,0.1,10.0,100.0];
  nee = mcmc_testmodel(x,truth);


  	
  obs = nee + (randn(1,length(nee)).*noise);
  



param.name={'r1';'r2';'p1';'p2'};
param.value=[2.0; 0.3; 10.0; 500.0];
param.max=[10.0; 1.0; 100.0; 10000.0];
param.min=[0.0;0.001;0.0;10.0];
param.knob=[1.0;1.0;1.0;1.0];
param.changeable=[1;1;1;1];


%;---OPTIONAL KEYWORDS---
% (1) validdata is % of datapoints in each interval that are valid (default is 100%)  This has the same dimension as data.
% (2) valid_frac is min % of datapoints to accept (default is 50%)

% (3) model is the name of the model function (string) - default is "mcmc_testmodel"  (be sure to ignore the .m extension)
% (4) likelihood is the name of the likelihood function - it should call the model and compute the likelihood - default is "mcmc_likelihood"

% (5) numatonce - how often to check for chain convergence (default is 10000)
% (6) random_start - start at init guess if 0, else randomly within prior (default is 0)
% (7) numchains - number of chains (default is 10)
% (8) numspinups - how much burn in on final iteration (default is 125000)
% (9) iter - max number of iterations both for chains and final (default is 375000)
% (10) ranacc - tuning parameter for % of "worse" likelihoods to accept (default is -1) 
% (11) 'fast', 'medium', and 'superfast' are keywords with different default settings for the above settings - all are faster than default (see below)
% (12) 'quiet' turns off all printing of messages


%%% Empty fields will be set to the defaults in mcmc.m
inputOpts.model = 'mcmc_testmodel';
inputOpts.validdata=[];
inputOpts.valid_frac=[];
inputOpts.likelihood = [];
inputOpts.numatonce=[];
inputOpts.random_start=[];
inputOpts.numchains=[];
inputOpts.numspinups=5;
inputOpts.iter=10;
inputOpts.ranacc=[];
inputOpts.mode='superfast';
%inputOpts.mode=[];
inputOpts.quietMode=[];

[outputll,outputvalue,outputy]=mcmc(x,obs,param,inputOpts);
%;Call MCMC
%;Plot outputs
  plot(nee,outputy,'o'); axis([-10 10 -10 10]);
   t=sprintf(['Parameter Name / True Parameter Values / Estimated Parameter Values:']);
   disp(t) 
  for nC=1:length(truth)
      t=sprintf([char(param.name(nC)) ': ' num2str(truth(nC)) '   ' num2str(outputvalue(nC,1)) ]);
      disp(t)
  end
  
 corrMatrix=corrcoef(nee',outputy');
  t=sprintf(['r2: ' num2str(corrMatrix(1,2).^2) ]);
  disp(t)
 


 
