用多元指数多项式去拟合一组数据，如何加上x在0到40，y在0到360积分为1这个约束

下图为多元指数多项式，去拟合一组数据（x,y,f(x,y)），约束条件为x在0到40，y在0到360（xy的范围最好可变），f(x,y)的积分等于1。
一维的指数多项式，拟合一组数据，积分为1的方法，可以用非线性方程组方法和比例系数法，为下面2个程序。然后二维数据的多元指数多项式不带积分为1的程序也在下面写了，如果把一维指数多项式的非线性方程组方法和比例系数法用在二维数据里，要如何编写，写了一半的程序和仿造的程序都在下面。

一维指数多项式积分为1，非线性方程组法

x = [0:1:40]; % x随意
mu = 20; sigma = 10;
y = 1/sqrt(2*pi)/sigma*exp(-(x-mu+rand(size(x))*2).^2/2/sigma^2);% y随意设置，在正弦函数上加随机浮动
y(y<=0) = 10*eps;
n = 5; % n随意
x = x(:);
y = y(:);
S = zeros(length(x), n+1);
q = ones(length(x),1);
S(:,n+1)=q;
for i = n:-1:1
    q = q.*x;
    S(:,i) = q;
end
logy = log(y);
A = S'*S;
b = S'*logy;
a = A\b;%最小二乘
ff = @(x,a) exp(polyval(a,x));
eq = @(a)[A*a-b; integral(@(x)ff(x,a),0,40)-1];%构建了方程组后面一项指的是积分为1
options = optimoptions('fsolve','Display','off','FunctionTolerance',1e-2,'algorithm','levenberg-marquardt');
[a, f, flag] = fsolve(eq, a, options);
while (flag>=4 || flag<=0)
    [aa, f, flag] = fsolve(eq, a+a.*(2*rand(size(a))-1), options);
end
a = aa;
xfit = linspace(min(x),max(x),1001);
plot(x,y,'r-',xfit,exp(polyval(a,xfit)),'b-')
legend('原值','拟合后')

一维指数多项式积分为1，比例系数法

x = [0:1:40]; % x随意
mu = 20; sigma = 10;
y = 1/sqrt(2*pi)/sigma*exp(-(x-mu+rand(size(x))*2).^2/2/sigma^2);% y随意设置，在正弦函数上加随机浮动
y(y<=0) = 10*eps;
n = 5; % n随意
x = x(:);
y = y(:);
S = zeros(length(x), n+1);
q = ones(length(x),1);
S(:,n+1)=q;
for i = n:-1:1
    q = q.*x;
    S(:,i) = q;
end
logy = log(y);
A = S'*S;
b = S'*logy;
a = A\b;%最小二乘
ff = @(x,a) exp(polyval(a,x));
S = integral(@(x)ff(x,a),0,40);
a(end) = a(end) - log(S);%整体除以S
xfit = linspace(min(x),max(x),1001);
plot(x,y,'r-',xfit,exp(polyval(a,xfit)),'b-')
legend('原值','拟合后')

多元指数多项式，不带积分为1的约束程序

%% 为了得到分布，假设了一个协方差矩阵
mu=[20,180];%数学期望
sigma=[20/3 0;0,60].^2;%协方差矩阵
r=mvnrnd(mu,sigma,100000);%生成100000个样本
x = r(:,1);
y = r(:,2);
nx = 40;
ny = 40;
xmax = 40;
ymax = 360;
dx = xmax/nx;
dy = ymax/ny;
[X, Y, C, xmid, ymid] = ef2(x,y,nx+1,ny+1,[0,xmax],[0,ymax]);
C = C/(dx*dy);%单位面积上的概率即概率密度
figure(1);clf;
bar3(C)
title('原先数据')
xtick = 1:size(C,2); xticklabel = xmid; % xtick和xticklabel一定要对应，长度相等
ytick = 1:size(C,1); yticklabel = ymid; % ytick和yticklabel一定要对应，长度相等
set(gca, 'xtick', xtick, 'xticklabel',xticklabel,...
    'ytick', ytick, 'yticklabel',yticklabel)
xlabel('X')
ylabel('Y')
x = X(:);
y = Y(:);
z = C(:);%根据x和y算出的频率，x，y，z为输入数据
n = 7; % n为多项式的阶数，参数一共为(n+1)²个
[i1, i2] = meshgrid(0:n);
F = arrayfun(@(i)x.^i1(i).*y.^i2(i), 1:numel(i1), 'uniform', 0);
F = cell2mat(F);

