怎么用MATLAB实现下述计算

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

T0＝44.73，19.97，17.39，17.9
T1＝44.93，19.48，17.74，17.85
T2＝46.35，19.29，17.3，17.07
T3＝44.84，19.25，18.26，17.65
上述数据分别为02，03，04，05年四个地区所占比例的数据，现在要计算一步转移概率矩阵P，即T0*P=T1,T1*P=T2,T2*P=T3，其约束条件是：
依据残差平方和最小的思想，
第一步：即t年的状态j的概率＝t-1年状态i的概率*从状态i向状态j的转移概率，然后对所有的i求和；
第二步：用t年的j状态的比例减去上一步所求的和，对这个差求平方，然后对所有的j求和，即求平方和，让其最小。
同时，要求所求概率矩阵的各行元素和为1，各元素大于0小于1。

用一个笨方法来解的话是这样：
设转移概率矩阵为：P11 P12 P13 P14
                                     P21 P22 P23 P24
                                     P31 P32 P33 P34
                                     P41 P42 P43 P44
t1年状态1的差为：44.93－(17.97*P21+17.39*P31+17.9*P41)
        状态2的差为：19.48－（44.73*P12+17.39*P32+17.9*P42）
    状态3的差为：17.74－（44.73*P13+19.97*P23+17.9*P43）
    状态4的差为：17.85－（44.73*P14+19.97*P24+17.39*P34）
t2年状态1的差为：46.35－(19.48*P21+17.74*P31+17.85*P41)
        状态2的差为：19.29－（44.93*P12+17.74*P32+17.85*P42）
    状态3的差为：17.3－（44.93*P13+19.48*P23+17.85*P43）
    状态4的差为：17.07－（44.93*P14+19.48*P24+17.74*P34）
t3年状态1的差为：44.84－(19.29*P21+17.3*P31+17.07*P41)
        状态2的差为：19.25－（46.35*P12+17.3*P32+17.07*P42）
    状态3的差为：18.26－（46.35*P13+19.29*P23+17.07*P43）
    状态4的差为：17.65－（46.35*P14+19.29*P24+17.3*P34）
然后对上面所有的差分别求平方，再求和，让这个和最小，即残差平方和最小。
用这个约束条件进行求解，应该怎么运用MATLAB写语言呢？
正确结果应该是
P( 1, 1)       0.2569220
P( 1, 2)       0.1734723
P( 1, 3)       0.3918807
p( 1, 4)       0.1777250
P( 2, 1)        0.000000
p( 2, 2)       0.5165667
P( 2, 3)        0.000000
P( 2, 4)       0.4834333
P( 3, 1)        1.000000
P( 3, 2)        0.000000
P( 3, 3)        0.000000
P( 3, 4)        0.000000
P( 4, 1)       0.9227927
P( 4, 2)       0.7720731E-01
P( 4, 3)        0.000000
P( 4, 4)        0.000000

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：matlab实现 MATLAB matla atlab Lab MATLAB

function piggin
x0 = ones(12,1)*0.25;
A = kron(eye(4), ones(1,3)); % inequality constraints
b = ones(4,1); % A x <= b
lb = 0*x0;  % lower bound
ub = 1+lb; % upper bound
[xopt,fval,exitflag] = fmincon(@myfun, x0,A,b,[],[],lb,ub);
xe = 1-A*xopt;  % pii = 1-sum_j(pij)
P = diag(xe);
xind =[2:5, 7:10, 12:15];
P(xind) = xopt

% object function  is here
function f = myfun(x)
xind =[2:5, 7:10, 12:15];
P = zeros(16,1);
P(xind) = x;
P1 = reshape(P,4,4);
T0=[44.73,19.97,17.39,17.9];
T1=[44.93,19.48,17.74,17.85];
T2=[46.35,19.29,17.3,17.07];
T3=[44.84,19.25,18.26,17.65];
TL =[ T1; T2; T3];
TR =[ T0; T1; T2];
f = sum(sum((TL-TR*P1).^2));

P =

         0    0.2047    0.3719    0.1453
    1.0000    0.2239    0.0462    0.5585
         0         0    0.5818         0
    0.0000    0.5714         0    0.2961

3# jyliao
请问解是不是不唯一？

4# wpcz_c
Yes!

5# jyliao
可是和LINGO做出的结果为什么不一样呢？
P( 1, 1)       0.2569220
P( 1, 2)       0.1734723
P( 1, 3)       0.3918807
p( 1, 4)       0.1777250
P( 2, 1)        0.000000
p( 2, 2)       0.5165667
P( 2, 3)        0.000000
P( 2, 4)       0.4834333
P( 3, 1)        1.000000
P( 3, 2)        0.000000
P( 3, 3)        0.000000
P( 3, 4)        0.000000
P( 4, 1)       0.9227927
P( 4, 2)       0.7720731E-01
P( 4, 3)        0.000000
P( 4, 4)        0.000000

1.  check your statement in "用一个笨方法来解",  it is not feasible for the answer you provided in #5.
2. If you include the Pii terms in residual, then the answer from LINGO is better than Matlab .

[xopt,fval,exitflag] = fmincon(@myfun, x0,A,b,[],[],lb,ub);
这个命令是什么？fval 是否应为feval？等号右边括号里面的那两个方框又应该是什么？

MATLAB
help fmincon

If there are five more Tk (T0...T4), then it will solve by least-squares approach.