% Free Fall Synthetic Data
% Bayesian Data Analysis
% Prof. Kevin H. Knuth
% 18 and 19 Oct 2018
% Set g = 9.8 (m/s^2)
% Set y = 2 m
g = 9.8; % m/s^2
y = 2; % m
s = 0.2; % s (reaction time)
N = 1000; % N=1000 trials
% compute the time to fall to the ground from y=2m up
% Plus a Gaussian distributed reaction time, which is kind of a hack
% since this is not the best model of reaction time.
t = sqrt(2*y/g) + 0.2 * randn(1,1000);
% In this problem, with a uniform prior probability for g,
% the posterior probability is proportional to the likelihood, which
% has been selected to be Gaussian.
%
% This loop tries many values of g ranging from 9 m/s^2 to 11 m/s^2.
% The (unnormalized) posterior probability is computed in p(i)
% The logarithm of the (unnormalized) posterior probability is computed in logp(i)
i = 1;
for gg = 9:0.01:11
g(i) = gg;
thesum = sum(t - sqrt(2*y/g(i)));
p(i) = (2*pi*s^2)^-(N/2) * exp( -1/(2*s^2) * thesum^2 );
logp(i) = log(p(i));
i = i + 1;
end
% Plot the posterior probability
figure;
plot(g,p)
title('Posterior Probability');
xlabel('Acceleration of Gravity (g in m/s^2)');
ylabel('Unnormalized Probability');
% Plot the log posterior probability
figure;
plot(g,logp)
title('Log Posterior Probability');
xlabel('Acceleration of Gravity (g in m/s^2)');
ylabel('Unnormalized Log Probability');
% Note that the peak of the Posterior Probability is the same as the peak
% of the Log Posterior Probability. However, the Log Posterior is a much
% broader and smoother function, which makes numerical peak-finding easier.
% Also the values of the Posterior are on the order of 10^300 (1e+300),
% which is getting close to the largest number that Matlab can handle:
% realmax = 1.8e+308.
% Whereas the Log Posterior's numeric range is more manageable.
%
% For these reasons, we oftem work with the Log Posterior, rather than the
% posterior.