I have this error and I got stoped, please help me: Attempted to access indicis(2); index out of bounds because numel(indicis)=1. Error in Dag_MCMC002 (line 120), index = indicis(2);

function [Run] = Dag_MCMC002(data,dist,steps, DAG, fan_in, vector, k_max, lambda_snr_vec, alpha_sig, beta_sig)

[n_parents, n_nodes] = size(DAG);
 DATA_ALL = SHIFT_DATA(data);
[n_plus, m]  = size(DATA_ALL); 
Counter=1;
for t=1:steps

for i_node=1:n_nodes

      DATA     = [];
      y2       = rand(1);

         if (y2<0.5) % Perform an addition/deletion move

              child_node    = i_node;

              old_parents = DAG(:,child_node);
              new_parents = old_parents;

              parents_card_old = sum(old_parents);

              parents_indicis_old = find(old_parents);

                  if (parents_card_old<fan_in);

                          neighbours_old = n_parents-1; 

                          indicis = randperm(n_parents);
                          x = indicis(1);

                          if (x==child_node)
                              x = indicis(2);
                          end

                          parent_index = x; % delete or add this parent node

                          new_parents(parent_index,1) = 1 - new_parents(parent_index,1); 

                  else % elseif (parent_card_old==fan_in)

                      x = randperm(fan_in);
                      x = x(1);

                      parent_index = parents_indicis_old(x); % delete this parent node

                      new_parents(parent_index,1) = 0;

                      neighbours_old = parents_card_old; % = fan_in
                  end

                      parents_card_new    = sum(new_parents);

                      if (parents_card_new<fan_in)         
                      neighbours_new = n_parents-1; 
                      else
                      neighbours_new = parents_card_new; % = fan_in
                      end

                      DAG_NEW = DAG;
                      DAG_NEW(:,child_node) = new_parents;

                          data = DATA_ALL{child_node};
                          k_i = max(vector);
                          for component=1:k_i  
                              DATA{child_node}{component} = data(:,find(vector==component));
                          end

                  local_score_new = COMPUTE_LOCAL_LOG_SCORE_REGRESS(DATA, DAG_NEW, vector, child_node, alpha_sig, beta_sig, lambda_snr_vec);
                  local_score_old = COMPUTE_LOCAL_LOG_SCORE_REGRESS(DATA, DAG,     vector, child_node, alpha_sig, beta_sig, lambda_snr_vec);

                  A = exp(local_score_new - local_score_old + log(neighbours_old) - log(neighbours_new));

                  u = rand(1); 

                  if u > min(1,A) % then reject the movw
                  else % accept the move:
                  DAG = DAG_NEW;
                  end

                  clear DATA;

         else % Perform an exchange move

                  child_node    = i_node;

                  old_parents = find(DAG(:,child_node));

                  parents_card_old = length(old_parents);

                      if (parents_card_old==0);

                          % do nothing

                      else % perform an exchange move

                          DAG_NEW = DAG;

                          indicis        = randperm(parents_card_old);
                          index          = indicis(1);

                          parent_old_index = old_parents(index); % delete this parent node

                          candidate_parents = find(DAG(:,child_node)==0);

                          candidates_card = length(candidate_parents);

                          indicis = randperm(candidates_card);
                          index   = indicis(1);

                          parent_new_index = candidate_parents(index);                        DAG_NEW(parent_new_index,child_node) = 1;

                              data = DATA_ALL{child_node};
                              k_i = max(vector);

                              for component=1:k_i      
                                  DATA{child_node}{component} = data(:,find(vector==component));
                              end

                          local_score_new = COMPUTE_LOCAL_LOG_SCORE_REGRESS(DATA, DAG_NEW, vector, child_node, alpha_sig, beta_sig, lambda_snr_vec);
                          local_score_old = COMPUTE_LOCAL_LOG_SCORE_REGRESS(DATA, DAG,     vector, child_node, alpha_sig, beta_sig, lambda_snr_vec);

                          log_hastings = 0;

                          A = exp(local_score_new - local_score_old + log_hastings);

                          u = rand(1); 

                          if u > min(1,A) % then reject the move:

                          else % accept the move:
                              DAG = DAG_NEW; 
                          end

                          clear DATA;
                      end

         end
  end
       if (mod(t,dist)==1)
       Run.Dag{Counter}       =  DAG;
       Counter = Counter+1;
       end 
      end        
  return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [log_score] = COMPUTE_LOCAL_LOG_SCORE_REGRESS(DATA, Dag, vector, child, alpha_sig, beta_sig, lambda_snr_vec)

global Prior;
m=6;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%log_prob_graph =  Prior(length(find(Dag(:,child)))+1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Compute the local score for child:

Kg =max(vector); 

parents = find(Dag(:,child));

lambda_snr  = lambda_snr_vec(child,1);

sum_log_det_Sigma_tilde = 0;
sum_Delta2              = 0;

     for component=1:Kg

      data = DATA{child}{component};
      [n_plus, obs] =  size(data);

     X{component}= [ones(1,length(find(vector==component)));data(parents,:)];
     y{component}= data(child,:)';
     sum_Delta2 = sum_Delta2 +y{component}'*inv(eye(length(find(vector==component)))+lambda_snr* X{component}'* X{component})*y{component};

     log_det_Sigma_tilde = log(det(eye(length(find(vector==component)))+lambda_snr* (X{component}'* X{component})));
     sum_log_det_Sigma_tilde = sum_log_det_Sigma_tilde + log_det_Sigma_tilde;

     end 
     log_gamma_ratio = gammaln((m-1)/2+alpha_sig) - gammaln(alpha_sig);
     sum_1 = (alpha_sig)*log(2*beta_sig) - ((m-2)/2)*log(pi) - 0.5 * sum_log_det_Sigma_tilde;
     sum_2 = -1* ((m-1)/2+ alpha_sig)*log(2*beta_sig+sum_Delta2);
     log_score= log_gamma_ratio + sum_1 +  sum_2;  

%log_prob_graph

       log_score =  log_score;
  return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [DATA_SHIFTED] = SHIFT_DATA(data)

[n, t] = size(data);

n_plus = n+1; 

for node_i=1:n

              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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Answer 1

[Walter Roberson]

If you pass in a row vector for the DAG then size() of it will be 1 x something, and n_parents will become 1. randperm(n_parents) is going to be of length 1, but you then attempt to access the second element of that.