Input argument error length of variables

I am trying to use the following script to process some of my data. I keep getting this error:

Error using nargin
You can only call nargin/nargout from within a MATLAB function.

   function [corrected, details] = pin(year,d13c,alpha,varargin)
   ↑
  Error: Function definitions are not permitted in this context.

] This code was written in 2009 so I don't know if it is outdated or what is going on. I need this to work so any help would be greatly appriciated!

[Massimo Zanetti]

You should upload your code otherwise how can we solve your problem? We cannot guess.......

[Geoff Hayes]

Elizabeth - the second error message suggests that you have some code that precedes the function definition/signature. If you are defining a function, then it must be the first line in your file.

Answer 1

[Marc Jakobi]

The problem is exactly what the error message says:

Somewhere in your script is the keyword "nargin". nargin is used to check the amount of inputs of a function and because of that it can only be used within a function, not within a script.

Example:

function numinputs = checkinputs(a, b, c)
   if nargin == 1
      numinputs = 1; %a
   elseif nargin == 2
      numinputs = 2; %a and b
   else
      numinputs = 3; %a, b and c
   end
end

[Elizabeth Olson]

So do I delete it from the script? Or do I run the function some other way?

[Elizabeth Olson]

This is the code: Appendix A: Supplementary data. Matlab Code to Perform the ‘pin’ Correction

function [corrected, details] = pin(year,d13c,alpha,varargin)

% PIN - PreINdustrial Correction for Carbon Isotope Tree-Ring Series for pC02
% 
%     [corrected, [details]] = pin(year,c,span)
%
%  A user-written Matlab function that performs a correction on carbon
%  isotope series from tree-rings to account for changes in atmospheric
%  pCO2 during the last ~150 years.
% 
%  The function is based on the physiological mechanisms and logical
%  constraints discussed in detail in McCarroll et al. [2009].  
%
%
% Inputs:   year = vector (m x 1) of years (required)
%           d13c    = vector (m x n) of atmospheric-corrected d13C values
%                  [there can be missing values (designated as 'NaN')]
%                  (required)  
%           alpha = [optional, default = 0.65]. This is the fraction of the
%                   time series used in the local (LOESS) regression (see
%                   Cleveland 1993). Must be between 0 and 1.
%
%           'verbose' = will print diagnostic figures and save results for debugging
%
%           'noplot' = do not create the final plot (batch mode)
%
% Outputs:  corrected = vector of PIN corrected d13C values
%           details   = optional structure for debugging and supporting information
%

% Revision History
% v0.1 - Original coding of loess regression
% v0.2 - full code for pin correction from initial d13C data
% v0.3 - heavy debugging to match R/Excel output and diagnose differences
% v0.4 - modified to accept multiple carbon isotope series
% v0.5 - minor corrections, extend ca data to 2005
% v0.9 - peer-review version 
% v0.91 - changes based on reviews and new testing
% v0.92 - heavy changes based on peer review
% v0.93 - heavy debugging to match R package
% v1.00 - final version for publication
%
% References: 
%
% McCarroll, D., Gagen, M. H., Loader, N. J., Robertson, I., Anchukaitis, K. J.,
% Los, S., Young, G. H. F., Jalkanen, R., Kirchhefer, A., and Waterhouse, J. S.
% (2009). Correction of tree ring stable carbon isotope chronologies for changes
% in the carbon dioxide content of the atmosphere. Geochimica et Cosmochimica 
% Acta, 73(6):1539?1547.
%
% Cleveland, W.S. and S. J. Devlin, 1988, Locally-Weighted Regression:
% An Approach to Regression Analysis by Local Fitting. Journal of the 
% American Statistical Association, 83:596-610.
%
% Cleveland, W.S. and E. Grosse, 1991, Computational Methods for Local
% Regression. Statistics and Computing, 1:47-62, 1991
%
% Cleveland, W.S. 1993, Visualizing Data, Hobart Press.

% Original Code Author: Kevin J Anchukaitis ([e-mail])
% Loess smoothing based on data visualization toolbox code by Cleveland 1993

