%BetaDemo
% This program calls betaval.m, a function to return a beta value
% from a Beta distribution with Mean = mn and standard
% deviation = sd. The other input is fx, a value for the
% cumulative distribution function. To pick a random value, make
% fx = rand. 
% This program calculates random values for each of several mean
% and st. deviation values

clear all;
%*****************Simulation Parameters***********************
means = [0.5,0.9] ; %  means to use;
psds = [0.1,0.5,0.90]; %percentage of max possible s.d. to use
reps = 1000; % how many values to make
%*************************************************************
results = [];
rand('state',sum(100*clock)); % seeds random number generator
meannum = length(means); % how many means
psdnum = length(psds);   % how many s.d. values

% now, 2 loops to get all the combinations of means and s.d.s
count = 1
for meani =1:meannum
 for psdi = 1:psdnum 
	mn = means(meani);
	sd = ((mn*(1-mn))^0.5)*psds(psdi); % calculate sd
	for fxi = 1:reps
		temp(fxi) = betaval(mn,sd,rand); % find random values
	end; % fxi
    results(:,count) = [mn; sd; temp']; % store the values
    count = count+1
 end; % psdi
end; % meani

% make a set of frequencies and bin means for the results
[freqs,bins] = hist(results(3:(reps+2),:),40);  
disp('In the following histogram results')
disp('the first two rows of each column are')
disp('the mean and standard deviation, and')
disp('the first column has the bin means')
finalresults = [zeros(2,1),results(1:2,:); bins,freqs]
hist(results(3:(reps+2),:),40) % make a histogram of all results