machine_learn/conditional_probability.m at main · MutluEray/machine_learn · GitHub

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
%% Conditional Probability
%% Generate two long-spike time series
N = 10000;
spikeDur  = 10; % a.u. but must be an even number
spikeNumA = .01; % in proportion of total number of points
spikeNumB = .05; % in proportion of total number of points

% Initialize to zeros
[spike_tsA,spike_tsB] = deal( zeros(N,1) );

% Populate time series 1
spiketimesA = randi(N,round(N*spikeNumA),1);

% Flesh out spikes (loop per spike)
for spikei=1:length(spiketimesA)
    % Find boundaries
    bnd_pre = max(1,spiketimesA(spikei)-spikeDur/2);
    bnd_pst = min(N,spiketimesA(spikei)+spikeDur/2);

    % Fill in with ones
    spike_tsA(bnd_pre:bnd_pst) = 1;
end

% Repeat for time series 2
spiketimesB = randi(N,round(N*spikeNumB),1);
% spiketimesB(1:length(spiketimesA)) = spiketimesA; % induce strong conditional probability

% Flesh out spikes (loop per spike)
for spikei=1:length(spiketimesB)
    % Find boundaries
    bnd_pre = max(1,spiketimesB(spikei)-spikeDur/2);
    bnd_pst = min(N,spiketimesB(spikei)+spikeDur/2);

    % Fill in with ones
    spike_tsB(bnd_pre:bnd_pst) = 1;
end

%% Let's see what they look like
figure(1)
plot(1:N,spike_tsA, 1:N,spike_tsB)
set(gca,'ylim',[0 1.2])

%% Compute their probabilities and intersection
% Probabilities
probA = sum(spike_tsA==1) / N;
probB = mean(spike_tsB);

% Joint probability
probAB = mean(spike_tsA+spike_tsB==2);

%% Compute the conditional probabilities
% p(A|B)
pAgivenB = probAB/probB;

% p(B|A)
pBgivenA = probAB/probA;

% Print a little report
disp('  ')
disp([ 'P(A)   = ' num2str(probA) ])
disp([ 'P(A|B) = ' num2str(pAgivenB) ])
disp([ 'P(B)   = ' num2str(probB) ])
disp([ 'P(B|A) = ' num2str(pBgivenA) ])

%% end.