University of South Florida

We’ll fix *your* carbon...

Matlab program for calculating kinetic isotope effects using Pitman estimators

 

This is the program cited in

 

Scott K. M., Fox G., and Girguis P. R. (2010).  Measuring isotope fractionation by autotrophic microorganisms and enzymes.  In “Methods in Methane Metabolism”, Methods in Enzymology, ed. S. Ragsdale.  Elsevier Inc., Cambridge. In press.

 

Scott, K.M., Lu, X., Cavanaugh, C.M., and Liu, J. (2004) Optimal methods for estimating kinetic isotope effects from different forms of the Rayleigh distillation equation. Geochimica et Cosmochimica Acta 68:  433-442.

 

 

 

 

 

Related Links

function [alpha, CI] = PitmanAlpha( C, R, n, range)

 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Pitman estimator for estimating alpha from multiple datasets

%

% Usage: [alpha, CI] = PitmanAlpha( C, R, n, range )

%

% Input Arguments:

%   C:   Reactant concentrations

%   R:   Isotope ratios

%        C and R should be cell arrays, with each cell being a vector

%        of the data from one experiment.

%        The number of vectors of C and R should be the same,

%        the number of measurement points of C and R in corresponding

%        vectors should be the same, too.

%        Different number of measurement points in different runs is allowed.

%   n: number of discretization steps in the digital integral

%        optional input argument, larger n values yield more precise results.

%      default: 4000

%   range: range for digital integral, optional input argument,.

%        This range should be sufficiently large, otherwise the digital integral

%        will not generate precise results.

%      default: center on the average of the slopes,

%        then choose a width that is 10 times the average standard deviation

%        or 5 times the range of the three slope values, which ever is greater.

%

% Output Arguments:

%   alpha: the alpha value to be estimated

%   CI: 95% confidence interval

%

% Example: Three experiments, one has 7 datapoints, one has 6, and one has 5.

%   C = {[4.3856    3.8050    3.1714    2.6484    2.1228    1.4681    0.8586],

%        [4.5339    3.7489    3.2564    2.6904    2.0016    1.4690],

%        [4.3238    3.8756    3.2132    2.5901    2.1444]};

%   R = {[0.0110    0.0111    0.0111    0.0112    0.0112    0.0113    0.0115],

%        [0.0110    0.0111    0.0111    0.0112    0.0112    0.0113],

%        [0.0111    0.0113    0.0112    0.0113    0.0112]};

%   [alpha, CI] = PitmanAlpha(C, R)

%

%   alpha =

%       1.0240

%   CI =

%       1.0191    1.0282

%  

% References:

%   Scott K. M., Lu X., Cavanaugh C. M., and Liu J. (2003)

%       Optimal methods for estimating kinetic isotope effects

%       from different forms of the Rayleigh distillation equation.

%       Geochimica Cosmochimica Acta (submitted).

%   Pitman E.J.G. (1939)

%       The estimation of the location and scale parameters of a continuous

%       population of any given form. biometrika, 30(3/4): 391-421

%

% Author:

%   Xin Lu. Dept. Statistics, Harvard University

%       Science Center 802, 1 Oxford St. Cambridge, MA, 02138

%      

%   Nov. 5, 2002

%

% Instruction to user:

%      Please save this file as a text file with “.m” as the extension (not “.txt”),

%      for example “PitmanAlpha.m”, under the current working directory

%      or Matlab search path.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 

%check number of input arguments, if less than 2, then print help information

if( nargin < 2 ) 

    help PitmanAlpha;

    return;

end

 

% check the number of runs

nRun1 = length(C); nRun2 = length(R);

if( nRun1 ~= nRun2 )

    error('Error input: number of runs in C and R is not equal');

end

nRun = nRun1;

if( nRun < 2 )

    error('Error input: C and R should contain at least 2 experiments');

end

 

if( nargin < 3 )    % set default n, if not inputed

    n = 4000;

end

 

for i=1:nRun

    CC = C{i}; RR = R{i};

    % check the number of measurement points in each run

    nP1 = length( CC );  nP2 = length( RR );

    if( nP1 ~= nP2 )

        s = sprintf('Number of measurement point is not equal');

        error(s);

    end

    CC = -1 * log(CC); RR = log(RR); nPoint = nP1;

 

    b(i) = (nPoint*sum(CC.*RR)-sum(CC)*sum(RR))/(nPoint*sum(CC.^2)-(sum(CC))^2);

    a(i) = sum( RR - b(i)*CC ) / nPoint;

    SumSqureError = sum( (RR - (a(i)+b(i)*CC)).^2);

    Db(i) = sqrt( SumSqureError / (nPoint-2) / sum( (CC-mean(CC)).^2) );

   

    if( Db(i) < 1e-10 )

        s = sprintf('The std of No.%d run is too small, data error?', i );

        warning(s);

    end

end

 

if( nargin < 4 )    % the range for digital integral

    w1 = 10*mean(Db); w2 = 5 * max( abs(diff(b))); w = max(w1,w2);

    m1 = mean(b)-w/2; m2 = mean(b)+w/2;

else

    m1 = range(1); m2 = range(2);

end

 

% the posterior distribution of true slope:

% product of the multiple t-distributions with degree of freedom v

step = (m2-m1)/(n-1); x = m1:step:m2; ft = ones(1,n); v = nPoint-2;

for i=1:nRun

    fv = tpdf( (b(i)-x)/Db(i), v)/Db(i);    ft = ft .* fv;

end

intft = sum(ft)*step;   ft = ft / intft;   

 

% integral the posterior distribution, and the 95%CI

cumft = cumsum(ft) * step;

CI(1) = x( min(find( cumft > 0.025) ) );   

CI(2) = x( max(find( cumft < 0.975) ) );

CI = 1./(1-CI);

 

% posterior mean (Expectation) of b

beta = x*ft'*step;

alpha = 1/(1-beta);

 

C = {[5.740574417  5.229865252  5.100165892  5.095393011  4.827003719  4.456000006       4.243410949],

[4.817738562       4.364664072  4.205424502  4.070863002  3.739159867  3.526208647       3.450679292],

[4.937908884       4.559116406  4.407069479  3.952054811  3.778851276  3.74923823   3.573326816       3.662250814]};

R = {[0.011215088  0.011233179  0.011262283  0.011298579  0.011286668  0.011323413       0.011334762],

[0.011210256       0.011246192  0.011268014  0.011288016  0.011292286  0.011337684       0.011355438],

[0.011219583 0.011235629 0.011261047 0.011296893 0.011324649 0.011332177 0.011339257 0.011340718]};