"Index exceeds 37" when 35 variables?

So, I've been trying to do 3 variable modeling in MATLAB. I have four sets of data (3 variables and one answer), and I'm trying to bet an equation modeling the data. When running the section of the code pasted below, I keep getting "Index exceeds number of array elements. Index must not exceed 37". There are only 35 variables in the code, so I'm not sure why this isn't working. As far as why I chose this format, I think x and z can be modeled by a quadratic individually, but I figured each variable needed to be the same power, so they're all 4th power. Any help would be greatly appreciated.

ft = fittype([ ...
    'p000 + p100*x + p010*y + p001*z + ' ...
    'p200*x.^2 + p020*y.^2 + p002*z.^2 + ' ...
    'p110*x*y + p101*x*z + p011*y*z + ' ...
    'p300*x.^3 + p030*y.^3 + p003*z.^3 + ' ...
    'p210*x.^2*y + p201*x.^2*z + p120*x*y.^2 + ' ...
    'p021*y.^2*z + p102*x*z.^2 + p012*y*z.^2 + ' ...
    'p111*x*y*z + ' ...
    'p400*x.^4 + p040*y.^4 + p004*z.^4 + ' ...
    'p310*x.^3*y + p301*x.^3*z + p130*x*y.^3 + ' ...
    'p031*y.^3*z + p103*x*z.^3 + p013*y*z.^3 + ' ...
    'p220*x.^2*y.^2 + p202*x.^2*z.^2 + p022*y.^2*z.^2 + ' ...
    'p211*x.^2*y*z + p121*x*y.^2*z + p112*x*y*z.^2'], ...
    'independent', {'x', 'y', 'z'}, ...
    'dependent', 'v');

[fitobject, gof, output] = fit([x, y, z], v, ft);

disp(fitobject)

[Matt J]

So, I've been trying to do 3 variable modeling in MATLAB. I have four sets of data (3 variables and one answer)

That's a non-starter. The fit() function doesn't support more than 2 independent variables.

It's also not a realistic model to attempt a fitting. A fully general 4th order 3D polynomial (35 coefficients) will be quite ill-conditioned.

Answer 1

[Matt J]

The root cause, as I commented above, is probably that you've attempted to construct a fit type with 3 independent variables, which is not allowed. Admittedly, the error message this returns is pretty uninformative, so it would be an appropriate bug report or enhancement request.

ft=fittype(@(a,b,c,x,y) a*x+b*y+c ,'independent', {'x', 'y'}, 'dependent','v')

ft = 
     General model:
     ft(a,b,c,x,y) = a*x+b*y+c

ft=fittype(@(a,b,c,x,y,z) a*x+b*y+c*z ,'independent', {'x', 'y','z'}, 'dependent','v')

Error using fittype/testCustomModelEvaluation (line 12)
Expression a*x+b*y+c*z is not a valid MATLAB expression, has non-scalar coefficients, or cannot be evaluated:
Index exceeds the number of array elements. Index must not exceed 5.

Error in fittype>iCreateFittype (line 367)
testCustomModelEvaluation( obj );
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in fittype (line 324)
obj = iCreateFittype( obj, varargin{:} );
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Caused by:
Index exceeds the number of array elements. Index must not exceed 5.

[Matt J]

As I mentioned also, your model needs to be reconsidered. Here is an AI-generated demo for 3D fitting with Chebyshev polynomials (well known to be more stable than conventional polynomials).

%% ===============================
%  Chebyshev 3D Fit Demo (Self-contained)
% ================================

% -----------------------------
% 1) Generate sample data
% -----------------------------
N = 2000;

x = 2*rand(N,1) - 1;
y = 2*rand(N,1) - 1;
z = 2*rand(N,1) - 1;

% True function
v_true = 1 + 2*x - y + 0.5*z.^2 + x.*y - 0.3*x.^2.*z;

% Add noise
v = v_true + 0.05*randn(size(v_true));

% -----------------------------
% 2) Fit model
% -----------------------------
degree = 4;
lambda = 1e-3;

model = fitCheb3D(x, y, z, v, degree, lambda);

% -----------------------------
% 3) Training error
% -----------------------------
v_pred = model.predict(x, y, z);

rmse = sqrt(mean((v - v_pred).^2));
fprintf('Training RMSE: %.4f\n', rmse);

Training RMSE: 0.0504

% -----------------------------
% 4) Test on new data
% -----------------------------
Nt = 500;

xt = 2*rand(Nt,1) - 1;
yt = 2*rand(Nt,1) - 1;
zt = 2*rand(Nt,1) - 1;

