demonstrating concept of confidence interval for the mean value?

Can somebody look into my codes. I am trying to design a small simulation experiment to demonstrate concept of confidence interval for the mean.

function [alpha,m] = CIMV(Me,Sigma,n,alpha)

%CIMV= Confidence Interval of Mean Value
%Me=mean, Sigma=standard Deviation, n=record length,
%aplha=confidence interval

X=normrnd(Me,Sigma,[n,100000]);

i=0; k=0; 
Sample=[];

for i=1:100000

      Me_sample(i)=mean(X(:,i),1);
      Sigma_sample(i)=std(X(:,i));
      Se(i)=Sigma_sample(i)/sqrt(n); % Standard Error
      V=tinv(1-0.5*alpha,n-1);
      Me0(i)=Me_sample(i)-V*Se(i);
      Me1(i)=Me_sample(i)+V*Se(i);
      if (Me0(i)<=Me & Me<=Me1(i))
          k=k+1;
      end
      Sample(1,i)=Me_sample(i);
      Sample(2,i)=Sigma_sample(i);
      Sample(3,i)=Se(i);
      Sample(4,i)=Me0(i);
      Sample(5,i)=Me1(i);
      Sample(6,i)=k;
  end
    Sample

end

I run by feeding values as mentioned CIMV(0,1,5,.05) however, instead of 5 sample values i am getting 100000.

Answer 1

[Star Strider]

You are asking it to produce 100000 sample values in your for loop. Change your for loop initial statement to:

for i=1:n

and it should do what you want.

[Hydro]

Exactly.Many thanks

[Star Strider]

My pleasure!

Answer 2

[Matt Tearle]

Can you clarify your intent? It looks like you're doing 100000 experiments where you take 5 (n) random values in each experiment and calculate 6 statistics (mean, std dev, std err, upper bound, lower bound, cumulative sum of how many times the population mean is between the sample bounds). In that case, you should be getting a 6-by-100000 result.

Perhaps you want it the other way around (5 experiments using 100000 values each)? If so, you need to flip a few indices around. But your std err and CI formulas use n (not 100000).

Either way, it would really make your life easier if you vectorized this code. The loop is not necessary, and, as you have it, is actually inefficient.

X = normrnd(Me,Sigma,[n 100000]);
Me_sample = mean(X);
Sigma_sample = std(X);
Se = Sigma_sample/sqrt(size(X,1));
V = tinv(1-0.5*alpha,size(X,1)-1);
Me0 = Me_sample - V*Se;
Me1 = Me_sample + V*Se;
mean_in_int = (Me0 <= Me) & (Me <= Me1);
Sample1 = [Me_sample;Sigma_sample;Se;Me0;Me1;cumsum(mean_in_int)];

You can use n where I've used size(X,1), but the way I have it makes it trivial to swap the roles of n and 100000, if that's what you want. Just change the first line to

X = normrnd(Me,Sigma,[100000 n]);

and everything else should just work.

[As an aside, your function returns alpha and m, but alpha is one of the inputs (which isn't modified) and m is never calculated. I assume you will calculate these from Sample once you have that part working correctly...?]

[Hydro]

Thanks Matt,

Valueable thoughts. Those out put Arguments (Alpha,m) are to Compare the frequency of cases falling within the confidence bounds with the specified confidence level. I still have to do it, however, am struggling a bit. Any hint would be appreciated.