You are an expert in writing code in many different languages. Your goal is to perform code completion for the following code, keeping in mind the rest of the code and the file meta data.
    Metadata details:
    filename: pvalue2.m
    path of file: ./repos/BRAPH-2/braph2genesis/src/analysis
    
    the code of the file until where you have to start completion: function P2 = pvalue2(observed_difference, random_differences)
% PVALUE2 Calculates the two-tailed p-value
%
% P2 = PVALUE2(OBSERVED_DIFFERENCE, RANDOM_DIFFERENCES) calculates the
% two-tailed P-value of OBSERVED_DIFFERENCE with respect to the
% distribution given by RANDOM_DIFFERENCES.
%
% This function admits various kinds of inputs:
%
% 1) OBSERVED_DIFFERENCE is a scalar, a vector or a matrix with the difference.
% RANDOM_DIFFERENCES must be a cell array of samples, where each cell contains a
% scalar, a vector or a matrix with random variables.
% P1 is a scalar, a vector or a matrix with the two-sided p-values.
%
% 2) OBSERVED_DIFFERENCE is a cell containing a scalar, a vector or a
% matrix with the difference.
% RANDOM_DIFFERENCES must be a cell array of samples, where each cell contains a
% scalar, a vector or a matrix with random variables.
% P1 is a cell containing a scalar, a vector or a matrix with the two-sided
% p-valu
    
function P2 = for_pvalue2(observed_difference, random_differences)
        
        M = numel(random_differences);  %#ok<NASGU> % number of samples (per variable)
        
        row_number = size(random_differences{1}, 1);
        column_number = size(random_differences{1}, 2);
        d_number = size(random_differences{1}, 3);  % add the 3rd dimension 
           
        Q = cell(row_number, column_number); %#ok<NASGU>
        P2 = zeros(row_number, column_number, d_number);
        for i = 1:1:row_number
            for j = 1:1:column_number
                for k = 1:1:d_number
                    current_observed_difference = observed_difference(i, j, k);
                    current_random_differences = cellfun(@(x) x(i, j, k), random_differences);

                    P2(i, j, k) =  ...
                        length(find(abs(current_random_differences) > abs(current_observed_difference))) ...
                        / ...
                        length(current_random_differences);
                    
                    if P2(i, j, k) == 0
                        P2(i, j, k) = .5 * length(current_random_differences)^-1;
                    end