v_true_test = 1 + 2*xt - yt + 0.5*zt.^2 + xt.*yt - 0.3*xt.^2.*zt;

v_pred_test = model.predict(xt, yt, zt);

rmse_test = sqrt(mean((v_true_test - v_pred_test).^2));
fprintf('Test RMSE: %.4f\n', rmse_test);

Test RMSE: 0.0074

% -----------------------------
% 5) Visualization (slice z=0 with data overlay)
% -----------------------------
[xg, yg] = meshgrid(linspace(-1,1,50), linspace(-1,1,50));
zg = zeros(size(xg));

% Model prediction
vg_pred = model.predict(xg(:), yg(:), zg(:));
vg_pred = reshape(vg_pred, size(xg));

figure;
surf(xg, yg, vg_pred)
hold on

% ---- Select data near z = 0
eps_z = 0.05;   % thickness of slice
idx = abs(z) < eps_z;

% Overlay data points
scatter3(x(idx), y(idx), v(idx), ...
    25, 'r', 'filled', 'MarkerEdgeColor','k');

xlabel('x'); ylabel('y'); zlabel('v')
title('Fitted surface with data overlay (z ≈ 0)')
shading interp
alpha(0.7)

legend('Fit surface','Data (z≈0)')
view(45,30)

%% ===============================
%  Function: fitCheb3D
% ===============================
function model = fitCheb3D(x, y, z, v, degree, lambda)

    if nargin < 5
        degree = 4;
    end
    if nargin < 6
        lambda = 1e-3;
    end

    x = x(:); y = y(:); z = z(:); v = v(:);

    % ---- Scale to [-1,1]
    scale.xmin = min(x); scale.xmax = max(x);
    scale.ymin = min(y); scale.ymax = max(y);
    scale.zmin = min(z); scale.zmax = max(z);

    x1 = scaleToUnit(x, scale.xmin, scale.xmax);
    y1 = scaleToUnit(y, scale.ymin, scale.ymax);
    z1 = scaleToUnit(z, scale.zmin, scale.zmax);

    % ---- Chebyshev basis
    Tx = chebpoly(x1, degree);
    Ty = chebpoly(y1, degree);
    Tz = chebpoly(z1, degree);

    % ---- Design matrix
    A = buildTensorBasis(Tx, Ty, Tz, degree);

    % ---- Ridge regression
    c = (A' * A + lambda * eye(size(A,2))) \ (A' * v);

    % ---- Store model
    model.c = c;
    model.degree = degree;
    model.scale = scale;

    model.predict = @(xq,yq,zq) predictCheb3D(xq,yq,zq,model);
end


%% ===============================
%  Prediction
% ===============================
function vpred = predictCheb3D(x, y, z, model)

    x = x(:); y = y(:); z = z(:);

    x1 = scaleToUnit(x, model.scale.xmin, model.scale.xmax);
    y1 = scaleToUnit(y, model.scale.ymin, model.scale.ymax);
    z1 = scaleToUnit(z, model.scale.zmin, model.scale.zmax);

    Tx = chebpoly(x1, model.degree);
    Ty = chebpoly(y1, model.degree);
    Tz = chebpoly(z1, model.degree);

    A = buildTensorBasis(Tx, Ty, Tz, model.degree);

    vpred = A * model.c;
end


%% ===============================
%  Helpers
% ===============================
function x1 = scaleToUnit(x, xmin, xmax)
    x1 = 2*(x - xmin)/(xmax - xmin) - 1;
end

function T = chebpoly(x, n)
    N = numel(x);
    T = zeros(N, n+1);

    T(:,1) = 1;
    if n >= 1
        T(:,2) = x;
    end

    for k = 2:n
        T(:,k+1) = 2*x.*T(:,k) - T(:,k-1);
    end
end

function A = buildTensorBasis(Tx, Ty, Tz, degree)
    A = [];

    for i = 0:degree
        for j = 0:degree
            for k = 0:degree
                if i + j + k <= degree
                    A = [A, Tx(:,i+1).*Ty(:,j+1).*Tz(:,k+1)];
                end
            end
        end
    end
end

[John]

Thanks for this, I was actually able to get some results this time. Just to make sure I'm reading this right, the coefficients that appear in model.c are in what order? My current assumption is that it reads as (1)+(2)*(z)+(3)*(z^2)+(4)*(z^3)+(5)*(z^4)+(6)*(y)... Is that correct?

[Matt J]

The AI offers the example below, where T() is the Chebyshev polynomial of the first kind, and x',y',z' are the normalized x,y,z coordinates. It shouldn't really be something you need to delve into though. You should never be directly manipulating the coefficients in your work. You should be using the predictCheb3D function to evaluate the polynomials.