How do I implement the standard deviation of the standard deviation in MATLAB?

Are there in MATLAB already predefined functions to calculate the standard deviation of the standard deviation or the standard deviation of the average standard deviation?

(use Matlab R2012a)

Answer 1

[Image Analyst]

You simply use std() again on the whole group (vector or matrix) of standard deviations that you have.

stdOfGroup = std(individualStDevs);

[Image Analyst]

The standard deviation of the standard deviations will give you the correct result. It will give you the standard deviation of the standard deviations. What Star was showing was that this does not give you the standard deviation of the combined populations. But after your explanation I don't think you're after that.

You have a Monte Carlo simulation and it gives some output. Let's call it y. Now if you run the simulations lots of times, y will be different each time. If you plot y versus "number of runs" you might find y to be highly variable (high st dev) for low number of runs but y settles down to a "steady state" as the number of runs increases. So, say for less than a thousand runs, y might be 1000 +/- 200, but for a million runs y might be 1000 +/- 50, and for 10 million runs you might find y = 1000 +/- 10. So stdev decreases from 200 to 50 to 10 as the number of experiments/runs/trials increases from a thousand to ten million.

So what you might want to do is to plot y vs. number of runs and plot std(y) vs. # runs. You might have a funnel shaped scatterplot and you will probably observe how the width of the scatter decreases with number of runs. The y funnel will settle down to the steady state mean of y while the stdev funnel will settle down to 0 eventually (hopefully but not necessarily). So how do you decide on the number of runs to use? You could do a moving mean with moveman() or conv() and see where the mean of this scatter plot levels out. Or you could use movstd() or stdfilt() to get a measure of the "width" of the funnel by using it's standard deviation. So if your stddev funnel starts at 50 and narrows down towards zero, you might somewhat arbitrarily say a std dev of std of 5 (a funnel envelope width of 5) is good for you, then you just see how many runs that corresponds to.

[Glazio]

@Image Analyst: Thank you, that's what I meant :-)

[Image Analyst]

I think it's worth looking up for the sampsizepwr() function Star mentioned. That's something I didn't know existed, and could be helpful to you and many other people.

[Glazio]

@Image Analyst: movstd() and movemean() don't work in MATLAB 2012 a.

Does filter() represent a useful alternative since conv() and stdfilt() provide unsuitable results?

[Image Analyst]

They could, but you state an incorrect statement, that conv() and stdfilt() don't work. They can work to provide a moving mean or std. Here's proof:

vec = rand(1, 100)
kernel = ones(1,3)/3;
f = conv(vec, kernel, 'same')
nhood = ones(1,3);
fsd = stdfilt(vec, nhood)

[Glazio]

@Image Analyst: Thanks for your answer. stdfilt(), conv() and filter() give the following results (for a test example):

stdfilt:
   1.0e-07 *

      0.8429
           0
           0
           0

This result looks strange.

conv:

      2.4814
      2.8835
      1.2061
      0.8041

The result of filter is of the same order of magnitude as conv.

filter:
    6.2382
    1.2061
    1.2061
    1.2061

[Glazio]

@Image Analyst: What can be the problem if both graphs always fluctuate by the same value and no trend can be detected?

[Image Analyst]

Glazio, I don't understand the question. You describe signals that are highly correlated. I don't know what else to say. Why is that a problem? Two graphs, one graph of y and one graph of y+10, always fluctuate by the same value as x moves along. That is not inherently a problem. Why do you say it is a problem?

[Glazio]

@Image Analyst: Sorry for this unclear question.

If I plot y vs. number of runs (in this example 1000 runs) than the width of the scatter doesn't decrease with number of runs and so the graph doesn't settle down to a steady state. There is always a fluctuation and no trend.

The same applies to the standard deviation: if I plot std(y) vs. number of runs, the stdev funnel doesn't settle down to 0. There is no trend within this graph, but only small fluctuations around a value.

Can this result be considered plausible or is there a data-internal problem (or a bug in the code)?

[Image Analyst]

That's completely reasonable for a "stationary stochastic process". See https://en.wikipedia.org/wiki/Stationary_process

For example if I took a snapshot of something every second for an hour, the mean intensity (your "y") might vary by +/- 5 gray levels. And it might always be like that no matter how long I take pictures as long as the scene or other image capture parameters don't change. Perfectly plausible. The standard deviation for a given image might also change. It could be 5 for one photo, 4.5 for another, and 5.6 for another, averaging out to 5 over hundreds of snapshots. There is no reason to think the variation would decrease over time unless there is some reason for it to do so. Why do you think your fluctuations should decrease over time?

[Glazio]

@Image Analyst: Thanks for your answer.

If I am not wrong, Monte Carlo optimization requires a minimum number of runs to achieve a good coverage of my solution space. Now it should be possible to assess visually the standard dev. + std of std with the number of runs. From this representation the minimum required number of runs can be identified.

In a similar approach, the settling down could be observed. If this does not occur, the minimum number of runs can not be visually identified!?

[Image Analyst]

Well that's a different situation. In that case, your y is a function of the number of experiments, not a value that you get from just one experiment. So i f y is the mean, then the "standard error of the mean" with a run of a million experiments will be smaller than a run with only 5 experiments. See https://en.wikipedia.org/wiki/Standard_error#Standard_error_of_the_mean

But if you had say a 6-sided die, and looked at the standard deviation of the face values after 100 rolls, then if you did that twice, or did it a thousand times, the standard deviation will be around 1.7 each time.

Basically I think you're getting confused between standard error of the mean and standard deviation, which are different concepts.

[Glazio]

@Image Analyst: Thank you, yes I think I have exchanged the terms with each other.

If I now plot the standard error of the mean against the number of runs, the graph settles down to 0 (after 500 runs).

If I now plot the mean value against the number of runs, the graph does not settle down to a steady state. There are always slight fluctuations.

[Image Analyst]

Yep, sounds right.