%% Error Propogation Example 1 (See Class Notes)
% Tau=1/(Beta*b)

% Number of repeat measurement to demonstrate the distributions
N=10000;

Beta_mean=3;    %Underlying true value of Beta
sigma_Beta=Beta_mean/10/sqrt(2);  %standard deviation of measured values
Beta=Beta_mean+sigma_Beta*randn(N,1); %Observed values of Beta

b_mean=2;    %Underlying true value of b
sigma_b=b_mean/10/sqrt(2);  %standard deviation of measured values
b=b_mean+sigma_b*randn(N,1); %Observed values of b

Tau=1./Beta./b;

mean(Tau)
std(Tau)

dTau_dBeta=-(1/(Beta_mean^2*b_mean));
dTau_db=-(1/(Beta_mean*b_mean^2));

sigma_Tau=sqrt(sigma_Beta^2*dTau_dBeta^2+sigma_b^2*dTau_db^2)
        
fprintf('\nObserved mean of Tau: %g\n',mean(Tau))
fprintf('Expected mean of Tau: %g\n\n',1/(Beta_mean*b_mean))

fprintf('Observed standard deviation of Tau: %g\n',std(Tau))
fprintf('Expected standard deviation of Tau: %g\n',sigma_Tau)
fprintf('Relative Error in Tau (sigma_Tau/Tau): %g\n\n',sigma_Tau/mean(Tau))


%% Correcting for Bias in Tau

%Correct for second order bias (SOB)
Tau_SOB=1./Beta./b-sigma_Beta.^2*(1./(Beta.^3.*b))...
    - sigma_b.^2*(1./(b.^3.*Beta));

fprintf('\nObserved mean of Bias Corrected Tau: %g\n',mean(Tau_SOB))
fprintf('Expected mean of Tau: %g\n\n',1/(Beta_mean*b_mean))



