按照如图的公式，如何用fmincon把约束加上

Question

我的代码是用极大似然估计来把下图的公式里面x,y,w参数算出来，w是权重，有约束，请问如何把打红框w的约束加进代码

如果i等于1的话，是如下代码，没有错误

X = xlsread('1.xlsx'); % X是个列向量
I0 = @(y) 1/(2*pi)*integral(@(x)exp(y.*cos(x)),0,2*pi);
% modelFun中, beta(1)对应y, beta(2)对应z
f = @(beta, x)1/(2*pi*I0(beta(1)))*exp(beta(1)*cos(x-beta(2)));
L = @(beta) -sum(log(f(beta, X))); % L取极小,-sum(log(f(beta, x)))就取极大
x0 = [1,1];
beta = fmincon(L, x0); % 这样beta就出来了
y = beta(1)
z = beta(2)

但当i等于4的时候，需要算w的约束，请问如何把w的约束加上去，下面是我暂时写出的代码，感觉有许多错误和可以改进的地方

X=xlsread('1.xlsx'); %x是一堆角度列向量

I0 = @(y) 1/(2*pi)*integral(@(x)exp(y(1).*cos(x)),0,2*pi);
I1 = @(y) 1/(2*pi)*integral(@(x)exp(y(2).*cos(x)),0,2*pi);
I2 = @(y) 1/(2*pi)*integral(@(x)exp(y(3).*cos(x)),0,2*pi);
I3 = @(y) 1/(2*pi)*integral(@(x)exp(y(4).*cos(x)),0,2*pi);%这一块感觉写的有问题

f = @(beta, x)beta(3)*(1/(2*pi*I0(beta(1)))*exp(beta(1)*cos(x-beta(2))))...
+beta(6)*(1/(2*pi*I1(beta(4)))*exp(beta(4)*cos(x-beta(5))))...
+beta(9)*(1/(2*pi*I2(beta(7)))*exp(beta(7)*cos(x-beta(8))))...
+beta(10)*(1/(2*pi*I3(beta(11)))*exp(beta(11)*cos(x-beta(12))));

L = @(beta) -sum(log(f(beta, X)));
% 如何在这里把w的约束加进去
x0 = [1,1];
beta = fmincon(L, x0);

y(1) = beta(1),y(2) = beta(4),y(3) = beta(7),y(4) = beta(10)
z(1) = beta(2),z(2) = beta(5),z(3) = beta(8),z(4) = beta(11)
w(1) = beta(3),w(2) = beta(6),w(3) = beta(9),w(4) = beta(12)

Answer 1

可以这么写，你看看

clc; clear
% X = rand(100,1); % 我是用随机数做测试的
X = xlsread('1.xlsx'); % X是个列向量
f = @(beta, x)fun(x, beta);
L = @(beta) -sum(log(f(beta, X))); % L取极小,-sum(log(f(beta, x)))就取极大
n = 5; % 比方说有y、z、w各n个
beta0 = rand(n,3);
% Aeq 和beq保证wi相加为1
Aeq = [zeros(1, 2*n), ones(1, n)];
beq = 1;
% lb ub保证wi在0~1之间
lb = [-inf*ones(2*n,1); zeros(n,1)];
ub = [ inf*ones(2*n,1);  ones(n,1)];
beta = fmincon(L, beta0(:),[],[],Aeq,beq,lb,ub); % 这样beta就出来了
beta = reshape(beta, n, 3);
y = beta(:,1)
z = beta(:,2)
w = beta(:,3)

function f = fun(x, yzw)
n = numel(yzw);
m = n/3;
y = yzw(1:m);
z = yzw(m+1:2*m);
w = yzw(2*m+1:3*m);
I0 = @(y) 1/(2*pi)*integral(@(x)exp(y.*cos(x)),0,2*pi);
f = zeros(size(x));
for i = 1:m
    f = f + w(i)*1/(2*pi*I0(y(i)))*exp(y(i)*cos(x-z(i)));
end
end