MATLAB中的单元阵列平均值

Question

i对线性模型y = x*beta eps进行仿真研究，大小(x(= [n d]。我认为基于两种方法的维度D的效果。我运行10个模拟数据并获得相应的beta估计，然后我想在10个模拟数据上计算Beta的平均值。

我的玩具Matlab代码如下：

        nsim=10;   %iteration number
        dd=[4 6];  %two dimension cases,beta=(beta_1,cdots,beta_d)^T
        ddlen=length(dd);
        nmethod=2; %two methods
        seednum=0;
        BH  = cell(nsim,ddlen,nmethod); %estimation of beta w.r.t two dimension cases and two methods
        for di = 1:ddlen
            d = dd(di);
            for ni = 1:nsim
                seednum = seednum + di*ni;
                randn('seed', seednum);
                betahat=randn(d,1); 
                for method = 1:nmethod
                    if method==1
                        BH{ni,di,method} = betahat;
                    else
                        BH{ni,di,method} = 10*betahat;
                    end
                end
            end
        end

然后我们可以获得

BH(:,:,1) = 
    [4x1 double]    [6x1 double]
    [4x1 double]    [6x1 double]
    [4x1 double]    [6x1 double]
    [4x1 double]    [6x1 double]
    [4x1 double]    [6x1 double]
    [4x1 double]    [6x1 double]
    [4x1 double]    [6x1 double]
    [4x1 double]    [6x1 double]
    [4x1 double]    [6x1 double]
    [4x1 double]    [6x1 double]

BH(:,:,2) = 
    [4x1 double]    [6x1 double]
    [4x1 double]    [6x1 double]
    [4x1 double]    [6x1 double]
    [4x1 double]    [6x1 double]
    [4x1 double]    [6x1 double]
    [4x1 double]    [6x1 double]
    [4x1 double]    [6x1 double]
    [4x1 double]    [6x1 double]
    [4x1 double]    [6x1 double]
    [4x1 double]    [6x1 double]

我想要10行(NSIM = 10(的平均值，并得到类似

的东西

mean(BH(:,:,1))= 
    [4x1 double]    [6x1 double]
mean(BH(:,:,2)) = 
    [4x1 double]    [6x1 double]

有什么想法吗？谢谢！

Answer 1

我不知道这是最有效的方法，但是您可以使用arrayfun：

% generate random array
BH = repmat({rand(4,1),rand(6,1)},[10 1 2]);
% generate indexes for the 2nd and 3rd dimensions
[n2,n1] = meshgrid(1:size(BH,2),1:size(BH,3));
% get mean across 1st (cell) dimension
[res] = arrayfun(@(n1,n2)mean([BH{:,n1,n2}],2),n1(:),n2(:),'UniformOutput',false);
% reshape to desired output
res = reshape(res,[1 size(BH,2) size(BH,3)]);

如果您想概括为n维单元格数组：

% generate random array
BH = repmat({rand(4,1),rand(6,1)},[10,1,2,2,5]);
sz = size(BH);
% generate indexes for the 2nd and 3rd dimensions
n = cell(1,numel(sz) - 1);
[n{:}] = ndgrid(1:sz(2),1:sz(3),1:sz(4),1:sz(5));
n = cell2mat(cellfun(@(x) {x(:)},n));
idx = 1:size(n,1);
% get mean across 1st (cell) dimension
[res] = arrayfun(@(idx)mean([BH{:,n(idx,1),n(idx,2),n(idx,3),n(idx,4)}],2),...
    idx,'UniformOutput',false);
% reshape to desired output
res = reshape(res,[1 sz(2:end)]);

Answer 2

替代，

% split into seperate cell arrays
BH_1 = BH(:,:,1);
BH_2 = BH(:,:,2);
% create matrix of compatible vectors, and take mean and put result back into cell array
BH_1_mean = cat(2,{mean(cell2mat(BH_1(:,1)'),2)}, {mean(cell2mat(BH_1(:,2)'),2)});
BH_2_mean = cat(2,{mean(cell2mat(BH_2(:,1)'),2)}, {mean(cell2mat(BH_2(:,2)'),2)});

Answer 3

如果我让您正确，您想在矢量中相同位置处的所有元素上均值。因此，从BH(:,1,1)中的所有向量中，我们获得了4均值的一个向量，每个向量都用于向量中的一个位置。BH(:,1,2)也是如此。对于BH(:,2,1)和BH(:,2,1)，我们执行了相同的操作，但是向量中有6个元素。

您可以使用以下代码：

% split BH to 2 arrays:
bh4 = reshape(cell2mat(BH(:,1,:)),[],nsim,2); % all the 4 elements vectors
bh6 = reshape(cell2mat(BH(:,2,:)),[],nsim,2); % all the 6 elements vectors
meanBH4 = squeeze(mean(bh4,2)); % mean over all 4 element vectors
meanBH6 = squeeze(mean(bh6,2)); % mean over all 6 element vectors

但是，以正确的方式做的一步是定义两个数组，一个用于每种方法：

BH1  = zeros(nsim,ddlen,dd(1));
BH2  = zeros(nsim,ddlen,dd(2));

然后在您的循环中为它们分配值：

if method==1
    BH1(ni,di,:) = betahat;
else
    BH2(ni,di,:) = 10*betahat;
end

最后，只需以每个：

的平均值

meanBH1 = mean(BH1,3)
meanBH2 = mean(BH1,3)

---

编辑：

要以一种更" matlabish"的方式写所有这些，我会建议以下内容：

nsim = 10;   % iteration number
dd = [4 6];  % two dimension cases,beta=(beta_1,cdots,beta_d)^T
methods = 2; % two methods
% preapering random seeds
s = bsxfun(@times,1:numel(dd),(1:nsim).');
seednum = cumsum(s(:));
% initialize results array
BH = nan(max(dd),nsim,numel(dd),methods);
counter = 1;
for k = 1:numel(dd)
    for n = 1:nsim
        % set a new random seed from the list:
        rng(seednum(counter));
        % compute all betahats with this seed:
        betahat = randn(max(dd),2).*repmat([1 10],[max(dd) 1]);
        % assign the values to BH by dimension:
        for m = 1:methods
            BH(1:dd(k),n,k,m) = betahat(1:dd(k),m);
        end
        counter = counter+1;
    end
end
% compute the means over iterations:
means = squeeze(mean(BH,2,'omitnan'))

，因此您将获得means作为结果。

---

P.S。我不知道为什么您在每次迭代中都调用randn('seed', seednum)，除了不是推荐的语法，但是如果您可以删除它，那么您可以将大多数循环矢量化，而您的代码挤压到：

% compute all betahats:
betahat = randn(nsim,max(dd),numel(dd),2);
% apply dimensions:
for k = dd
    betahat(:,k+1:end,1,:) = nan;
end
% apply methos 2:
betahat(:,:,:,2) = betahat(:,:,:,2)*10;
% compute the means over iterations:
means = squeeze(mean(betahat,1,'omitnan'))

希望现在看起来更清晰...