-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmain_SMC.m
153 lines (146 loc) · 6.71 KB
/
main_SMC.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
clear;
close all;
%% Generate simulated data
load data_simulated; % load dataset that has been generated or generate the data as follows
% P = 2; % number of covariates
% N = 1e4; % number of test statistcs
% betatrue=[-3.5 repmat(1/sqrt(P),1,P)]; % the regression coefficient \beta
% X=randn(N,P); % covariates
% psi=[ones(N,1) X]*betatrue';
% wsuccess=1./(1+exp(-psi)); % the prior probability that a test statistic is a signal (alternative hypothesis)
% % Some theta's are signals, most are noise
% gammatrue=binornd(ones(N,1),wsuccess,N,1); % generate the true hypotheses in the simulation. 1: alternative; 0: null
% % Density of signals
% thetatrue = normrnd(3,0.5,N,1);
% ind=find(gammatrue==0);
% thetatrue(ind)=zeros(length(ind),1);
% z=normrnd(thetatrue,1); % simulate the test statistics for use
%% Now the generation of simulated data has been completed
disp('number of 0s in gammatrue');
sum(gammatrue==0)
disp('number of 1s in gammatrue');
sum(gammatrue==1)
%% Preparation for performing SMC
ps=1e4; % particle size of importance sampling
std_new_comp=1; % initial standard error for a new added component in f_{1,m}
Prob_alt=.1; % inital probability of the alternative hypothesis
Prob_null=1-Prob_alt; % inital probability of the null hypothesis
N_alt=1;% inital number of data points included in the alternative distribution
N_null=9;% inital number of data points included in the null distribution
Q2=.5;
beta.mean=betatrue;%[-3 1/ 1];%mean(bfdr1.coefficients(1:10:end,:),1);
beta.cov=diag([10 10 10]);%cov(bfdr1.coefficients(1:10:end,:));
betapar=[unifrnd(-10,10,ps,1) unifrnd(-10,10,ps,1) unifrnd(-10,10,ps,1)];%betapar=mvnrnd(beta.mean,beta.cov,ps);
for i=1:ps
M0(i).mu=0; % initialize the mean of the null distribution M0
M0(i).sig=1.5; % initialize the standard error of the null distribution M0
% initialize parameters for the alternative distribution M1
M1(i).weights=1;
M1(i).means=3;
M1(i).variance=20;
M1(i).ncomps=1;
N_alt(i)=1; % inital number of data points included in the alternative distribution
N_null(i)=9; % inital number of data points included in the null distribution
end
lik=ones(ps,1); % initial partice weights
b=(4/((P+1+2)*ps))^(1/(P+1+4)); % kernel width for use in particle rejuvenation operations for \beta
a=sqrt(1-b^2);
%% SMC sampling process begin
for t=1:N
t/N
%% weighting step
for i=1:ps
f0(i) = dnorm(z(t), M0(i).mu, M0(i).sig);
f1(i) = densitynormix(z(t), M0(i).sig.^2, M1(i).weights, M1(i).means, M1(i).variance);
Psi =cbind(1,X(t,:))*betapar(i,:)';
W = ilogit(Psi);
lik(i) = W*f1(i).mix+(1-W)*f0(i);
end
lik=lik/sum(lik);
%% now weighting step is completed
ess(t)=1/sum(lik.^2)/ps; % effective sample size
[~,ind]=max(lik);
beta_est(t,:)=betapar(ind,:); % record the maximum a posterior estimate of \beta
M0_est=M0(ind); % record the maximum a posterior estimate of M0
M1_est=M1(ind); % record the maximum a posterior estimate of M1
%% Resampling step
outIndex = residualR(1:ps,lik);
betapar = betapar(outIndex,:);
M0=M0(outIndex);
M1=M1(outIndex);
N_alt=N_alt(outIndex);
N_null=N_null(outIndex);
f0=f0(outIndex);
f1=f1(outIndex);
%% now the Resampling step is completed
%% Particle rejuvenation operations for beta
mean_particle=mean(betapar,1);
error=betapar-repmat(mean_particle,ps,1);
cov_particle=error'*error/ps;
betapar_tmp=a*betapar+repmat((1-a)*mean_particle,ps,1); % Ref:Fig.3 in A One-Pass Sequential Monte Carlo Method for Bayesian Analysis of Massive Datasets
betapar=betapar_tmp+mvnrnd(zeros(1,P+1),b^2*cov_particle,ps);
%% Particle rejuvenation operations for mixture model parameters
for i=1:ps
% online k-means approximation
Psi2 =cbind(1,X(t,:))*betapar(i,:)';
W2 = ilogit(Psi2);
Pp = W2*f1(i).mix/(W2*f1(i).mix+(1-W2)*f0(i)); % Posterior probability of the alternative
if Pp<1-Q2 % allocate z(t) to the null distribution
rou=1/(1+N_null(i));
M0(i).mu=0;
M0(i).sig=sqrt((1-rou)*M0(i).sig^2+rou*(z(t)-M0(i).mu)^2);
N_null(i)=N_null(i)+1;
else % allocate z(t) to the alternative distribution
alpha=1/(1+N_alt(i));
M1(i).weights=(1-alpha)*M1(i).weights;
dis=zeros(1,M1(i).ncomps);
for j=1:M1(i).ncomps
dis(j)=abs(z(t)-M1(i).means(j))/sqrt(M1(i).variance(j));
end
[min_dis,ind]=min(dis);
if min_dis>2.5 && M1(i).ncomps<3 % z(Ni) does not belong to any component
M1(i).ncomps=M1(i).ncomps+1;
M1(i).means(M1(i).ncomps)=z(t);
M1(i).weights(M1(i).ncomps)=alpha;
M1(i).variance(M1(i).ncomps)=std_new_comp^2;
else
rou=alpha/(M1(i).weights(ind)+alpha);
M1(i).means(ind)=(1-rou)*M1(i).means(ind)+rou*z(t);
M1(i).weights(ind)=M1(i).weights(ind)+alpha;
M1(i).variance(ind)=(1-rou)*M1(i).variance(ind)+rou*(z(t)-M1(i).means(ind))^2;
end
N_alt(i)=N_alt(i)+1;
end
end
end
%% SMC sampling process is completed
%% Bayes testing
Psi =cbind(1,X)*beta_est(end,:)';
W = ilogit(Psi);
f0_full = dnorm(z, M0_est.mu, M0_est.sig*ones(N,1));
f1_full = marnormix(z, M0_est.sig^2, M1_est.weights, M1_est.means, M1_est.variance);
PostProb = W.*f1_full./((1-W).*f0_full + W.*f1_full);
Res.BF=(W.*f1_full)./((1-W).*f0_full);
guess=zeros(N,1);
ind =find(Res.BF>(1-Q2)/Q2);
guess(ind)=ones(length(ind),1);
table(1,1)=length(find(gammatrue==0 & guess==0 ));
table(1,2)=length(find(gammatrue==0 & guess==1 ));% false positives
table(2,1)=length(find(gammatrue==1 & guess==0 ));
table(2,2)=length(find(gammatrue==1 & guess==1 ));% correct detections
correct_findings=find(gammatrue==1 & guess==1 );
wrong_findings=find(gammatrue==0 & guess==1 );
%% show the number of detections and errors
table
%% Plot inference results
figure,
histogram(betapar(:,1), 20,'Normalization','pdf');hold on;
plot([betatrue(1) betatrue(1)],[0 5.5],'r','LineWidth',4);hold off;axis([-4.2 -2.4 0 5.5]);
xlabel('\beta_0'); ylabel('pdf estimate');grid on; hold off;
figure,plot(1:t,beta_est(:,1),1:t,betatrue(1)*ones(1,t),'r','LineWidth',2);grid on;
xlabel('t');ylabel('\beta_0');
figure,plot(1:t,beta_est(:,2),1:t,betatrue(2)*ones(1,t),'r','LineWidth',2);grid on;
xlabel('t');ylabel('\beta_1');
figure,plot(1:t,beta_est(:,3),1:t,betatrue(3)*ones(1,t),'r','LineWidth',2);grid on;
xlabel('t');ylabel('\beta_2');
figure,plot(1:t,ess);xlabel('t');ylabel('NESS');axis([0 10000 0 1.1]);grid on;