Estimation of parameters Problem

clc;
clear all

t = [1 2 4 6 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262];
c = [0 81 71 10 0 00 106 89 17 0 00 115 96 19 0 00 163 142 21 0 00 267 240 27 0 00 273 244 29 0 00 274 245 29 0 00 276 247 29 0 00 287 257 30 0 00 290 259 31 0 00 298 266 32 0 00 302 269 33 0 00 313 279 34 0 00 326 291 35 0 00 394 345 49 0 00 462 398 64 0 00 501 429 72 0 00 571 480 91 0 00 579 487 92 0 00 603 505 98 0 00 608 508 100 0 00 615 511 104 0 00 620 515 105 0 00 633 525 108 0 00 634 526 108 0 00 635 527 108 0 00 639 528 111 0 00 642 531 111 0 00 652 541 111 0 00 669 549 120 0 00 704 577 127 0 00 716 586 130 0 00 727 593 134 0 00 734 600 134 0 00 737 603 134 0 00 746 609 137 0 00 753 615 138 0 00 758 618 140 0 00 768 625 143 0 00 845 696 149 0 00 862 707 155 0 00 939 773 166 0 00 969 797 172 0 00 981 806 175 0 00 990 813 177 0 00 994 817 177 0 00 998 820 178 0 00 1009 830 179 0 00 1018 838 180 0 00 1021 841 180 0 00 1026 846 180 0 00 1036 855 181 0 00 1040 857 183 0 00 1049 866 183 0 00 1053 870 183 0 00 1061 878 183 0 00 1074 887 187 0 00 1075 888 187 0 00 1088 893 195 0 00 1169 965 204 0 00 1170 965 205 0 00 1178 973 205 0 00 1183 976 207 0 00 1192 982 210 0 00 1197 986 211 0 00 1211 998 213 0 00 1259 1040 219 0 00 1368 1141 227 0 00 1369 1141 228 0 00 1376 1145 231 0 00 1381 1150 231 0 00 1387 1156 231 0 00 1393 1162 231 0 00 1403 1170 233 0 00 1413 1180 233 0 00 1450 1207 243 0 00 1496 1244 252 0 00 1548 1285 263 0 00 1586 1309 277 0 00 1592 1315 277 0 00 1610 1332 278 0 00 1618 1337 281 0 00 1623 1342 281 0 00 1630 1349 281 0 00 1635 1354 281 0 00 1636 1355 281 0 00 1640 1358 282 0 00 1652 1362 290 0 00 1716 1414 302 0 00 1735 1429 306 0 00 1745 1439 306 0 00 1749 1443 306 0 00 1757 1449 308 0 00 1760 1452 308 0 00 1765 1457 308 0 00 1770 1462 308 0 00 1775 1467 308 0 00 1785 1474 311 0 00 1922 1593 329 0 00 1941 1609 332 0 00 2032 1679 353 0 00 2064 1707 357 0 00 2093 1730 363 0 00 2095 1731 364 0 00 2102 1737 365 0 00 2106 1740 366 0 00 2112 1746 366 0 00 2117 1750 367 0 00 2120 1753 367 0 00 2140 1771 369 0 00 2145 1775 370 0 00 2150 1778 372 0 00 2159 1787 372 0 00 2161 1789 372 0 00 2168 1796 372 0 00 2172 1799 373 0 00 2189 1813 376 0 00 2240 1855 385 0 00 2248 1863 385 0 00 2258 1871 387 0 00 2271 1882 389 0 00 2273 1883 390 0 00 2290 1897 393 0 00 2364 1960 404 0 00 2479 2057 422 0 00 2484 2059 425 0 00 2494 2066 428 0 00 2500 2071 429 0 00 2501 2072 429 0 00 2507 2077 430 0 00 2515 2085 430 0 00 2528 2098 430 0 00 2545 2114 431 0 00 2570 2133 437 0 00 2629 2182 447 0 00 2652 2201 451 0 00 2691 2238 453 0 00 2698 2243 455 0 00 2716 2260 456 0 00 2722 2265 457 0 00 2732 2274 458 0 00 2736 2277 459 0 00 2738 2279 459 0 00 2740 2281 459 0 00 2761 2300 461 0 00 2767 2305 462 0 00 2768 2306 462 0 00 2777 2314 463 0 00 2790 2325 465 0 00 2794 2329 465 0 00 2808 2336 472 0 00 2826 2344 482 0 00 2842 2356 486 0 00 2851 2362 489 0 00 2857 2368 489 0 00 2864 2374 490 0 00 2877 2386 491 0 00 2886 2391 495 0 00 2891 2395 496 0 00 2896 2400 496 0 00 2909 2411 498 0 00 2912 2413 499 0 00 2941 2439 502 0 00 3058 2540 518 0 00 3148 2618 530 0 00 3155 2622 533 0 00 3177 2640 537 0 00 3180 2641 539 0 00 3184 2644 540 0 00 3196 2655 541 0 00 3202 2660 542 0 00 3205 2662 543 0 00 3215 2669 546 0 00 3220 2672 548 0 00 3229 2681 548 0 00 3236 2686 550 0 00 3246 2695 551 0 00 3249 2698 551 0 00 3277 2723 554 0 00 3287 2733 554 0 00 3290 2734 556 0 00 3300 2742 558 0 00 3307 2746 561 0 00 3321 2758 563 0 00 3363 2794 569 0 00 3465 2883 582 0 00 3471 2886 585 0 00 3482 2894 588 0 00 3486 2897 589 0 00 3488 2899 589 0 00 3496 2907 589 0 00 3508 2919 589 0 00 3531 2942 589 0 00 3571 2977 594 0 00 3586 2987 599 0 00 3609 3006 603 0 00 3617 3011 606 0 00 3625 3016 609 0 00 3634 3023 611 0 00 3643 3031 612 0 00 3647 3034 613 0 00 3657 3042 615 0 00 3660 3045 615 0 00 3666 3049 617 0 00 3669 3052 617 0 00 3672 3053 619 0 00 3694 3072 622 0 00 3716 3089 627 0 00 3811 3165 646 0 00 3865 3209 656 0 00 3883 3221 662 0 00 3891 3229 662 0 00 3898 3235 663 0 00 3901 3238 663 0 00 3911 3247 664 0 00 3915 3250 665 0 00 3943 3275 668 0 00 3970 3291 679 0 00 4076 3382 694 0 00 4101 3406 695 0 00 4158 3450 708 0 00 4161 3453 708 0 00 4171 3459 712 0 00 4180 3468 712 0 00 4183 3469 714 0 00 4190 3475 715 0 00 4196 3481 715 0 00 4202 3486 716 0 00 4217 3499 718 0 00 4221 3502 719 0 00 4226 3506 720 0 00 4230 3509 721 0 00 4239 3516 723 0 00 4255 3530 725 0 00 4263 3537 726 0 00 4266 3539 727 0 00 4268 3540 728 0 00 4281 3551 730 0 00 4289 3557 732 0 00 4305 3569 736 0 00 4313 3576 737 0 00 4332 3594 738 0 00 4338 3600 738 0 00 4342 3602 740 0 00 4344 3604 740 0 00 4355 3614 741 0 00 4361 3618 743 0 00 4373 3629 744 0 00 4383 3639 744 0 00 4389 3644 745 0 00 4393 3647 746 0 00 4400 3654 746 0 00 4404 3658 746 0 00 4462 3710 752 0 00 4465 3711 754 0 00 4478 3724 754 0 00 4491 3736 755 0 00 4495 3740 755 0 00 4512 3757 755 0 00 4542 3786 756 0 00 4565 3806 759 0 00 4571 3812 759 0 0]; 

