%In this script, we will fit the data extracted from the 1958 paper on
%myoglobin O2 binding. This data was never digital, so we had to play some
%tricks to get it. We extracted the data from the plot using the software
%DigitizeIt. We saved the points as a csv file (comma separated values)
%which is composed of two columns which correspond to x and y coordinates
%relative to the axis of the plot. We can now extract the data from the csv
%file in MATLAB.
data = csvread('Rossi-Fanelli1958fig2.csv');
%The xy coordinates are now in a 15x2 matrix. Before we plot the data, we'll
%need to separate the x and y coordinates. Matlab indexes matrices in the
%order (row, column). For example, if I want to pull the y coordinate of
%the third point, I would do the following.
data(3,2);
%We want to pull all values out, however. We can do this as follows.
dataX = data(:, 1);
dataY = data(:, 2);
%Now we can plot the data!
figure(1)
plot(dataX, dataY, 'ko'); %Will show as dots with radius=5.
xlabel('PO2 (mmHg)');
ylabel('P bound');
%Let's make a street fighting theory curve. By eye, we can estimate that
%the Kd should be around 1mmHG O2. To generate the curve, we'll first make
%a range of values from 0 to 8 mmHg linearly spaced.
xRange = linspace(0, 8); %100 linearly spaced points between 0 and 8.
Kd = 0.6; %Guess for Kd in mmHg. Try changing this!
Pstreetfighting = (xRange/Kd) ./ (1 + xRange/Kd); %Using our equation.
%Now we can plot our estimate.
hold on
plot(xRange, Pstreetfighting, 'r-');
legend('Points', 'Prediction');
hold off
%Now let's do the actual curve fit. Matlab has some built in functions for
%curvefitting, but it is more instructive to do this by hand. If you have
%the curve fitting package installed, you could use the 'fittype' command
%and specify the equation you would like MATLAB to fit. Some students do
%not have the the curve fitting package, so we'll have to do this another
%way by writing an external function. We'll make our fit by minimize the
%squares of the distances from the point to the line. The best fitting
%curve will have the smallest distance from the points. See Hernan's notes
%on the course website for a more thorough description.
%As an introduction to functions, let's make a new function called
%'BindingChiSq'. This will be in a separate .m file. Once the function is
%written, we can call it from any script! Let's plot the chi sq. as a
%function of Kd. We need to run a for loop because our function only takes
%a specific value of Kd, not an array.
figure(2);
KdRange = linspace(0, 8);
for i=1:length(KdRange)
chisq(i) = BindingChiSq(dataX, dataY, KdRange(i));
end
plot(KdRange, chisq, '.-k', 'Markersize', 10); %Marker size of 10 px.
xlabel('Kd (mmHg)');
ylabel('chi^2');
%Now we have an idea of how Chi sq helps us optimize parameters to generate
%a good data fit. We'll define a new function that will do the same
%procedure, except it will incrementally change the Kd to find the minimum
%Chi sq value.
StartingParameters = 0.6; %Our best guess for Kd.
FitResults = lsqnonlin(@(KdFit) BindingChiSq(dataX, dataY, KdFit),...
StartingParameters, [], [])
%Now let's plot the theory curve and our data together to see the results
%of our true fit.
figure(3)
plot(dataX, dataY, 'ok');
hold on
plot(xRange, (xRange/Kd) ./ (1 + xRange/Kd), 'b-');
plot(xRange, (xRange/FitResults) ./ (1 + xRange/FitResults), '-r');
xlabel('PO2 (mmHg)');
ylabel('Pbound');
legend('Data', 'Prediction', 'Theory', 'Location', 'Northeast');
hold off