%%Screen display
disp([' '])
disp(['PreINdustrial Correction for d13C tree-ring time series [PIN]'])
disp(['-------------------------------------------------------------'])
disp(['Version 1.00  March 2009   McCarroll et al. 2009'])
disp([' '])
disp(['McCarroll et al. 2009, Correction of tree ring stable carbon'])
disp(['isotope chronologies for changes in the carbon dioxide'])
disp(['content of the atmosphere, Geochim. Cosmochim. Acta, 73(6):1539?1547'])
disp([' '])

%%Check for verbose mode
if nargin > 3;
 if strcmp(varargin{1}, 'verbose')
  verboseMode = 1; % Extra graphics and output for debugging
 else
  verboseMode = 0; % No verbose debugging
 end
end

%%Check for plotting mode
if nargin > 3;
 if strcmp(varargin{1}, 'noplot')
  plotMode = 0; % Final graphic file
 else
  plotMode = 1; % no final graphic file
 end
end

%%Error Checking 
% Is the number of years equal to the number of isotope records?
if length(d13c) ~= length(year)
    error('Length of year vector and length of carbon isotope data must be the same')
end

% get the span value, or set to the default span 
% some intermingling here of of 'span' and 'alpha' is for consistency with R
if nargin < 3; alpha = 0.65; else alpha = alpha; end % a relatively wide 'local' span (alpha)
if any(size(alpha))~=1 | alpha < 0 | alpha > 1; 
  error('Alpha should be a single number between 0 and 1')
end

%%STEP 1. Begin Correction Procedure
% Sort seamlessly such that the data is by ascending year and arrange data vectors
data = sortrows([year d13c],1); year = data(:,1); d13c = data(:,2:end);
disp('Input data are automatically reordered prior to correction if necessary')

% save the original data before we start to operate on it
original.d13c.year = year; original.d13c.ratio = d13c;
d13co = d13c;

% How many individual series are there (nseries)?
[mseries,nseries] = size(d13c);

% Preallocate the output structure for speed and error catching
corrected.d13c.year = year;  corrected.d13c.ratio = NaN * ones(mseries,nseries);

% loop over for the number of individual series
for iter = 1:nseries

   graphname = ['pinDiagnostics',num2str(iter),'.ps'];
   disp([' '])
   disp(['Correcting carbon isotope series ', num2str(iter,0), ' ... '])
   full = ~isnan(d13co(:,iter)); % find the continous series
   yyear = year(full); % this is span of years for the series we'll work on
   cc = d13co(full,iter); % this is the series we'll work on

%%STEP 1 Retrieve Annual Values of ca 
if min(yyear) >= 1850  % do not operate on preindustrial values, result sensitive to this date
 firstyear = min(yyear);
else
 firstyear = 1850;
end ; 

if max(yyear) > 2005;
 lastyear = 2005;
else
 lastyear = max(yyear);
end

ca = co2(firstyear,lastyear); 
ic = find(yyear >= firstyear & yyear <= lastyear);
iyear = yyear(ic);
c = cc(ic);  % c and iyear go together now