theta0 = rand(8,1);

[theta, Rsdnrm, Rsd, ExFlg, OptmInfo, Lmda, Jmat] = lsqcurvefit(@kinetic, theta0, t, c, zeros(size(theta0))); 

fprintf(1, '\tRate Constants:\n')

for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%2d) = %8.5f\n', k1, theta(k1))
end 

tv = linspace(min(t), max(t));
Cfit = kinetic(theta, tv);
figurehd = plot(t, c, '--p');

for k1 = 1:size(c,2)
    Cv(k1,:) = hd(k1).Color; 
    hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);

for k1 = 1:size(c,2)
    hlp(k1).Color = CV(k1,:);
end
hold off

grid on
xlabel('Time')
ylabel('Concentration')
legend(hlp, compose('C_%d',1:size(c,2)), 'Location','best') 

% theta0 0.5948960740086140.2622117477808450.6028430893820830.7112157804336830.2217467340172400.1174176508558060.2966758732183270.318778301925882

function C = kinetic(theta, t) 
    c0      = theta(end-3:end); 
    [T, Cv] = ode45(@DifEq, t, c0);
        function dC = DifEq(t, c)
            dcdt    = zeros(6,1); 
            dcdt(1) = 1 - theta(1).*(c(1)*c(3)) - theta(2).*(c(1)*c(5))- theta(3).*c(1) + theta(4).*c(4); 
            dcdt(2) = theta(1).*(c(1)*c(3)) + theta(2).*(c(1)*c(5)) - (theta(5) + theta(3)).*c(2); 
            dcdt(3) = theta(5).*c(2) - (theta(6) + theta(7) + theta(3)).*c(3); 
            dcdt(4) = theta(6).*c(3) - (theta(4) + theta(3)).*c(4); 
            dcdt(5) = 1 - theta(2).*(c(5)*c(6)) - theta(8).*c(5); 
            dcdt(6) = theta(2).*(c(5)*c(6)) - theta(8).*c(6); 
            dC      = dcdt;
        end 
    C       = Cv;
end

Error???

Error using ==> lsqncommon at 101

