-
Notifications
You must be signed in to change notification settings - Fork 5
/
Figure4.m
57 lines (42 loc) · 1.7 KB
/
Figure4.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
% Matlab code to make Figure 4 for "Spike-centered jitter can mistake temporal structure" (Platkiewicz, Stark, Amarasingham)
% Illustrates non-uniformity of randomized p-values induced by Poisson approximation.
% (C) Asohan Amarasingham, 6/5/2016
N=500; % Binomial
p=.1; % parameters
lambda=N*p; % for Poisson approximation
m=50000; % number of samples
binw=0.01;
clf
for j=1:m
% Sample Binomial(N,p)
S = sum( rand(N,1)<=p );
% compute p-value (poiss approx.)
pvalp(j)=1-poisscdf(S-1,lambda);
% compute randomized p-value (poiss approx.)
pvalp_rand(j)=[1-poisscdf(S-1,lambda)] - rand(1)*poisspdf(S,lambda);
% compute p-value (exact binomial)
pvalb(j)=1-binocdf(S-1,N,p);
% compute randomized p-value (binomial xact)
pvalb_rand(j)=[1-binocdf(S-1,N,p)] - rand(1)*binopdf(S,N,p);
%if mod(j,1000)==0, j, end
end
subplot(2,2,1)
histogram(pvalp,0:binw:1,'Normalization','probability')
hold on
plot(0:.005:1,binw*ones( size(0:.005:1)),'r-.') % raw poiss approx. pvalues
title('(a) Raw p-values (Poisson approximation method)')
subplot(2,2,3)
histogram(pvalp_rand,0:binw:1,'Normalization','probability')
hold on
plot(0:.005:1,binw*ones( size(0:.005:1)),'r-.') % randomized poiss approx. pvalues
title('(b) Randomized p-values (Poisson approximation)')
subplot(2,2,2)
histogram(pvalb,0:binw:1,'Normalization','probability')
hold on
plot(0:.005:1,binw*ones( size(0:.005:1)),'r-.') % raw binomial (exact) pvalues
title('(c) Exact p-values (Binomial)')
subplot(2,2,4)
histogram(pvalb_rand,0:binw:1,'Normalization','probability')
hold on
plot(0:.005:1,binw*ones( size(0:.005:1)),'r-.') % randomized binomial exact pvalues
title('(d) Randomized exact p-values (Binomial)')