
%
% BUGS example
%
close all
clear all

curdir=pwd;
n=512; level=4;

wf=[  -0.075765714789341,  -0.029635527645954, ...
  0.497618667632458,   0.803738751805216, ...
  0.297857795605542,  -0.099219543576935, ...
 -0.012603967262261,   0.032223100604071];
ip=5;
alphatau1 = 0.05;  betatau1=0.5;    %mean= 0.05/0.5=0.1; var=0.05/0.25=0.2
alphatau2=0.5;  betatau2 = 0.05;    %mean= 0.5/0.05=10; var=0.5/0.025=20

s=MakeSignal('Blocks',n);
subplot(3,1,1);plot(s) 
axis tight
sn=s + 0.50 * randn(size(s));
subplot(3,1,2);plot(sn) 
axis tight

%DD=FWT_PO(sn, level, wf);
DD = dwtr(sn, log2(n)-level, wf);
    finest = log2(n)-1;
    is = std( DD(dyad(finest)));
    
dataStruct = struct( ...
    'DD', DD , ...
    'inputp', ip,...
    'inputs', is,...
    'alphatau1', alphatau1,...
    'betatau1', betatau1,...
    'alphatau2', alphatau2,...
    'betatau2', betatau2,...
    'N', n);

init0 = struct( 'tau', 1, 'tau1', 0.05,'tau2', 20 );
% ===============bugs.txt========================
% model{
%   for (j in 1:N){
%     DD[j] ~ dnorm(theta[j], tau);
%     theta[j] <- p[j] * xi1[j] + (1-p[j]) * xi2[j];
%     xi1[j] ~ ddexp(0, tau1);
%     xi2[j] ~ ddexp(0, tau2);
%     p[j] ~ dbern(r[j]);
%     r[j] ~ dbeta(1, inputp)
%    }
%   tau  ~  dgamma(1, inputs);
%   tau1 ~ dgamma(alphatau1, betatau1);
%   tau2 ~ dgamma(alphatau2, betatau2);
% }
%==============================================
cd('C:\Brani\MyBugs\matbugs\')
[samples, stats, structArray] = matbugs(dataStruct, ...
		fullfile(pwd, 'waveshrink.txt'), ...
		'init', init0, ...
        'openBugs', 1, ...
		'nChains', 1, ...
		'view', 0, 'nburnin', 1000, 'nsamples', 5000, ...
        'monitorParams', {'theta'}, ...
		'Bugdir', 'C:/Program Files/BUGS');
 cd(curdir)
%==============================================
%ss=IWT_PO(stats.mean.theta, level, wf);
ss=idwtr(stats.mean.theta, log2(n)-level, wf);

subplot(3,1,3);plot(ss) 
 hold on
 plot(s,'r-')
axis tight

%stats.DIC.total