LSQCURVEFIT cannot continue because user supplied objective function failed with the following error:Attempted to access c(5); index out of bounds because numel(c)=4.

Error in ==> lsqcurvefit at 182

[x,Resnorm,FVAL,EXITFLAG,OUTPUT,LAMBDA,JACOB] = ...

Error in ==> kl at 522

[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetic,theta0,t,c,zeros(size(theta0)));

Please I need help to fix the error and to get the following 1. Theta(1)-theta(8) 2. Display the graph 3. hd 4. Rsdnrm

[Dyuman Joshi]

Please format your code properly.

[Sunday Aloke]

I don't know how to go about it again

[Sam Chak]

@Sunday Aloke, try fixing the rest by pressing the "Enter" key to break the line.

Since you are in possession of the code, it is peanut and butter for you to compare each line of the code.

[Star Strider]

@Sunday Aloke — Thank you for quoting my code.

There are two serious problems —

(1) The ‘t’ variable makes no sense

(2) The ‘c’ variable should be a matrix

There needs to be one value of ‘t’ (that ideally should start from 0) for each row of ‘c’ for my code to work correctly.

[Sunday Aloke]

I've followed the instruction but am still having errors

?? Error using ==> lsqncommon at 101LSQCURVEFIT cannot continue because user supplied objective function failed with the following error:Subscript indices must either be real positive integers or logicals.Error in ==> lsqcurvefit at 182[x,Resnorm,FVAL,EXITFLAG,OUTPUT,LAMBDA,JACOB] = ...Error in ==> kl at 24[theta,Resnorm,FVAL,EXITFLAG,OUTPUT,LAMBDA,JACOB] =lsqcurvefit(@kinetic,theta0,t,c,zeros(size(theta0)));

[Torsten]

"t" has to be a vector of increasing values (usually starting with 0 as first element). "c" has to be a vector with as many elements as "t".

If you changed your code accordingly, please include it here so that we are able to execute it after your modifications.

[Walter Roberson]

t = [1246891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262];

t

t = Inf

That is a scalar. MATLAB has no idea where you want to divide those digits into separate numbers.

[Sam Chak]

@Sunday Aloke, Can you attach (click on the paperclip icon ) the m-file where the MATLAB code is stored? I think you want to estimate the coefficients in the ODE.

[Sunday Aloke]

I have a reported case for only the infected, recovered and death, which was collected for a period of 262 weeks. C(1)=susceptible humans C(2)=exposed humans C(3)=infected humans C(4)=recovered humans C(5)=susceptible rodents C(6)=infected rodents. I wish to estimate the parameters in the models which are theta(1)-theta(8).

I need to be guided on how to use either of the data to estimate the parameters.

[Sunday Aloke]

I have a system of 6 differential equations, which includes the susceptible humans, exposed humans, infected humans, recovered humans, susceptible rodents and infected rodents. I also have data on Infected humans, Death humans and recovered humans. Please, I need MATLAB code to write a program that will help me to estimate the parameters.

[Star Strider]

Seems like you’re modeling Yersinia pestis (although it could be any number of zoönoses, such as Hantavirus) and the ‘t’ variable still doesn’t make sense.

[Sam Chak]

@Star Strider, after studying the pattern of the values, I think t is a time vector.

t = [1 2 4 6 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262]

t = 1×256

1 2 4 6 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33

@Sunday Aloke, I have editted t in the question. If this is not true, please correct it.

[Star Strider]

@Sam Chak — Unfortunately, even with that change, it doesn’t work because there is no way to reshape ‘c’ to fit both ‘t’ and the 6 differential equations.

The data matrix has to have a row size the same size as ‘t’ and a column size the same size as the number of differential equations, here 6. So, in order for ‘c’ be be able to be reshaped correctly into a (256 x 6) matrix, it has to have:

256 * 6 = 1536

elements. It doesn’t.

There is nothing to do here until the data make sense.

[Sunday Aloke]

The reported data are only 3 which I've stated above.

[Torsten]

@Sunday Aloke

If you have 256 times when your three data (Infected humans, Death humans and recovered humans) are reported (the vector t is of length 256, not 262), you should get a vector c of length 3*256 = 768 or a matrix of size 3x256. Aren't you able to arrange your data this way or explain what your c-vector of length 1311 is supposed to represent (especially these strange 00 entries) ?

And you write C(1)=susceptible humans C(2)=exposed humans C(3)=infected humans C(4)=recovered humans C(5)=susceptible rodents C(6)=infected rodents. Where are the deads or how are they computed from C(1),...,C(6) ?

[Sunday Aloke]

@Torsen the death is not in c(1)-c(6). The main data is c(3)= infected and c(4)=recovered

[Torsten]

Then generate a time vector of length 262 and a matrix of size 2x262 where row 1 is the c(3) and row 2 is c(4) from your measurements at the times from your time vector, put them here and we can start answering your question.