Updating mean and standard deviation

2 次查看（过去 30 天）
Michael Madelaire
评论： dpb 2018-4-10
Hi
I am working on a project which includes Bayesian inversion. And am working on an idea about updating the prior model iteratively, in a loop.
Here is a simplification of the code that troubles me.
xAll = zeros(4, 400000);
for i = 1:400000
x = randi(5, 4, 1);
xAll(:, i) = x;
mu = mean(xAll, 2);
sigma = std(xAll, 0, 2);
fvec = (mu-x)./sigma;
logprior = -1/2*sum(fvec.^2);
% Then comes lots of other stuff that uses the above, but is not relevant for the time being.
% But normally this section would determine a new "x".
end
When determining the "logprior"-value I need the mean and std of "xAll". This takes time to do over and over again. Do any of you know a way of updating "mu" and "sigma", without recomputing them?
Thanks in advance

请先登录，再回答此问题。

采纳的回答

dpb
dpb 2018-4-9
编辑：dpb 2018-4-9
You can keep the partial sums
mu=SUM(x)/N --> SUMX/N
SUMX=SUMX+xnew-xold;
xold=xnew;
in each iteration. Similarly for STD excepting need SUMXSQ as well;
v = [sum(x.^2)-1/N*(sum(x))^2]/(N-1)
sd=sqrt(v)
or, simplifying the sums to variables;
v=[SUMXSQ-N*MU^2]/(N-1)
IOW, you keep the two running sums of sum(x) and sum(x^2) and update them for each iteration by removing the previous and adding the new terms.
From what's shown it doesn't appear that there's any reason you couldn't vectorize the i loop away, but perhaps what's not shown precludes doing so...
  6 个评论
显示 4更早的评论
Michael Madelaire
Michael Madelaire 2018-4-10
Thanks for the clarification.
I will accept the answer since it is a solution. Although it is not specifically what I am looking for.
dpb
dpb 2018-4-10
What, then, specifically, are you looking for? The only way I know to compute mean/var on the fly is as shown; keep the intermediary terms and update. If one talks of exceedingly long sequences it becomes rather impractical, unfortunately. You could make approximations; but that's different and if the difference in computed value once you've gotten several K of observations is significant enough to bother to recompute then the approximation probably isn't good enough, either, as compared to just using the previously calculated value.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Parametric Spectral Estimation 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by