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]};