function [p,stddev]=PvalMean(samp,nullmean,N)
if nargin<3 N=20000; end
% samp is a vector containing numerical data.
% The test statistic here is xbar, the mean of samp, and that mean is to be compared to a number nullmean. The p-values calculated here estimate the probabilities of sample mean differences from nullmean greater than or equal to those actually observed assuming the null hypothesis which is the statement that the true population mean is nullmean.
% If the probability is very small it is unlikely that the null hypothesis that generated it is true, and this counts as evidence in favor of the alternative hypothesis.
% Let diff1=abs(xbar-nullmean).
% p(1) estimates the probability that you would get a sample whose mean is to the left of nullmean-diff1 provided nullmean is the true population mean.
% p(2) estimates the probability that you would get a sample whose mean is to the right of nullmean+diff1 provided nullmean is the true population mean.
% p(3) is the sum of these two, and corresponds to the 2-sided alternative hypothesis.
% Heed the warning provided as part of the function outout: if you have a 1-sided hypothesis and mean(samp) is ``on the wrong side'' these probabilities are not relevant to your test: the data in samp automatically constitutes evidence against your alternative hypothesis, no test required!
% Here is a list of numbers to test.
% samp=[72,65,81,73,53,55,78,43, 45,12,33,24,49,87,10,69];
n=length(samp);
xbar=mean(samp);
diff1=abs(xbar-nullmean);
rng('shuffle') %assign time-based seed to random number generator
counttwosided=1;
countlesser=0;
countgreater=0;
if nullmean<=xbar countlesser=1;
'The sample mean is to the right of the null hypothesis mean. If you think it should be to the left (a 1-sided alternative hypothesis) this data automatically constitutes evidence against your alternative hypothesis, no calculation required!'
else countgreater=1;
'The sample mean is to the left of the null hypothesis mean. If you think it should be to the right (a 1-sided alternative hypothesis) this data automatically constitutes evidence against your alternative hypothesis, no calculation required!'
end
for i=1:N
y(i) = mean(datasample(samp,n));
% The histogram of y should have (nearly) the SAME SHAPE AS THE SAMPLING DISTRIBUTION OF xbar, but centered at xbar instead of the true mean.
w(i) = y(i)-xbar+nullmean;
% The histogram of w should have (nearly) the SAME SHAPE AS THE SAMPLING DISTRIBUTION OF xbar, but centered at the null hypothesis mean nullmean instead of the unknown population mean. This is the distribution we use to calculate the P-values.
% We count the fraction of w-values to the right of nullmean+diff1 and the fraction to the left of nullmean-diff1.
if w(i)>=nullmean+diff1
countgreater=countgreater+1;
end
if w(i)<=nullmean-diff1
countlesser=countlesser+1;
end
end
ToTheLeft=countlesser/(N+1)
ToTheRight=countgreater/(N+1)
p=[ToTheLeft,ToTheRight,ToTheLeft+ToTheRight];
% (1 is added because the original sample automatically satisfies one of the conditions.)
figure(74)
histogram(w,100);
hold on
plot([nullmean-diff1 nullmean-diff1 nullmean+diff1 nullmean+diff1],[0 N/40 N/40 0])
plot([nullmean nullmean],[0 N/30])
hold off
stddev=std(w);
end