F1=[z,F]%这里是因为z数据有时候频率为0，必须把含0的数据去掉
[r,c]=size(F1);
index=1:r; %一维矢量, F1的行指标
all(F1') %将F1转置一下, 返回一个矢量,它的每个元素进行判断, F1所在行的元素全不为0则是1, 否则为0
F1=F1(index(all(F1')),:) %取出F1不含0元素的行,构成新的矩阵
F1(:,1)=[]
z(z==0)=[]
logz=log(z)

A = F1'*F1;
b = F1'*logz;
a = A\b;%最小二乘法
zfit = polyfunval(X,Y,a,n);
zfit1 = exp(zfit);
surf(X,Y,zfit1);%这之后的非线性方程组法和比例系数法如何仿造上面的程序继续

function f = polyfunval(x,y,a,n)%指数多元多项式的公式
[i1, i2] = meshgrid(0:n);
Cfit = arrayfun(@(i)a(i)*x.^i1(i).*y.^i2(i), 1:numel(i1), 'uniform', 0);
f = Cfit{1};
for i = 2:numel(Cfit)
    f = f + Cfit{i};
end
end

function [X, Y, CDF, xmid, ymid] = ef2(x,y,nx,ny,xminmax,yminmax)%输入数据，只知道x和y，算出频率z
num = length(x);
if(num~=length(y))
   error('输入的x和y长度必须相等') 
end
if(nargin>6) % 如果变量个数大于6个，太多了
    error('太多输入变量')
elseif(nargin<2) % 如果变量个数小于2个，太少了
    error('输入变量数目不足!!')
end
if(nargin==6) % 如果变量个数等于6个，赋值给ymin和ymax
    ymin = yminmax(1);
    ymax = yminmax(2);
end
if(nargin>=5)% 如果变量个数大于等于5个，赋值给xmin和xmax
    xmin = xminmax(1);
    xmax = xminmax(2);
end
if(nargin<=4)% 如果变量个数小于等于4个，自定义xmin和xmax
    xmin = min(x);
    xmax = max(x)+eps;
    ymin = min(y);
    ymax = max(y)+eps;
end
if(nargin<=3)% 如果变量个数小于等于3个，自定义y方向划分段数ny
    ny = 30;
end
if(nargin==2)% 如果变量个数等于2个，自定义x方向划分段数nx
    nx = 30;
end
xg = linspace(xmin, xmax, nx);%x方向的点
yg = linspace(ymin, ymax, ny);%y方向的点
xmid = (xg(1:end-1)+xg(2:end))/2;
ymid = (yg(1:end-1)+yg(2:end))/2;
[X,Y] = meshgrid(xmid, ymid);%形成网格
[I,J] = meshgrid(2:nx, 2:ny);%下标网格
CDF = arrayfun(@(i,j)sum(x>=xg(i-1)&x<xg(i)&y<yg(j)&y>=yg(j-1))/num,I,J);%形成经验分布
end

你好，这里面直接构造最小二乘法然后用比例系数比较方便，建议n不要取太大，容易造成过拟合，使得拟合误差过大：

%% 为了得到分布，假设了一个协方差矩阵
mu=[20,180];%数学期望
sigma=[20/3 0;0,60].^2;%协方差矩阵
r=mvnrnd(mu,sigma,100000);%生成100000个样本
x = r(:,1);
y = r(:,2);
nx = 40;
ny = 40;
xmax = 40;
ymax = 360;
dx = xmax/nx;
dy = ymax/ny;
[X, Y, C, xmid, ymid] = ef2(x,y,nx+1,ny+1,[0,xmax],[0,ymax]);
C = C/(dx*dy);%单位面积上的概率即概率密度
figure(1);clf;
bar3(C)
title('原先数据')
xtick = 1:size(C,2); xticklabel = xmid; % xtick和xticklabel一定要对应，长度相等
ytick = 1:size(C,1); yticklabel = ymid; % ytick和yticklabel一定要对应，长度相等
set(gca, 'xtick', xtick(1:4:end), 'xticklabel',xticklabel(1:4:end),...
    'ytick', ytick(1:4:end), 'yticklabel',yticklabel(1:4:end))
xlabel('X')
ylabel('Y')
x = X(:);
y = Y(:);
z = C(:);%根据x和y算出的频率，x，y，z为输入数据
n = 4; % n为多项式的阶数，参数一共为(n+1)²个, 注意n不要取太大，容易产生过拟合
[i1, i2] = meshgrid(0:n);
F = arrayfun(@(i)x.^i1(i).*y.^i2(i), 1:numel(i1), 'uniform', 0);
F = cell2mat(F);
%% 标注------------改动的地方-------------------
z(z==0) = eps; %如果z等于零，赋值给一个很小的数字即可
logz=log(z);
A = F'*F;
b = F'*logz;
a = A\b;%最小二乘法
zfit = polyfunval(X,Y,a,n);
zfun = @(x,y)exp(polyfunval(x,y,a,n)); % 构建z的函数获取整体积分值
V = integral2(zfun, 0, 40, 0, 360);
a(1) = a(1) - log(V);
zfun = @(x,y)exp(polyfunval(x,y,a,n)); % 构建z的函数获取整体积分值,修正了a，使得积分为1
V_after_correct = integral2(zfun,0, 40, 0, 360) %这里可以显示积分值
zfit1 = exp(zfit);
figure(2);clf
surf(X,Y,zfit1, 'facealpha', 0.8);%这之后的非线性方程组法和比例系数法如何仿造上面的程序继续

function f = polyfunval(x,y,a,n)%指数多元多项式的公式
[i1, i2] = meshgrid(0:n);
Cfit = arrayfun(@(i)a(i)*x.^i1(i).*y.^i2(i), 1:numel(i1), 'uniform', 0);
f = Cfit{1};
for i = 2:numel(Cfit)
    f = f + Cfit{i};
end
end
function [X, Y, CDF, xmid, ymid] = ef2(x,y,nx,ny,xminmax,yminmax)%输入数据，只知道x和y，算出频率z
num = length(x);
if(num~=length(y))
   error('输入的x和y长度必须相等') 
end
if(nargin>6) % 如果变量个数大于6个，太多了
    error('太多输入变量')
elseif(nargin<2) % 如果变量个数小于2个，太少了
    error('输入变量数目不足!!')
end
if(nargin==6) % 如果变量个数等于6个，赋值给ymin和ymax
    ymin = yminmax(1);
    ymax = yminmax(2);
end
if(nargin>=5)% 如果变量个数大于等于5个，赋值给xmin和xmax
    xmin = xminmax(1);
    xmax = xminmax(2);
end
if(nargin<=4)% 如果变量个数小于等于4个，自定义xmin和xmax
    xmin = min(x);
    xmax = max(x)+eps;
    ymin = min(y);
    ymax = max(y)+eps;
end
if(nargin<=3)% 如果变量个数小于等于3个，自定义y方向划分段数ny
    ny = 30;
end
if(nargin==2)% 如果变量个数等于2个，自定义x方向划分段数nx
    nx = 30;
end
xg = linspace(xmin, xmax, nx);%x方向的点
yg = linspace(ymin, ymax, ny);%y方向的点
xmid = (xg(1:end-1)+xg(2:end))/2;
ymid = (yg(1:end-1)+yg(2:end))/2;
[X,Y] = meshgrid(xmid, ymid);%形成网格
[I,J] = meshgrid(2:nx, 2:ny);%下标网格
CDF = arrayfun(@(i,j)sum(x>=xg(i-1)&x<xg(i)&y<yg(j)&y>=yg(j-1))/num,I,J);%形成经验分布
end

直接求解非线性方程组无法搜索到合理解，在这里给出，只作为参考

%% 为了得到分布，假设了一个协方差矩阵
mu=[20,180];%数学期望
sigma=[20/3 0;0,60].^2;%协方差矩阵
r=mvnrnd(mu,sigma,100000);%生成100000个样本
x = r(:,1);
y = r(:,2);
nx = 40;
ny = 40;
xmax = 40;
ymax = 360;
dx = xmax/nx;
dy = ymax/ny;
[X, Y, C, xmid, ymid] = ef2(x,y,nx+1,ny+1,[0,xmax],[0,ymax]);
C = C/(dx*dy);%单位面积上的概率即概率密度
figure(1);clf;
bar3(C)
title('原先数据')
xtick = 1:size(C,2); xticklabel = xmid; % xtick和xticklabel一定要对应，长度相等
ytick = 1:size(C,1); yticklabel = ymid; % ytick和yticklabel一定要对应，长度相等
set(gca, 'xtick', xtick(1:4:end), 'xticklabel',xticklabel(1:4:end),...
    'ytick', ytick(1:4:end), 'yticklabel',yticklabel(1:4:end))
xlabel('X')
ylabel('Y')
x = X(:);
y = Y(:);
z = C(:);%根据x和y算出的频率，x，y，z为输入数据
n = 4; % n为多项式的阶数，参数一共为(n+1)²个, 注意n不要取太大，容易产生过拟合
[i1, i2] = meshgrid(0:n);
F = arrayfun(@(i)x.^i1(i).*y.^i2(i), 1:numel(i1), 'uniform', 0);
F = cell2mat(F);
%% 标注------------改动的地方-------------------
z(z==0) = eps; %如果z等于零，赋值给一个很小的数字即可
logz=log(z);
A = F'*F;
b = F'*logz;
a = A\b;%最小二乘法

ff = @(x,y,a)exp(polyfunval(x,y,a,n)); % 构建z的函数获取整体积分值
eq = @(a)[A*a-b; integral2(@(x,y)ff(x,y,a),0,40,0,360)-1];%构建了方程组后面一项指的是积分为1
options = optimoptions('fsolve','Display','off','FunctionTolerance',1e-2,'algorithm','levenberg-marquardt');
[aa, f, flag] = fsolve(eq, a, options);
count = 0;
while (flag>=4 || flag<=0)
    count = count + 1;
    fprintf('第%d次寻根\n',count);
    [aa, f, flag] = fsolve(eq, a+a.*(0.0002*rand(size(a))-1), options);
    if(count>100)
        fprintf('寻根未成功\n')
        aa = a;
        break;
    end
end
a = aa;
zfit = polyfunval(X,Y,a,n);
V_after_correct = integral2(zfun, min(x), max(x), min(y), max(y)) %这里可以显示积分值
zfit1 = exp(zfit);
figure(2);clf
surf(X,Y,zfit1, 'facealpha', 0.8);%这之后的非线性方程组法和比例系数法如何仿造上面的程序继续

function f = polyfunval(x,y,a,n)%指数多元多项式的公式
[i1, i2] = meshgrid(0:n);
Cfit = arrayfun(@(i)a(i)*x.^i1(i).*y.^i2(i), 1:numel(i1), 'uniform', 0);
f = Cfit{1};
for i = 2:numel(Cfit)
    f = f + Cfit{i};
end
end
function [X, Y, CDF, xmid, ymid] = ef2(x,y,nx,ny,xminmax,yminmax)%输入数据，只知道x和y，算出频率z
num = length(x);
if(num~=length(y))
   error('输入的x和y长度必须相等') 
end
if(nargin>6) % 如果变量个数大于6个，太多了
    error('太多输入变量')
elseif(nargin<2) % 如果变量个数小于2个，太少了
    error('输入变量数目不足!!')
end
if(nargin==6) % 如果变量个数等于6个，赋值给ymin和ymax
    ymin = yminmax(1);
    ymax = yminmax(2);
end
if(nargin>=5)% 如果变量个数大于等于5个，赋值给xmin和xmax
    xmin = xminmax(1);
    xmax = xminmax(2);
end
if(nargin<=4)% 如果变量个数小于等于4个，自定义xmin和xmax
    xmin = min(x);
    xmax = max(x)+eps;
    ymin = min(y);
    ymax = max(y)+eps;
end
if(nargin<=3)% 如果变量个数小于等于3个，自定义y方向划分段数ny
    ny = 30;
end
if(nargin==2)% 如果变量个数等于2个，自定义x方向划分段数nx
    nx = 30;
end
xg = linspace(xmin, xmax, nx);%x方向的点
yg = linspace(ymin, ymax, ny);%y方向的点
xmid = (xg(1:end-1)+xg(2:end))/2;
ymid = (yg(1:end-1)+yg(2:end))/2;
[X,Y] = meshgrid(xmid, ymid);%形成网格
[I,J] = meshgrid(2:nx, 2:ny);%下标网格
CDF = arrayfun(@(i,j)sum(x>=xg(i-1)&x<xg(i)&y<yg(j)&y>=yg(j-1))/num,I,J);%形成经验分布
end

总体来讲，比例系数法会更加合理且简单易操作，但是要注意n的取值不要过大（建议5以内），根据自己的数据调整n的取值会好一些

拟合问题：
$$
f(x,y)=exp(\sum\limits_{i,j}{a_{ij}x^iy^j})
$$
$$
\int_{x_{min}}^{x_{max}}\int_{y_{min}}^{y_{max}}{f(x,y)dxdy}=1
$$

抱歉没有细看你的第三段程序。参考比例系数法，如果你的第三段程序没问题的话，在你的第三个程序43行后面加上

S = integral2(@(x,y)polyfunval(x,y,a,n),0,40,0,360); % 求此二元函数在定义域内的二重积分
a(end) = a(end) - log(S);% 指数多项式中的常数项的系数更新，相当于整体函数除以S，积分因此变成1

注意这两行代码中第二行，需要更新a中对应多项式系数中的常数项。我不确定你的代码中是a(1)还是a(end)。

仅供参考，希望对你有所帮助。