function [Run] = PROC_cpBGe_regularized(data, steps, p, k)
global DATA_ALL;
[DATA_ALL] = SHIFT_DATA(data);
global T_0;
global my_0;
global v;
global alpha;
global Prior;
global g_vec;
global g0_vec;
global G_vec;
global G0_vec;
[T_0, my_0, v, alpha, Prior, g_vec, g0_vec, G_vec, G0_vec] = SET_HYPERPARAMETERS(DATA_ALL, p,k);
[Run] = INITIALISE_cpBGe_regularized(DATA_ALL, steps);
[Run] = START_cpBGe_regularized(DATA_ALL, steps, Run);
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [DATA_SHIFTED] = SHIFT_DATA(data)
[n, t] = size(data);
n_plus = n+1; % Covariance matrices will be of size (n+1)x(n+1)
for node_i=1:n
% For each variable X_i consider
% data_1[X_1(t-1),...,X_(n)(t-1),...,X_1(t-slice),...,X_(n)(t-slice), X_(i)(t)]
obs = 1:(t-1);
data_new = zeros(n_plus,t-1);
data_new(1:n,obs) = data(1:n,1:(t-1));
data_new(n+1,obs) = data(node_i,2:t);
DATA_SHIFTED{node_i} = data_new;
end
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [T_0, my_0, v, alpha, Prior, g_vec, g0_vec, G_vec, G0_vec] = SET_HYPERPARAMETERS(DATA_ALL, p_transition, k_success)
[n_plus, n_obs] = size(DATA_ALL{1});
n_nodes = n_plus-1;
Prior = zeros(n_nodes+1,1);
% THESE VALUES MAY BE MODIFIED BY A USER: v, alpha, my_0, and B
v = 1;
alpha = n_plus+2;
my_0 = 0 * ones(n_plus,1);
v_vec = 1 * ones(n_plus,1);
B = zeros(n_plus,n_plus);
% COMPUTE PRECISION-MATRIX W
W = 1/v_vec(1);
for i = 1:(n_plus-1)
b_vec = B(1:i,(i+1));
W = [W + (b_vec * b_vec')/v_vec(i+1) , (-1)*b_vec/v_vec(i+1) ; (-1)*b_vec'/v_vec(i+1), 1/v_vec(i+1)];
end
% COMPUTE T_0:
T_0 = (v*(alpha-n_plus-1))/(v+1) * inv(W);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
g_vec = zeros(1,100);
g0_vec = zeros(1,100);
for t=k_success:100
g_vec(t) = prod(1:(t-1))/(prod(1:(k_success-1))*prod(1:(t-k_success))) * p_transition^(k_success) * (1-p_transition)^(t-k_success);
end
for t=1:100
g0_vec_t = 0;
for i=1:min([t,k_success])
g0_vec_t = g0_vec_t + (1/k_success) * prod(1:(t-1))/(prod(1:(i-1))*prod(1:(t-i))) * p_transition^i * (1-p_transition)^(t-i);
end
g0_vec(1,t) = g0_vec_t;
end
G_vec = cumsum(g_vec);
G0_vec = cumsum(g0_vec);
return;
function [Run] = INITIALISE_cpBGe_regularized(DATA_ALL, steps)
[n_plus, m]= size(DATA_ALL{1});
n_nodes = length(DATA_ALL);
DAG = zeros(n_nodes,n_nodes);
MATRIX = ones(n_nodes,m);
NODE_VEC = 1:n_nodes;
for node_i=1:n_nodes
data = DATA_ALL{node_i};
node_i_comp = NODE_VEC(node_i);
for component=1:max(MATRIX(node_i_comp,:))
DATA{node_i}{component} = data(:,find(MATRIX(node_i_comp,:)==component));
end
end
[log_score] = COMPUTE_LOG_SCORE_MIX(DATA, DAG, MATRIX, NODE_VEC);
for i=1:steps+1
Run.dag{i} = 0;
Run.matrix{i} = 0;
Run.log_score(i) = 0;
Run.Node_Vec{i} = 0;
end
Run.dag{1} = DAG;
Run.matrix{1} = MATRIX;
Run.Node_Vec{1} = NODE_VEC;
Run.steps(1) = 1;
Run.log_score = log_score;
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [log_score] = COMPUTE_LOG_SCORE_MIX(DATA, DAG, MATRIX, NODE_VEC)
% The hyperparameters:
global T_0;
global v;
global alpha;
global my_0;
global g_vec;
global g0_vec;
global G_vec;
global G0_vec;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n_nodes = length(DAG);
counts = hist(NODE_VEC,1:n_nodes);
ind_non_empty = find(counts);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
G_vec_zero = [0,G_vec];
[n_nodes, m]= size(MATRIX);
log_prob_breaks = 0;
for ind=ind_non_empty
break_vec = find(MATRIX(ind,2:end) - MATRIX(ind,1:(end-1)));
if (length(break_vec)==0)
log_prob_breaks_i = log(1-G0_vec(m-1));
elseif (length(break_vec)==1)
log_prob_breaks_i = log(g0_vec(break_vec(1))) + log(1-G_vec_zero(m-break_vec(1)));
else
log_prob_breaks_i = log(g0_vec(break_vec(1))) + log(1-G_vec_zero(m-break_vec(end)));
for i=2:(length(break_vec))
log_prob_breaks_i = log_prob_breaks_i + log(g_vec(break_vec(i)-break_vec(i-1)));
end
end
log_prob_breaks = log_prob_breaks + log_prob_breaks_i ;
end
log_prob_data = 0;
% Compute the local score for each node:
for i_node=1:n_nodes
k_i = length(DATA{i_node}); % k_i is the number of mixture components for the i-th node
log_score_node_i = 0;
for component=1:k_i
[n_plus, n_obs_in_component] = size(DATA{i_node}{component});
if(n_obs_in_component==0)
% do nothing
else
data = DATA{i_node}{component};
% Compute T_m for the current combination of node and
% component:
parents_of_node = find(DAG(:,i_node));
parents_and_node = [parents_of_node; n_plus];
S_m = (n_obs_in_component-1) * cov(data(parents_and_node,:)');
T_m = T_0(parents_and_node,parents_and_node) + S_m + (v*n_obs_in_component)/(v+n_obs_in_component) * (my_0(parents_and_node,1) - mean(data(parents_and_node,:)',1)') * (my_0(parents_and_node,1) - mean(data(parents_and_node,:)',1)')';
[log_score_local_node_nominator] = Gauss_Score_complete(length(parents_and_node), n_obs_in_component, v, alpha, T_0(parents_and_node,parents_and_node), T_m);
[log_score_local_node_denominator] = Gauss_Score_complete(length(parents_of_node), n_obs_in_component, v, alpha, T_0(parents_of_node,parents_of_node), T_m(1:(end-1),1:(end-1)));
log_score_node_i = log_score_node_i + (log_score_local_node_nominator - log_score_local_node_denominator);
end
end
log_prob_data = log_prob_data + log_score_node_i;
end
log_score = log_prob_breaks + log_prob_data;
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [log_Score_i] = Gauss_Score_complete(n_nodes_in_sub, n_obs, v, alpha, T_0_sub, T_m_sub)
if n_nodes_in_sub==0
log_Score_i = 0;
else
sum_1 = (-1)*(n_nodes_in_sub)*(n_obs)/2 * log(2*pi) + (n_nodes_in_sub/2) * log(v/(v+n_obs));
sum_3 = (alpha/2) * log(det(T_0_sub)) + (-1)*(alpha+n_obs)/2 * log(det(T_m_sub));
log_value_a = (alpha *n_nodes_in_sub/2)*log(2) + (n_nodes_in_sub*(n_nodes_in_sub-1)/4)*log(pi);
log_value_b = ((alpha+n_obs) *n_nodes_in_sub/2)*log(2) + (n_nodes_in_sub*(n_nodes_in_sub-1)/4)*log(pi);
for i=1:n_nodes_in_sub
log_value_a = log_value_a + gammaln(( alpha +1-i)/2);
log_value_b = log_value_b + gammaln(((alpha+n_obs)+1-i)/2);
end
log_Score_i = sum_1 + (-1)*(log_value_a - log_value_b) + sum_3;
end
return
function [Run] = START_cpBGe_regularized(DATA_ALL, steps, Run)
last = Run.steps(1); % index of last sampled dags
DAG = Run.dag{last}; % the current DAG
MATRIX = Run.matrix{last}; % the current allocation vector
NODE_VEC = Run.Node_Vec{last};
% Output on screen:
fprintf('\n###########################################################\n')
fprintf('The Gibbs simulation for the regularized cpBGe model has been started \n')
fprintf('###########################################################\n')
CONSTANT = exp(300);
n_nodes = length(DAG);
counts = hist(NODE_VEC,1:n_nodes);
ind_non_empty = find(counts);
DAG_full = ones(n_nodes,n_nodes);
for ind=ind_non_empty
nodes_in_comp = find(NODE_VEC==ind);
[Q_VEC] = Dynamic_Q_MIX(DATA_ALL, nodes_in_comp, DAG_full, CONSTANT);
[Vector] = Simulate_changepoints_MIX(Q_VEC, DATA_ALL, nodes_in_comp, DAG_full, CONSTANT);
MATRIX(ind,:) = Vector;
end
%%% Start of the MCMC simulation
for i = last:(steps)
i
[DAG, MATRIX, NODE_VEC, log_score] = MCMC_INNER_MIX(DAG, MATRIX, NODE_VEC, DATA_ALL, CONSTANT);
Run.dag{i+1} = DAG;
Run.matrix{i+1} = MATRIX;
Run.Node_Vec{i+1} = NODE_VEC;
Run.steps(1) = i+1;
Run.log_score(i+1) = log_score;
end
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [DAG, MATRIX, NODE_VEC, log_score] = MCMC_INNER_MIX(DAG, MATRIX, NODE_VEC, DATA_ALL, CONSTANT)
[thrash, n_nodes] = size(DAG);
for child_node=1:n_nodes
data = DATA_ALL{child_node};
clear DATA;
child_node_comp = NODE_VEC(child_node);
for component=1:max(MATRIX(child_node_comp,:))
DATA{child_node}{component} = data(:,find(MATRIX(child_node_comp,:)==component));
end
[SCORES] = SCORE_ALL_PARENTS(DATA, child_node, n_nodes);
cumulative = cumsum(SCORES.exp_scores)/sum(SCORES.exp_scores);
x = rand;
indicis = find(cumulative>=x);
index = indicis(1);
new_parents_of_node = nonzeros(SCORES.parents(index,:));
DAG(:,child_node)=0;
DAG(new_parents_of_node,child_node)=1;
clear DATA;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
counts = hist(NODE_VEC,1:n_nodes);
ind_non_empty = find(counts);
for ind=ind_non_empty
nodes_in_comp = find(NODE_VEC==ind);
[Q_VEC] = Dynamic_Q_MIX(DATA_ALL, nodes_in_comp, DAG, CONSTANT);
[Vector] = Simulate_changepoints_MIX(Q_VEC, DATA_ALL, nodes_in_comp, DAG, CONSTANT);
MATRIX(ind,:) = Vector;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BIRTH/DEATH MOVES
P_Death=0.5;
for steps=1:10
counts = hist(NODE_VEC,1:n_nodes);
ind_eq1 = find(counts==1);
ind_gr1 = find(counts>1);
ind_geq1 = find(counts>=1);
n_eq1 = length(ind_eq1);
n_gr1 = length(ind_gr1);
n_geq1 = length(ind_geq1);
x = rand;
if (x<=P_Death) % -> Death move
if (n_eq1==0)
% do nothing
else
indicis_1 = randperm(n_eq1); % the single node with label_1 WILL GET label_2 -> i.e. the component with label_1 "dies"
label_1 = ind_eq1(indicis_1(1));
indicis_2 = randperm(n_geq1);
label_2 = ind_geq1(indicis_2(1));
if (label_1==label_2)
label_2 = ind_geq1(indicis_2(2)); % this ensures that label_1 is different from label_2
end
single_node_with_label_1 = find(NODE_VEC==label_1); % the node with label_1
nodes_with_label_2 = find(NODE_VEC==label_2); % the node(s) with label_2
NODE_VEC_NEW = NODE_VEC;
NODE_VEC_NEW(single_node_with_label_1) = label_2; % the single node with label_1 has just obtained label_2
nodes_with_label_2_new = find(NODE_VEC_NEW==label_2); % the nodes with label_2 after death move
[Q_new] = Dynamic_Q_MIX(DATA_ALL, nodes_with_label_2_new, DAG, CONSTANT);
[Vector_2] = Simulate_changepoints_MIX(Q_new, DATA_ALL, nodes_with_label_2_new, DAG, CONSTANT);
MATRIX_NEW = MATRIX;
MATRIX_NEW(label_2,:) = Vector_2;
counts_new = hist(NODE_VEC_NEW,1:n_nodes);
ind_gr1_new = find(counts_new>1);
n_gr1_new = length(ind_gr1_new);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[Q_old1] = Dynamic_Q_MIX(DATA_ALL, single_node_with_label_1, DAG, CONSTANT);
[Q_old2] = Dynamic_Q_MIX(DATA_ALL, nodes_with_label_2, DAG, CONSTANT);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
log_forward = log(1/n_eq1) + log(1/(n_geq1-1)) - log(Q_new(1));
log_backward = log(1/n_gr1_new) + log(1/length(nodes_with_label_2_new)) - log(Q_old1(1)) - log(Q_old2(1)) + log(CONSTANT);
log_hastings = log_backward - log_forward;
A = exp(log_hastings);
y = rand;
if (y Birth move
if(n_gr1==0)
% do nothing
else
indicis_1 = randperm(n_gr1); % one node with label_1, namely node_out, WILL GET the label_2
label_1 = ind_gr1(indicis_1(1));
ind_empty = find(counts==0);
label_2 = ind_empty(1);
nodes_with_label_1_old = find(NODE_VEC==label_1);
n_nodes_with_label_1_old = length(nodes_with_label_1_old);
indicis = randperm(n_nodes_with_label_1_old);
node_out = nodes_with_label_1_old(indicis(1));
NODE_VEC_NEW = NODE_VEC;
NODE_VEC_NEW(node_out) = label_2;
nodes_with_label_1_new = find(NODE_VEC_NEW==label_1);
[Q_new1] = Dynamic_Q_MIX(DATA_ALL, nodes_with_label_1_new, DAG, CONSTANT);
[Vector_1] = Simulate_changepoints_MIX(Q_new1, DATA_ALL, nodes_with_label_1_new, DAG, CONSTANT);
[Q_new2] = Dynamic_Q_MIX(DATA_ALL, node_out, DAG, CONSTANT);
[Vector_2] = Simulate_changepoints_MIX(Q_new2, DATA_ALL, node_out, DAG, CONSTANT);
MATRIX_NEW = MATRIX;
MATRIX_NEW(label_1,:) = Vector_1;
MATRIX_NEW(label_2,:) = Vector_2;
counts_new = hist(NODE_VEC_NEW,1:n_nodes);
ind_eq1_new = find(counts_new==1);
n_eq1_new = length(ind_eq1_new);
ind_geq1_new = find(counts_new>=1);
n_geq1_new = length(ind_geq1_new);
[Q_old1] = Dynamic_Q_MIX(DATA_ALL, nodes_with_label_1_old, DAG, CONSTANT);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
log_forward = log(1/n_gr1) + log(1/n_nodes_with_label_1_old) - log(Q_new1(1)) - log(Q_new2(1)) + log(CONSTANT);
log_backward = log(1/n_eq1_new) + log(1/n_geq1_new) - log(Q_old1(1));
log_hastings = log_backward - log_forward;
A = exp(log_hastings);
y = rand;
if (y1);
ind_geq1 = find(counts>=1);
n_gr1 = length(ind_gr1);
n_geq1 = length(ind_geq1);
if (n_gr1==0 | n_geq1==1)
% do nothing
else
indicis_1 = randperm(n_gr1); % one node with label_1, namely node_out, WILL GET the label_2
label_1 = ind_gr1(indicis_1(1));
indicis_2 = randperm(n_geq1);
label_2 = ind_geq1(indicis_2(1));
if (label_1==label_2)
label_2 = ind_geq1(indicis_2(2));
end
nodes_with_label_1_old = find(NODE_VEC==label_1);
nodes_with_label_2_old = find(NODE_VEC==label_2);
n_nodes_with_label_1_old = length(nodes_with_label_1_old);
indicis = randperm(n_nodes_with_label_1_old);
node_out = nodes_with_label_1_old(indicis(1));
NODE_VEC_NEW = NODE_VEC;
NODE_VEC_NEW(node_out) = label_2;
nodes_with_label_1_new = find(NODE_VEC_NEW==label_1);
nodes_with_label_2_new = find(NODE_VEC_NEW==label_2);
[Q_new1] = Dynamic_Q_MIX(DATA_ALL, nodes_with_label_1_new, DAG, CONSTANT);
[Vector_1] = Simulate_changepoints_MIX(Q_new1, DATA_ALL, nodes_with_label_1_new, DAG, CONSTANT);
[Q_new2] = Dynamic_Q_MIX(DATA_ALL, nodes_with_label_2_new, DAG, CONSTANT);
[Vector_2] = Simulate_changepoints_MIX(Q_new2, DATA_ALL, nodes_with_label_2_new, DAG, CONSTANT);
MATRIX_NEW = MATRIX;
MATRIX_NEW(label_1,:) = Vector_1;
MATRIX_NEW(label_2,:) = Vector_2;
counts_new = hist(NODE_VEC_NEW,1:n_nodes);
ind_gr1_new = find(counts_new>1);
n_gr1_new = length(ind_gr1_new);
ind_geq1_new = find(counts_new>=1);
n_geq1_new = length(ind_geq1_new);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[Q_old1] = Dynamic_Q_MIX(DATA_ALL, nodes_with_label_1_old, DAG, CONSTANT);
[Q_old2] = Dynamic_Q_MIX(DATA_ALL, nodes_with_label_2_old, DAG, CONSTANT);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
log_forward = log(1/n_gr1) + log(1/(n_geq1-1)) - log(Q_new1(1)) - log(Q_new2(1));
log_backward = log(1/n_gr1_new) + log(1/(n_geq1_new-1))- log(Q_old1(1)) - log(Q_old2(1));
log_hastings = log_backward - log_forward;
A = exp(log_hastings);
y = rand;
if (y [1,2,3,...,m_new]
%2nd column: cp at t=1 -> [2,3,4,...,m_new]
%(m_new)-th column: cp at t=m_new-1 -> [m_new]
Vector(:,1) = 1; % changpoint at t=0
for t=1:(m_new-1) % corresponds to possible cps at 0:(m_new-2)
% Note that m_new-1 is the LAST possible cp anyway
n_t = 0;
cps_at_t = find(Vector(1,:)==1);
last_cp = cps_at_t(end); % last_cp out of {1,...,m_new}
if (last_cp==t)
n_t = n_t+1;
end
if (n_t==1)
[Prob] = compute_prob_MIX(Q_vec, DATA_ALL, nodes_in_comp, DAG, t, CONSTANT);
Prob = Prob/sum(Prob); % I guess this normalisation is required!!!!
index_no_cp = length(Prob); % check will refer to the last probability in Prob -> no further changepoint
[OUTPUT] = CARPENTER(1,Prob); % A sample of size 1 from the discrete distribution Prob with realisations in {1,...,m_new)
% 1 means cp at 1 -> Vector(1,2)
% 2 means cp at 2 -> Vector(1,3)
% m_new-1 means cp at m_new-1 ->
% Vector(1,m_new)
% m_new means no further changepoint
indicis_for_cp = find(OUTPUT~=index_no_cp);
OUTPUT = OUTPUT(indicis_for_cp);
for j=1:length(OUTPUT)
Vector(1,OUTPUT(j)+1)=1; % OUTPUT(j)+1 because the first column of Matrix corresponds to cp at 0
end
end
end
Vector = cumsum(Vector,2);
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Prob] = compute_prob_MIX(Q_vec, DATA_ALL, nodes_in_comp, DAG, t, CONSTANT)
% t = 1,2,3,....,m_new-1
% corresponds to the old changepoint
% tau_old = 0,1,2,...,m_new-2
global g_vec;
global g0_vec;
global G_vec;
global G0_vec;
tau_old = t-1;
data_1 = DATA_ALL{1};
[n,m_new] = size(data_1);
Prob = zeros(1,m_new); % success time 1,2,3,...,m_new-1 AND m_new means no further changepoints
if (tau_old>0)
for tau=(tau_old+1):(m_new-1)
Prob(1,tau) = exp(COMPUTE_LOCAL_LOG_SCORE_MIX(DATA_ALL, nodes_in_comp, (tau_old+1):tau, DAG)+log(Q_vec(tau+1))+log(g_vec(tau-tau_old))-log(Q_vec(tau_old+1)));
end
Prob(1,m_new) = exp(COMPUTE_LOCAL_LOG_SCORE_MIX(DATA_ALL, nodes_in_comp, (tau_old+1):m_new, DAG)+log(1-G0_vec(m_new-tau_old-1))+log(CONSTANT)-log(Q_vec(tau_old+1)));
else % tau_old==0
for tau=1:(m_new-1)
Prob(1,tau) = exp(COMPUTE_LOCAL_LOG_SCORE_MIX(DATA_ALL, nodes_in_comp, 1:tau, DAG)+log(Q_vec(tau+1))+log(g0_vec(tau))-log(Q_vec(1)));
end
Prob(1,m_new) = exp(COMPUTE_LOCAL_LOG_SCORE_MIX(DATA_ALL, nodes_in_comp, 1:m_new, DAG)+log(1-G0_vec(m_new-1))+log(CONSTANT)-log(Q_vec(1)));
end
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [log_score] = COMPUTE_LOCAL_LOG_SCORE_MIX(DATA_ALL, nodes_in_comp, interval, DAG)
global T_0;
global v;
global alpha;
global my_0;
log_score = 0;
for node = nodes_in_comp
data_node = DATA_ALL{node};
data_node = data_node(:,interval);
[n_plus,m_new] = size(data_node);
parents = find(DAG(:,node));
parents_and_node = [parents; n_plus];
data_node = data_node(parents_and_node,:);
S_m = (m_new-1) * cov(data_node');
T_m = T_0(parents_and_node,parents_and_node) + S_m + (v*m_new)/(v+m_new) * (my_0(parents_and_node,1) - mean(data_node',1)') * (my_0(parents_and_node,1) - mean(data_node',1 )')';
[log_score_local_nominator] = Gauss_Score_complete(length(parents_and_node), m_new, v, alpha, T_0(parents_and_node,parents_and_node), T_m);
[log_score_local_denominator] = Gauss_Score_complete(length(parents), m_new, v, alpha, T_0(parents,parents) , T_m(1:(end-1),1:(end-1)));
log_score = log_score + (log_score_local_nominator - log_score_local_denominator);
end
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Q_VEC] = Dynamic_Q_MIX(DATA_ALL, nodes_in_comp, DAG, CONSTANT)
data_1 = DATA_ALL{1};
[n,m_new] = size(data_1); % NOTE: m_new = m-1
Q_VEC = zeros(1,m_new); % Q_VEC(1,m_new) corresponds to Q(m_new)=P(m_new:m_new|cp at m_new-1)
% Q_VEC(1,1) corresponds to Q(1)=P(1:m_new| cp at 0)
% Determine Q(m_new)
[BN_score] = COMPUTE_LOCAL_LOG_SCORE_MIX(DATA_ALL, nodes_in_comp, m_new:m_new, DAG);
Q_VEC(1,m_new) = exp(BN_score +log(CONSTANT));
for i=1:(m_new-1)
t = m_new-i;
% i=1 gives t=m_new-1
% i=m_new-1 gives t=1
Q_t = RECURSION_MIX(Q_VEC, t, DATA_ALL, nodes_in_comp, DAG, CONSTANT);
Q_VEC(1,t) = Q_t;
end
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Q_t] = RECURSION_MIX(Q_VEC, t, DATA_ALL, nodes_in_comp, DAG, CONSTANT)
global g_vec;
global g0_vec;
global G_vec;
global G0_vec;
% t ranges from: m_new-1,...,1
if (t>1)
data_1 = DATA_ALL{1};
[n, m_new] = size(data_1);
Q_t = 0;
for s=t:(m_new-1)
Q_t = Q_t + exp(COMPUTE_LOCAL_LOG_SCORE_MIX(DATA_ALL, nodes_in_comp, t:s, DAG)+log(Q_VEC(1,s+1))+log(g_vec(s+1-t)));
end
Q_t = Q_t + exp(COMPUTE_LOCAL_LOG_SCORE_MIX(DATA_ALL, nodes_in_comp, t:m_new, DAG)+log(1-G_vec(m_new-t))+log(CONSTANT));
else % t==1
data_1 = DATA_ALL{1};
[n, m_new] = size(data_1);
Q_t = 0;
for s=1:(m_new-1)
Q_t = Q_t + exp(COMPUTE_LOCAL_LOG_SCORE_MIX(DATA_ALL, nodes_in_comp, 1:s, DAG)+log(Q_VEC(1,s+1))+log(g0_vec(s)));
end
Q_t = Q_t + exp(COMPUTE_LOCAL_LOG_SCORE_MIX(DATA_ALL, nodes_in_comp, 1:m_new, DAG)+log(1-G0_vec(m_new-1))+log(CONSTANT));
end
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [OUTPUT] = CARPENTER(n,Prob)
x_vec = exprnd(1,1,n+1);
S = sum(x_vec);
u_vec = cumsum(x_vec/S);
Q = 0;
U = u_vec(1,1);
j = 1;
i = 1;
OUTPUT = [];
while (i