% STEP 2 convert d13c to smoothed d13c using loess
[m,n] = size(iyear);

   % ... using loess regression (Cleveland and Grosse 1991)
   % Create the loess weights
   span = alpha*(lastyear-firstyear)/(2002-1850); % modify alpha but should be close to input 
   q = floor(span*m);  % alpha expressed as percentage of series length
   q = max(q,1); q = min(q,m);

   datails(iter).alpha = alpha;
   details(iter).span = span; % store the span size
   details(iter).spanN = q; % store the sample size

   for i = 1:m ; % for every year ...
     delta = abs(iyear(i) - iyear(:)); % get the distance from current point ...
     dsorted = sort(delta); % ... sort the values from lowest to highest (closest to furthest) ...
     dmax = dsorted(q); % ... find the maximum value within the desired span ...
     w = (1-abs(delta./(dmax)).^3).^3; % ... and calculate the tricubic weights.

    % index into the span of interest only
    gtz = find(w > 0); % find only the weights > 0, which will limit the least squares regression to the alpha proportion    
    xc = iyear(gtz); yc = c(gtz); wc = w(gtz); % cull the span from the data
    WC = spdiags(wc,0,length(wc),length(wc));

    % Construct the Vandermonde matrix for a second order (degree 2) polynomial
    V = [xc.^2 xc ones(length(xc),1)]; % simpler than built in function vender() for this operation 

    % QR decomposition for polynomial coefficients
    warning off;
    [Q,R] = qr((V'*WC*V),0); % qr is a function in both Matlab and Octave
    p = R\(Q'*(V'*WC*yc));  % this is most accurate result

    % use polyval.m (in both Matlab and Octave) to get predicted ci values from ca
    yhat(i) = polyval(p,iyear(i));
    cLow(i) = yhat(i);
    pp(i,:) = p;
   end

d13c.low = cLow(:); 
d13c.high = c - cLow(:);

warning on

if verboseMode == 1
    figure(1); clf;
    plot(iyear,c,'k'); ylabel('d13c'); xlabel('year'); hold on;
    plot(iyear,d13c.low,'b','linewidth',1.5);
    title('high and low frequency d13C')
    legend('Observed d13C','Low Frequency d13C','location','best')
    legend boxoff
    print('-dpsc',graphname)
end

clear m n p

%%STEP 3 Convert corrected low frequency d13C into ci
ci.actual = ca.ratio .* ((-10.8 - d13c.low)/22.6);

%%STEP 4 calculate ci for constant ci/ca ratio
ca.initial = ca.ratio(1); % set the initial ca
ci.initial = ci.actual(1); % set the initial ci to the actual ci
d13c.initial = d13c.low(1); 
ci.rcica.constant = ca.ratio * (ci.initial/ca.initial); % ci expected for constant ratio (contraint 1)

%%STEP 5 calculate ci for constant difference between ci - ca 
ci.dcica.constant = ca.ratio + (ci.initial - ca.initial); % ci expected for constant difference (constrait 2)

%%STEP 6 convert ci for constant ci/ca ratio and constant ci-ca difference to d13C units
d13c.rcica.ts = -10.8 - (22.6.*(ci.rcica.constant./ca.ratio)); % this should to be (nearly) constant
d13c.dcica.ts = -10.8 - (22.6.*(ci.dcica.constant./ca.ratio));

%%STEP 7 plot low frequency d13 c and constraints
if verboseMode
 figure(2); clf;
 plot(iyear,d13c.low,'k','linewidth',1.5); hold on; % ylim=range(c(d13c.low, d13c.rcica, d13c.dcica)) , type="l" , xlab="Year", ylab="d13C")
 plot(iyear,d13c.rcica.ts,'y','linewidth',1.5)
 plot(iyear,d13c.dcica.ts,'b','linewidth',1.5)
 legend('low frequency d13C','Constant Ratio Constraint','Constant Difference Contraint','location','best')
 legend boxoff
 grid on
 title('Comparison of low frequency d13c and constraints')
 print('-dpsc','-append',graphname)
end

%%STEP 8 Calculate increments for low frequency d13c (d13c.low) and for its 2 constraints
inc.d13c.low = [0; diff(d13c.low)];
inc.d13c.rcica = [0; diff(d13c.rcica.ts)];
inc.d13c.dcica = [0; diff(d13c.dcica.ts)];

% Preallocated final increment and set to zero
fin.inc1 = zeros(size(inc.d13c.low)); 
fin.inc2 = zeros(size(inc.d13c.low));

%%STEP 9 compare d13c increments against logical contraints; retain d13c increments only if within contraints
fin.inc1(find(inc.d13c.low > 0)) = inc.d13c.low(find(inc.d13c.low > 0)); % Constraint 1
fin.inc2(find((inc.d13c.low - inc.d13c.dcica) < 0)) = inc.d13c.low(find((inc.d13c.low - inc.d13c.dcica) < 0)) - inc.d13c.dcica(find((inc.d13c.low - inc.d13c.dcica) < 0)); % Constraint 2

%%STEP 10 add d13c increments after contraints are applied
fin.inc = fin.inc1 + fin.inc2;

if verboseMode
    figure(3); clf
    plot(iyear,inc.d13c.low,'k--','linewidth',2); hold on;
    plot(iyear,inc.d13c.rcica,'y--','linewidth',2);
    plot(iyear,inc.d13c.dcica,'b--','linewidth',2);
end

%%STEP 11 calculate cumulative correction and add to low frequency d13c; implement no change prior to 1850
fin.inc(find(iyear<1850)) = 0;

if verboseMode
  figure(3); 
  plot(iyear,fin.inc,'color',[0.75 0.75 0.75]);
  grid on
  legend('increment d13C low','increment d13c constant ratio constraint','increment d13c constant difference restraint','increment final','location','best');
  legend boxoff
  title('final, total, and cumulative correction')
 print('-dpsc','-append',graphname)
end

corrected.d13c.low = d13c.initial + cumsum(fin.inc);
details(iter).cumsum = cumsum(fin.inc);

% corrected.d13c.low(find(iyear<1850)) = d13c.low(find(iyear < 1850))

%%STEP 12 add high frequency back
[cyear,ai,bi] = intersect(corrected.d13c.year,iyear);
corrected.d13c.ratio(:,iter) = d13co(:,iter); % the full thing
corrected.d13c.ratio(ai,iter) = corrected.d13c.low + d13c.high; % now overwrite with corrected values

if verboseMode
 figure(4);
 plot (iyear,-c+corrected.d13c.ratio(ai,iter),'k','linewidth',1)
 legend('difference between original and corrected d13C series','location','best')
 title('Magnitude of PIN correction')
 print('-dpsc','-append',graphname)
end  

if plotMode & verboseMode
    figure(5); clf;
    plot(iyear,c,'k'); hold on
    plot(iyear,corrected.d13c.low,'b','linewidth',2)
    plot(iyear,d13c.low,'k','linewidth',2)
    plot(corrected.d13c.year(ai,1), corrected.d13c.ratio(ai,iter),'b');
    xlim([min(corrected.d13c.year(ai,1)) max(corrected.d13c.year(ai,1))]);
    legend('Measured','Corrected','location','best')
    legend boxoff
else
    figure; clf;
    plot(iyear,c,'k'); hold on
    plot(iyear,corrected.d13c.low,'b','linewidth',2)
    plot(iyear,d13c.low,'k','linewidth',2)
    plot(corrected.d13c.year(ai,1), corrected.d13c.ratio(ai,iter),'b');
    legend('Measured','Corrected','location','best')
    legend boxoff
end
 print('-dpsc','-append',graphname)

% add the first, uncorrected part back to the time series
if (firstyear < min(year))
 corrected.d13c.ratio = [cc(1:max(ic)-1)' corrected.d13c.ratio];
end

if verboseMode
    figure(6); clf
    plot(corrected.d13c.year,corrected.d13c.ratio(:,iter),'b'); hold on;
    plot(original.d13c.year,original.d13c.ratio(:,iter),'k');
    xlabel('year'); ylabel('d13C');
    xlim([min(corrected.d13c.year) max(corrected.d13c.year)]);
    legend('Full corrected','Full Original','location','best')
    legend boxoff
end
 print('-dpsc','-append',graphname)

if verboseMode
    details(iter).original = original;
    details(iter).fin = fin;
    details(iter).d13c = d13c;
    details(iter).p = pp;
    details(iter).ci = ci;
    details(iter).ca = ca;
end

clear delta cLow yhat pp

end  % termination of looping over different series, in columns

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Nested Subfunction for carbon dioxide concentration of the atmosphere

function [ca] = co2(iyear,jyear);

% CO2 The smoothed \emph{ca} curve for the appropriate time period
%
% Data are from Robertson et al. (2001) and McCarroll et al. (2007) and 
% based on Mauna Loa data from http://cdiac.ornl.gov/ftp/trends/co2/maunaloa.co2

caData = [...
1820    285 
1821    285 
1822    285 
1823    285 
1824    285 
1825    285 
1826    285 
1827    285 
1828    285 
1829    285 
1830    285 
1831    285 
1832    285 
1833    285 
1834    285 
1835    285 
1836    285 
1837    285 
1838    285 
1839    285 
1840    285 
1841    285 
1842    285 
1843    285 
1844    285 
1845    285 
1846    285.038334  
1847    285.102512  
1848    285.1726304 
1849    285.2487062 
1850    285.3307562 
1851    285.4187974 
1852    285.5128465 
1853    285.6129206 
1854    285.7190366 
1855    285.8312112 
1856    285.9490543 
1857    286.0722687 
1858    286.2010108 
1859    286.3354374 
1860    286.4757049 
1861    286.6219698 
1862    286.7743887 
1863    286.9331181 
1864    287.0983146 
1865    287.2701347 
1866    287.448735  
1867    287.637783  
1868    287.8398015 
1869    288.0532294 
1870    288.2765058 
1871    288.5080696 
1872    288.74636   
1873    288.9898159 
1874    289.2368764 
1875    289.4859805 
1876    289.7355671 
1877    289.9840754 
1878    290.2299444 
1879    290.480424  
1880    290.7425957 
1881    291.0146454 
1882    291.2947592 
1883    291.5811232 
1884    291.8719234 
1885    292.1653459 
1886    292.4595766 
1887    292.7528017 
1888    293.0432072 
1889    293.3289792 
1890    293.6160173 
1891    293.9106839 
1892    294.2118587 
1893    294.5184218 
1894    294.829253  
1895    295.1432322 
1896    295.4592392 
1897    295.776154  
1898    296.0928564 
1899    296.4082263 
1900    296.7211436 
1901    297.0304882 
1902    297.3414448
1903    297.6590506 
1904    297.9819645 
1905    298.3088451 
1906    298.6383513 
1907    298.9691419 
1908    299.2998756 
1909    299.6292112 
1910    299.9558075 
1911    300.2783233 
1912    300.5954172 
1913    300.9096966 
1914    301.2245886 
1915    301.539982  
1916    301.8557655 
1917    302.1718278 
1918    302.4880576 
1919    302.8043437 
1920    303.1205749 
1921    303.4366397 
1922    303.7524271 
1923    304.0678256 
1924    304.3827241 
1925    304.6864203 
1926    304.9711949 
1927    305.2414103 
1928    305.5014292 
1929    305.7556141 
1930    306.0083276 
1931    306.2639323 
1932    306.5267908 
1933    306.8012656 
1934    307.0917194 
1935    307.4025146 
1936    307.7380141 
1937    308.0644948 
1938    308.3520276 
1939    308.6106653 
1940    308.8504606 
1941    309.0814663 
1942    309.3137352 
1943    309.55732   
1944    309.8222737 
1945    310.1186489 
1946    310.4564984 
1947    310.845875  
1948    311.2968316 
1949    311.7636189 
1950    312.2012998 
1951    312.621145  
1952    313.034425  
1953    313.4524106 
1954    313.8863723 
1955    314.3475809 
1956    314.847307  
1957    315.3968211 
1958    316.007394  
1959    316.6902964 
1960    317.4253819 
1961    318.1854243 
1962    318.9710742 
1963    319.7829819 
1964    320.6217981 
1965    321.4881731 
1966    322.3827576 
1967    323.306202  
1968    324.2591568 
1969    325.2422725 
1970    326.2561997 
1971    327.3015888 
1972    328.375894  
1973    329.4764591 
1974    330.6037693 
1975    331.7583097 
1976    332.9405655 
1977    334.1510219 
1978    335.3901641 
1979    336.6584771 
1980    337.9564461 
1981    339.2845563 
1982    340.6432929 
1983    342.031623  
1984    343.4482469 
1985    344.8932493 
1986    346.3667151 
1987    347.868729  
1988    349.3993757 
1989    350.9587401 
1990    352.5469069 
1991    354.1639609 
1992    355.8099868 
1993    357.4850694 
1994    359.1892934 
1995    360.9227222 
1996    362.685303  
1997    364.476947  
1998    366.2975655 
1999    368.1470695 
2000    370.0253702 
2001    371.9323788 
2002    373.8680065
2003    375.83
2004    377.82
2005    379.85];

iyear = find(caData(:,1) >= iyear & caData(:,1) <= jyear); 
warning off
ca.year = caData(iyear,1); ca.ratio = caData(iyear,2);

[Marc Jakobi]

Please post Matlab code formatted as code. It is currently not readable.

The problem is that it is a function and not a script. You cannot run a function like you would run a script by pressing the run button. You should call the function from within a script or from the command window as:

[corrected, details] = pin(year,d13c,alpha,varargin)

where year, d13c and alpha are your data.

Type

doc pin

in the command window to find out how to use it.

[Walter Roberson]

All of the code starting from 'function' should be stored in pin.m . You can then call upon pin from the command line or from a script. You cannot paste the code in at the command prompt.

If you are using any version before R2016b then it is not permitted to store a function after a script in the same file.

If you are using the new R2016b then it is permitted to store a function after a script in the same file. In such a case the .m file name must not be the name of any function that is stored in the .m file.

[Marc Jakobi]

Nice, I was not aware of that new feature.