8
$\begingroup$

I know that there are already questions about fitting multiple datasets and about NDSolve and about shared and non shared parameters, but I tried to apply them and some things are still not clear.

Here is my equation :

l = 10^(-5)
k = 1/l
chic = 0.5
T = 100

eq = {R'[t] == -a[t]*R[t] + b[t], 
  b'[t] == beta/2*(Tanh[(chi[t] - chic)*k] - 1), 
  a'[t] == -alpha/2*(Tanh[(chi[t] - chic)*k] - 1), 
  chi'[t] == -kappa*R[t]*(chi[t] - 2*chic), a[0] == a0, b[0] == b0, 
  R[0] == R0, chi[0] == 0}

I want to fit with regards to the variables : $alpha, beta, kappa, a0, b0$ as shared parameters and $R0$ as non shared parameter, meaning it would be different fr each one.

The joined data is given as an appendix just afterwards.

The non-joined data (meaning the 5 data-sets separately) looks like that :

enter image description here

So I tried to change $R0$ as a variable, and I got inspired by the answer of @JimB in Finding NonlinearModelFit of multiple data sets with the same parameters and in two dimensions : 

model[alpha_?NumberQ, beta_?NumberQ, kappa_?NumberQ, a0_?NumberQ, 
  b0_?NumberQ] :=  (model[alpha, beta, kappa, a0, b0] = 
   Module[{R, chi, b, a, t, R0}, 
    First[R /. 
      NDSolve[{D[R[t, R0], t] == -a[t, R0]*R[t, R0] + 
          b[t, R0], 
        D[b[t, R0], t] == beta/2*(Tanh[(chi[t, R0] - chic)*k] - 1), 
        D[a[t, R0], t] == -alpha/2*(Tanh[(chi[t, R0] - chic)*k] - 1), 
        D[chi[t, R0], t] == -kappa*(chi[t, R0] - 2*chic), 
        a[0, R0] == a0, b[0, R0] == b0, R[0, R0] == R0, chi[0,R0] == 0}, {R, b, 
        a, chi}, {t, 0, T}, {R0, 0, 300}]]]);
nlm = NonlinearModelFit[data, 
   {model[alpha, beta, kappa, a0, b0][t, 
    R0], alpha >= 0, beta >= 0, kappa >= 0, a0 >= 0, b0 >= 0}, {{alpha, 0.1}, { beta, 0.1}, { kappa, 0.05}, {a0, 0.01}, {b0,
      3}}, {t, R0}];
nlm["BestFitParameters"]

The parameters are believed to be around :

alpha = 0.1
beta= 0.1
kappa = 0.05
a0 = 0.01
b0 = 3

But it didn't work... :

NonlinearModelFit::nrnum: The function value 1/2 ((-22.6124+R$3721[3.,22.])^2+(-119.51+R$3721[3.,119.])^2+(-24.738+R$3721[6.,22.])^2+(-60.1536+R$3721[6.,60.])^2+(-126.123+R$3721[6.,119.])^2+(-16.8895+R$3721[9.,17.])^2+(-25.4959+R$3721[9.,22.])^2+(-57.9807+R$3721[9.,60.])^2+(-110.446+R$3721[9.,119.])^2+(-17.3404+R$3721[12.,17.])^2+(-26.1946+R$3721[12.,22.])^2+(-60.9089+R$3721[12.,60.])^2+(-110.332+R$3721[12.,119.])^2+<<25>>+(-200.187+R$3721[27.,185.])^2+(-20.6519+R$3721[30.,17.])^2+(-34.5678+R$3721[30.,22.])^2+(-68.705+R$3721[30.,60.])^2+(-111.198+R$3721[30.,119.])^2+(-199.25+R$3721[30.,185.])^2+(-19.4591+R$3721[33.,17.])^2+(-35.9263+R$3721[33.,22.])^2+(-68.2107+R$3721[33.,60.])^2+(-109.903+R$3721[33.,119.])^2+(-198.411+R$3721[33.,185.])^2+(-20.6855+R$3721[36.,17.])^2+<<819>>) is not a real number at {alpha,beta,kappa,a0,b0} = {0.1,0.1,0.05,0.01,3.}.

I assume there is an issue with $R0$, but I don't get where.

How could I proceed ?

Also, I don't know how I could fix a priori the initial conditions for each fit in order to extract only the shared parameters.

DATA

MathematicaStackExchange doesn't give the possibility to enter to much characters. I can give only the joined data.

1. joined data with R0 as a variable

Here is the joined data.

data={{9., 17., 16.8895}, {12., 17., 17.3404}, {15., 17., 17.1633}, {18., 
  17., 19.3417}, {21., 17., 17.9899}, {24., 17., 19.9677}, {27., 17., 
  19.4362}, {30., 17., 20.6519}, {33., 17., 19.4591}, {36., 17., 
  20.6855}, {39., 17., 20.1952}, {42., 17., 21.9949}, {45., 17., 
  21.0234}, {48., 17., 22.7408}, {51., 17., 22.3908}, {54., 17., 
  25.0918}, {57., 17., 23.5989}, {60., 17., 26.0703}, {63., 17., 
  24.5605}, {66., 17., 27.2539}, {69., 17., 26.1619}, {72., 17., 
  28.4762}, {75., 17., 27.5854}, {78., 17., 29.8393}, {81., 17., 
  28.3553}, {84., 17., 30.3221}, {87., 17., 29.675}, {90., 17., 
  31.5653}, {93., 17., 30.5337}, {96., 17., 33.3734}, {99., 17., 
  31.6876}, {102., 17., 34.1503}, {105., 17., 33.3065}, {108., 17., 
  35.3291}, {111., 17., 33.9209}, {114., 17., 36.773}, {117., 17., 
  35.4094}, {120., 17., 41.5902}, {123., 17., 36.1305}, {126., 17., 
  37.971}, {129., 17., 36.402}, {132., 17., 39.1158}, {135., 17., 
  38.0177}, {138., 17., 40.8558}, {141., 17., 39.6065}, {144., 17., 
  40.9749}, {147., 17., 39.8896}, {150., 17., 41.8237}, {153., 17., 
  40.5802}, {156., 17., 42.3858}, {159., 17., 40.6619}, {162., 17., 
  44.4442}, {165., 17., 45.4162}, {168., 17., 46.1884}, {171., 17., 
  44.6008}, {174., 17., 47.1647}, {177., 17., 45.3808}, {180., 17., 
  46.5859}, {183., 17., 45.3035}, {186., 17., 47.6604}, {189., 17., 
  46.6771}, {192., 17., 45.9242}, {195., 17., 46.767}, {198., 17., 
  44.6899}, {201., 17., 46.6628}, {204., 17., 46.1571}, {207., 17., 
  46.5555}, {210., 17., 44.835}, {213., 17., 45.1423}, {216., 17., 
  45.1954}, {219., 17., 45.309}, {222., 17., 47.7791}, {225., 17., 
  46.7777}, {228., 17., 48.135}, {231., 17., 45.6493}, {234., 17., 
  45.8933}, {237., 17., 46.1803}, {240., 17., 46.7285}, {243., 17., 
  46.8063}, {246., 17., 47.1679}, {249., 17., 46.8787}, {252., 17., 
  47.2715}, {255., 17., 47.5362}, {258., 17., 48.9234}, {261., 17., 
  47.5456}, {264., 17., 53.5554}, {267., 17., 52.5704}, {270., 17., 
  49.6049}, {273., 17., 49.1189}, {276., 17., 48.9498}, {279., 17., 
  49.6024}, {282., 17., 49.7491}, {285., 17., 53.1681}, {288., 17., 
  51.7124}, {291., 17., 50.8069}, {294., 17., 50.0237}, {297., 17., 
  50.5922}, {300., 17., 50.6518}, {303., 17., 50.8827}, {306., 17., 
  51.2245}, {309., 17., 51.0911}, {312., 17., 52.3379}, {315., 17., 
  52.5112}, {318., 17., 53.9182}, {321., 17., 53.7082}, {324., 17., 
  54.9239}, {327., 17., 53.7369}, {330., 17., 51.7204}, {333., 17., 
  55.993}, {336., 17., 56.8489}, {339., 17., 53.3037}, {342., 17., 
  52.0201}, {345., 17., 52.6267}, {348., 17., 52.5615}, {351., 17., 
  55.4133}, {354., 17., 55.5549}, {357., 17., 52.2672}, {360., 17., 
  54.2202}, {363., 17., 50.3245}, {366., 17., 54.0435}, {369., 17., 
  51.0724}, {372., 17., 51.2091}, {375., 17., 51.6602}, {378., 17., 
  51.3684}, {381., 17., 51.5346}, {384., 17., 51.9204}, {387., 17., 
  52.3952}, {390., 17., 52.9114}, {393., 17., 54.3833}, {396., 17., 
  55.1898}, {399., 17., 51.3853}, {402., 17., 55.048}, {405., 17., 
  50.8574}, {408., 17., 51.9619}, {411., 17., 52.5775}, {414., 17., 
  52.5676}, {417., 17., 51.0891}, {420., 17., 54.3895}, {423., 17., 
  54.7591}, {426., 17., 53.9934}, {429., 17., 53.8877}, {435., 17., 
  55.4067}, {441., 17., 56.0656}, {447., 17., 57.4607}, {453., 17., 
  51.6628}, {456., 17., 54.3568}, {459., 17., 57.6827}, {465., 17., 
  54.8474}, {468., 17., 51.0797}, {471., 17., 53.1862}, {474., 17., 
  53.3921}, {477., 17., 54.468}, {480., 17., 54.1083}, {483., 17., 
  50.7948}, {486., 17., 53.3431}, {489., 17., 48.8646}, {492., 17., 
  53.3906}, {495., 17., 51.6016}, {498., 17., 54.1742}, {501., 17., 
  54.6549}, {504., 17., 50.0598}, {507., 17., 53.849}, {510., 17., 
  52.6431}, {513., 17., 54.3103}, {516., 17., 50.5004}, {519., 17., 
  50.8213}, {522., 17., 50.8512}, {525., 17., 52.4319}, {528., 17., 
  55.2716}, {3., 22., 22.6124}, {6., 22., 24.738}, {9., 22., 
  25.4959}, {12., 22., 26.1946}, {15., 22., 27.6091}, {18., 22., 
  29.1024}, {21., 22., 30.6462}, {24., 22., 32.9126}, {27., 22., 
  34.1471}, {30., 22., 34.5678}, {33., 22., 35.9263}, {36., 22., 
  37.4284}, {39., 22., 38.5027}, {42., 22., 39.5611}, {45., 22., 
  40.743}, {48., 22., 41.9482}, {51., 22., 42.7558}, {54., 22., 
  43.5064}, {57., 22., 44.43}, {60., 22., 45.7449}, {63., 22., 
  47.0524}, {66., 22., 48.0848}, {69., 22., 48.8836}, {72., 22., 
  49.6807}, {75., 22., 50.6801}, {78., 22., 51.6959}, {81., 22., 
  52.6475}, {84., 22., 53.5902}, {87., 22., 54.4008}, {90., 22., 
  54.774}, {93., 22., 55.6085}, {96., 22., 56.3299}, {99., 22., 
  56.4428}, {102., 22., 56.7936}, {105., 22., 57.4926}, {108., 22., 
  58.2406}, {111., 22., 59.1169}, {114., 22., 59.5766}, {117., 22., 
  59.7909}, {120., 22., 61.6917}, {123., 22., 62.4342}, {126., 22., 
  61.5979}, {129., 22., 61.8203}, {132., 22., 62.5629}, {135., 22., 
  63.4556}, {138., 22., 63.688}, {141., 22., 63.9159}, {144., 22., 
  63.9802}, {147., 22., 64.1833}, {150., 22., 64.3304}, {153., 22., 
  64.3847}, {156., 22., 64.6173}, {159., 22., 64.9009}, {162., 22., 
  65.1622}, {165., 22., 65.4684}, {168., 22., 65.5182}, {171., 22., 
  66.1171}, {174., 22., 66.4103}, {177., 22., 66.2592}, {180., 22., 
  66.185}, {183., 22., 65.8147}, {186., 22., 65.733}, {189., 22., 
  65.6618}, {192., 22., 64.7882}, {195., 22., 64.8274}, {198., 22., 
  64.9444}, {201., 22., 63.1305}, {204., 22., 62.3995}, {207., 22., 
  63.0431}, {210., 22., 62.2181}, {213., 22., 62.5286}, {216., 22., 
  62.1711}, {219., 22., 60.8353}, {222., 22., 60.7586}, {225., 22., 
  60.7004}, {228., 22., 59.5638}, {231., 22., 59.1517}, {234., 22., 
  58.9346}, {237., 22., 59.0493}, {240., 22., 59.5229}, {243., 22., 
  58.0876}, {246., 22., 56.247}, {249., 22., 56.173}, {252., 22., 
  56.1419}, {255., 22., 55.2417}, {258., 22., 56.2456}, {261., 22., 
  57.9169}, {264., 22., 60.728}, {267., 22., 63.6912}, {270., 22., 
  61.9647}, {273., 22., 57.0852}, {276., 22., 54.2803}, {279., 22., 
  55.3487}, {282., 22., 58.0208}, {285., 22., 60.8749}, {288., 22., 
  61.029}, {291., 22., 59.3053}, {294., 22., 56.7078}, {297., 22., 
  53.8873}, {300., 22., 55.2545}, {303., 22., 56.5482}, {306., 22., 
  56.0664}, {309., 22., 55.2537}, {312., 22., 55.3196}, {315., 22., 
  55.8909}, {318., 22., 55.6318}, {321., 22., 56.213}, {324., 22., 
  55.4207}, {327., 22., 54.2877}, {330., 22., 55.1178}, {333., 22., 
  51.193}, {336., 22., 48.5713}, {339., 22., 49.5028}, {342., 22., 
  49.4166}, {345., 22., 50.0304}, {348., 22., 50.9326}, {351., 22., 
  52.014}, {354., 22., 50.2956}, {357., 22., 49.8529}, {360., 22., 
  50.8205}, {363., 22., 51.376}, {366., 22., 50.6679}, {369., 22., 
  51.6815}, {372., 22., 53.5813}, {375., 22., 53.7359}, {378., 22., 
  54.6252}, {381., 22., 55.2786}, {384., 22., 53.4308}, {387., 22., 
  54.5401}, {390., 22., 57.9795}, {393., 22., 55.2026}, {396., 22., 
  55.386}, {399., 22., 59.8766}, {402., 22., 58.1028}, {405., 22., 
  57.129}, {408., 22., 56.9853}, {411., 22., 57.2221}, {414., 22., 
  56.9648}, {417., 22., 55.586}, {420., 22., 56.7903}, {423., 22., 
  56.2825}, {426., 22., 53.8012}, {429., 22., 52.6652}, {432., 22., 
  54.2455}, {435., 22., 56.3002}, {438., 22., 56.2343}, {441., 22., 
  56.7575}, {444., 22., 56.7977}, {447., 22., 56.3049}, {450., 22., 
  54.6538}, {453., 22., 52.5136}, {456., 22., 52.3433}, {459., 22., 
  52.828}, {462., 22., 54.0433}, {465., 22., 51.5131}, {468., 22., 
  50.4781}, {471., 22., 52.6831}, {474., 22., 52.4475}, {477., 22., 
  52.6825}, {480., 22., 52.5579}, {483., 22., 52.8213}, {486., 22., 
  53.6997}, {489., 22., 53.3714}, {492., 22., 52.3218}, {495., 22., 
  52.3176}, {498., 22., 53.8036}, {501., 22., 53.7502}, {504., 22., 
  55.6969}, {507., 22., 56.1864}, {510., 22., 52.9824}, {513., 22., 
  55.2477}, {516., 22., 54.727}, {519., 22., 54.0447}, {522., 22., 
  56.1034}, {525., 22., 53.0694}, {528., 22., 51.3001}, {6., 60., 
  60.1536}, {9., 60., 57.9807}, {12., 60., 60.9089}, {15., 60., 
  59.4291}, {18., 60., 61.3227}, {21., 60., 61.8788}, {24., 60., 
  67.2192}, {27., 60., 66.2767}, {30., 60., 68.705}, {33., 60., 
  68.2107}, {36., 60., 70.8731}, {39., 60., 68.7269}, {42., 60., 
  73.2306}, {45., 60., 72.3068}, {48., 60., 74.8006}, {51., 60., 
  72.1975}, {54., 60., 76.577}, {57., 60., 75.5894}, {60., 60., 
  76.342}, {63., 60., 75.5134}, {66., 60., 77.47}, {69., 60., 
  76.6854}, {72., 60., 78.7422}, {75., 60., 78.6074}, {78., 60., 
  81.0158}, {81., 60., 82.8521}, {84., 60., 85.1395}, {87., 60., 
  85.211}, {90., 60., 84.5157}, {93., 60., 83.622}, {96., 60., 
  88.1703}, {99., 60., 85.6195}, {102., 60., 86.8345}, {105., 60., 
  86.5568}, {108., 60., 87.5942}, {111., 60., 88.3053}, {114., 60., 
  88.3475}, {117., 60., 89.3993}, {120., 60., 91.7091}, {123., 60., 
  89.7268}, {126., 60., 90.6704}, {129., 60., 89.7999}, {132., 60., 
  90.369}, {135., 60., 88.7787}, {138., 60., 90.3022}, {141., 60., 
  89.8267}, {144., 60., 91.2241}, {147., 60., 91.2859}, {150., 60., 
  92.992}, {153., 60., 91.0079}, {156., 60., 93.0784}, {159., 60., 
  90.8868}, {162., 60., 92.7426}, {165., 60., 92.757}, {168., 60., 
  94.4202}, {171., 60., 92.2914}, {174., 60., 90.3876}, {177., 60., 
  89.3376}, {180., 60., 89.814}, {183., 60., 88.9134}, {186., 60., 
  89.7058}, {189., 60., 91.642}, {192., 60., 90.3205}, {195., 60., 
  87.8566}, {198., 60., 87.6065}, {201., 60., 87.0403}, {204., 60., 
  87.3344}, {207., 60., 87.2313}, {210., 60., 87.3705}, {213., 60., 
  86.9135}, {216., 60., 87.2684}, {219., 60., 87.2989}, {222., 60., 
  85.4766}, {225., 60., 85.3534}, {228., 60., 86.535}, {231., 60., 
  86.1929}, {234., 60., 86.089}, {237., 60., 85.9466}, {240., 60., 
  85.1389}, {243., 60., 85.0242}, {246., 60., 84.4313}, {249., 60., 
  83.7604}, {252., 60., 81.9419}, {255., 60., 83.773}, {258., 60., 
  82.7046}, {261., 60., 84.7331}, {264., 60., 86.0393}, {267., 60., 
  84.7472}, {270., 60., 79.1677}, {273., 60., 80.9426}, {276., 60., 
  79.9624}, {279., 60., 75.5272}, {282., 60., 79.3103}, {285., 60., 
  80.8015}, {288., 60., 81.3927}, {291., 60., 80.1678}, {294., 60., 
  80.268}, {297., 60., 79.9067}, {300., 60., 76.9766}, {303., 60., 
  81.8132}, {306., 60., 73.6449}, {309., 60., 76.4059}, {312., 60., 
  76.4056}, {315., 60., 81.7311}, {318., 60., 80.8468}, {321., 60., 
  80.958}, {324., 60., 86.9248}, {327., 60., 78.3434}, {330., 60., 
  74.8752}, {333., 60., 78.0912}, {336., 60., 81.5165}, {339., 60., 
  72.7919}, {342., 60., 74.2966}, {345., 60., 79.2233}, {348., 60., 
  81.9791}, {351., 60., 74.3276}, {354., 60., 85.1221}, {357., 60., 
  78.6944}, {360., 60., 75.8183}, {363., 60., 75.6696}, {366., 60., 
  75.9147}, {369., 60., 76.3326}, {372., 60., 80.0048}, {375., 60., 
  79.8311}, {378., 60., 79.0427}, {381., 60., 81.8084}, {384., 60., 
  73.5742}, {387., 60., 84.2291}, {390., 60., 84.9122}, {393., 60., 
  82.6657}, {396., 60., 78.2888}, {399., 60., 90.0235}, {402., 60., 
  83.3667}, {405., 60., 81.7737}, {408., 60., 81.19}, {411., 60., 
  82.3131}, {414., 60., 79.8072}, {417., 60., 74.4822}, {420., 60., 
  75.6291}, {423., 60., 82.2655}, {426., 60., 73.704}, {429., 60., 
  81.4184}, {432., 60., 72.1127}, {435., 60., 74.7053}, {438., 60., 
  79.4664}, {441., 60., 86.4491}, {444., 60., 79.5096}, {447., 60., 
  77.1761}, {450., 60., 83.082}, {453., 60., 80.3418}, {456., 60., 
  85.3873}, {459., 60., 85.7409}, {462., 60., 73.3735}, {465., 60., 
  72.2276}, {468., 60., 82.7752}, {471., 60., 71.6917}, {474., 60., 
  78.5233}, {477., 60., 82.4042}, {480., 60., 83.8073}, {483., 60., 
  91.5845}, {486., 60., 82.8906}, {489., 60., 87.3935}, {492., 60., 
  89.9856}, {495., 60., 74.1819}, {498., 60., 77.5752}, {501., 60., 
  82.6796}, {504., 60., 79.2659}, {507., 60., 81.5865}, {510., 60., 
  82.709}, {513., 60., 88.4083}, {516., 60., 81.7317}, {519., 60., 
  76.2638}, {522., 60., 86.2863}, {525., 60., 93.2163}, {528., 60., 
  82.6943}, {3., 119., 119.51}, {6., 119., 126.123}, {9., 119., 
  110.446}, {12., 119., 110.332}, {15., 119., 110.478}, {18., 119., 
  111.335}, {21., 119., 109.536}, {24., 119., 111.901}, {27., 119., 
  110.46}, {30., 119., 111.198}, {33., 119., 109.903}, {36., 119., 
  110.72}, {39., 119., 109.635}, {42., 119., 110.643}, {45., 119., 
  109.528}, {48., 119., 110.348}, {51., 119., 110.117}, {54., 119., 
  109.117}, {57., 119., 108.536}, {60., 119., 108.615}, {63., 119., 
  109.495}, {66., 119., 111.304}, {69., 119., 111.139}, {72., 119., 
  114.285}, {75., 119., 113.627}, {78., 119., 114.77}, {81., 119., 
  114.544}, {84., 119., 115.304}, {87., 119., 114.895}, {90., 119., 
  115.859}, {93., 119., 114.357}, {96., 119., 115.038}, {99., 119., 
  114.305}, {102., 119., 115.09}, {105., 119., 114.815}, {108., 119., 
  113.203}, {111., 119., 113.46}, {114., 119., 114.573}, {117., 119., 
  113.339}, {120., 119., 114.354}, {123., 119., 112.285}, {126., 119.,
   112.695}, {129., 119., 112.032}, {132., 119., 112.253}, {135., 
  119., 108.945}, {138., 119., 109.271}, {141., 119., 108.654}, {144.,
   119., 104.336}, {147., 119., 103.609}, {150., 119., 
  105.778}, {153., 119., 105.077}, {156., 119., 104.868}, {159., 119.,
   103.945}, {162., 119., 104.039}, {165., 119., 101.727}, {168., 
  119., 97.6562}, {171., 119., 99.6703}, {174., 119., 96.6503}, {177.,
   119., 98.3032}, {180., 119., 98.8859}, {183., 119., 
  97.9825}, {186., 119., 94.8383}, {189., 119., 93.4101}, {192., 119.,
   88.9132}, {195., 119., 91.7409}, {198., 119., 93.2425}, {201., 
  119., 86.1268}, {204., 119., 84.9263}, {207., 119., 86.3445}, {210.,
   119., 84.4667}, {213., 119., 85.9353}, {216., 119., 
  85.7998}, {219., 119., 85.2672}, {222., 119., 86.3356}, {225., 119.,
   86.7423}, {228., 119., 86.1353}, {231., 119., 84.8631}, {234., 
  119., 84.7305}, {237., 119., 83.385}, {240., 119., 87.5174}, {243., 
  119., 83.3014}, {246., 119., 86.9219}, {249., 119., 78.3219}, {252.,
   119., 78.9197}, {255., 119., 74.785}, {258., 119., 67.8261}, {261.,
   119., 75.8036}, {264., 119., 86.2339}, {267., 119., 
  87.3689}, {270., 119., 88.1322}, {273., 119., 86.1332}, {276., 119.,
   89.9111}, {279., 119., 90.5619}, {282., 119., 88.4012}, {285., 
  119., 85.5809}, {288., 119., 76.692}, {291., 119., 80.0753}, {294., 
  119., 90.1118}, {297., 119., 91.8565}, {300., 119., 85.0882}, {303.,
   119., 89.1269}, {306., 119., 96.8869}, {309., 119., 
  75.4618}, {312., 119., 96.3013}, {315., 119., 89.4435}, {318., 119.,
   103.21}, {321., 119., 94.6233}, {324., 119., 102.48}, {327., 119., 
  96.7664}, {330., 119., 84.2408}, {333., 119., 97.3822}, {336., 119.,
   74.2619}, {339., 119., 87.2886}, {342., 119., 118.024}, {345., 
  119., 113.648}, {348., 119., 112.4}, {351., 119., 107.295}, {354., 
  119., 111.618}, {357., 119., 112.181}, {360., 119., 112.119}, {363.,
   119., 90.6252}, {366., 119., 106.837}, {369., 119., 
  99.7227}, {372., 119., 97.5255}, {375., 119., 108.211}, {378., 119.,
   117.211}, {381., 119., 97.9301}, {384., 119., 104.567}, {387., 
  119., 117.343}, {390., 119., 121.622}, {393., 119., 106.117}, {396.,
   119., 116.022}, {399., 119., 118.856}, {402., 119., 
  106.854}, {405., 119., 112.418}, {408., 119., 112.79}, {411., 119., 
  112.225}, {414., 119., 116.686}, {417., 119., 111.297}, {420., 119.,
   115.404}, {423., 119., 117.563}, {426., 119., 116.243}, {429., 
  119., 119.805}, {432., 119., 112.863}, {435., 119., 103.505}, {438.,
   119., 116.846}, {441., 119., 115.508}, {444., 119., 
  115.579}, {447., 119., 101.756}, {450., 119., 102.848}, {453., 119.,
   112.506}, {456., 119., 113.93}, {459., 119., 116.386}, {462., 119.,
   108.138}, {465., 119., 108.635}, {468., 119., 110.514}, {471., 
  119., 108.217}, {474., 119., 110.008}, {477., 119., 95.7788}, {480.,
   119., 92.8073}, {483., 119., 104.382}, {486., 119., 98.77}, {489., 
  119., 112.527}, {492., 119., 94.6092}, {495., 119., 89.2861}, {498.,
   119., 92.0002}, {501., 119., 98.7618}, {504., 119., 
  105.274}, {507., 119., 96.7057}, {510., 119., 93.5207}, {513., 119.,
   90.5992}, {516., 119., 87.1486}, {519., 119., 103.466}, {522., 
  119., 100.133}, {525., 119., 120.605}, {528., 119., 125.717}, {12., 
  185., 185.791}, {15., 185., 199.035}, {18., 185., 197.796}, {21., 
  185., 185.256}, {24., 185., 199.576}, {27., 185., 200.187}, {30., 
  185., 199.25}, {33., 185., 198.411}, {36., 185., 198.288}, {39., 
  185., 194.506}, {42., 185., 189.658}, {45., 185., 191.203}, {48., 
  185., 185.757}, {51., 185., 183.642}, {54., 185., 183.513}, {57., 
  185., 186.524}, {60., 185., 182.793}, {63., 185., 182.218}, {66., 
  185., 182.045}, {69., 185., 176.614}, {72., 185., 182.432}, {75., 
  185., 181.409}, {78., 185., 182.438}, {81., 185., 179.939}, {84., 
  185., 182.435}, {87., 185., 181.521}, {90., 185., 176.654}, {93., 
  185., 175.39}, {96., 185., 179.446}, {99., 185., 173.541}, {102., 
  185., 176.645}, {105., 185., 176.715}, {108., 185., 173.915}, {111.,
   185., 173.14}, {114., 185., 173.045}, {117., 185., 160.089}, {120.,
   185., 165.306}, {123., 185., 165.906}, {126., 185., 
  165.712}, {129., 185., 159.285}, {132., 185., 163.219}, {135., 185.,
   156.287}, {138., 185., 150.445}, {141., 185., 153.388}, {144., 
  185., 138.083}, {147., 185., 137.152}, {150., 185., 133.003}, {153.,
   185., 130.634}, {156., 185., 131.832}, {159., 185., 
  136.142}, {162., 185., 133.906}, {165., 185., 130.929}, {168., 185.,
   136.717}, {171., 185., 129.749}, {174., 185., 148.377}, {177., 
  185., 133.068}, {180., 185., 149.921}, {183., 185., 134.802}, {186.,
   185., 150.543}, {189., 185., 138.678}, {192., 185., 147.06}, {195.,
   185., 143.604}, {198., 185., 143.368}, {201., 185., 
  140.587}, {204., 185., 138.171}, {207., 185., 140.699}, {210., 185.,
   137.346}, {213., 185., 126.241}, {216., 185., 131.743}, {219., 
  185., 134.835}, {222., 185., 134.086}, {225., 185., 137.185}, {228.,
   185., 135.892}, {231., 185., 141.62}, {234., 185., 135.963}, {237.,
   185., 133.382}, {240., 185., 134.258}, {243., 185., 
  141.568}, {246., 185., 137.642}, {249., 185., 131.681}, {252., 185.,
   132.635}, {255., 185., 134.506}, {258., 185., 136.089}, {261., 
  185., 138.973}, {264., 185., 141.048}, {267., 185., 133.785}, {270.,
   185., 133.245}, {273., 185., 116.408}, {276., 185., 123.9}, {279., 
  185., 120.251}, {282., 185., 116.984}, {285., 185., 135.753}, {288.,
   185., 123.026}, {291., 185., 112.116}, {294., 185., 
  134.164}, {297., 185., 134.548}, {300., 185., 129.032}, {303., 185.,
   116.97}, {306., 185., 113.993}, {309., 185., 99.4695}, {312., 185.,
   97.4854}, {315., 185., 100.422}, {318., 185., 117.461}, {321., 
  185., 99.4758}, {324., 185., 106.366}, {327., 185., 108.271}, {330.,
   185., 104.738}, {333., 185., 117.487}, {336., 185., 
  101.704}, {339., 185., 101.32}, {342., 185., 112.97}, {345., 185., 
  96.6092}, {348., 185., 99.2531}, {351., 185., 120.19}, {354., 185., 
  124.284}, {357., 185., 130.082}, {360., 185., 121.699}, {363., 185.,
   108.539}, {366., 185., 103.98}, {369., 185., 100.293}, {372., 185.,
   94.7848}, {375., 185., 103.281}, {378., 185., 114.4}, {381., 185., 
  94.8752}, {384., 185., 101.51}, {387., 185., 104.285}, {390., 185., 
  107.424}, {393., 185., 112.506}, {396., 185., 104.061}, {399., 185.,
   113.713}, {402., 185., 136.378}, {405., 185., 134.92}, {408., 185.,
   139.111}, {411., 185., 143.397}, {414., 185., 139.998}, {417., 
  185., 137.19}, {420., 185., 143.812}, {423., 185., 133.346}, {426., 
  185., 141.8}, {429., 185., 136.171}, {432., 185., 137.842}, {435., 
  185., 147.509}, {438., 185., 140.488}, {441., 185., 142.855}, {444.,
   185., 151.992}, {447., 185., 145.348}, {450., 185., 
  138.757}, {453., 185., 135.964}, {456., 185., 140.381}, {459., 185.,
   143.697}, {462., 185., 136.854}, {465., 185., 129.477}, {468., 
  185., 138.181}, {471., 185., 142.726}, {474., 185., 143.633}, {477.,
   185., 133.913}, {480., 185., 157.635}, {483., 185., 
  147.941}, {486., 185., 142.015}, {489., 185., 130.545}, {492., 185.,
   141.941}, {495., 185., 142.863}, {498., 185., 135.462}, {501., 
  185., 139.637}, {504., 185., 128.002}, {507., 185., 140.211}, {510.,
   185., 140.209}, {513., 185., 132.36}, {516., 185., 141.088}, {519.,
   185., 142.756}, {522., 185., 152.256}, {525., 185., 
  164.725}, {528., 185., 153.737}}
Improve this question
$\endgroup$
13
  • 1
    $\begingroup$ Observe the output of the following code piece and think about the reason: expr = x; Module[{x}, {x, expr}] $\endgroup$
    – xzczd ♦
    Commented May 26, 2019 at 11:36
  • 1
    $\begingroup$ No, this isn't that related to evaluation order. (Though it's not impossible to fix by tackling evaluation order. ) The main issue here is, 1. Module localizes variable by renaming it as …$…. (In the case above, x becomes x$24532); 2. When localizing, Module only sees explicit variable i.e. the x inside expr is not noticed by Module so it's not localized and not renamed. And the same thing happens on your R. $\endgroup$
    – xzczd ♦
    Commented May 26, 2019 at 12:39
  • 1
    $\begingroup$ Only those explicit Rs are localized. You have Rs inside eq. $\endgroup$
    – xzczd ♦
    Commented May 26, 2019 at 12:48
  • 1
    $\begingroup$ Please show the specific code. Also, the {R0, 0, 300} in NDSolve only makes the code unnecessarily complicated, and you've missed the initial condition for chi in your code, and D[chi[t, R0], t][t, R0] is obviously wrong. As to the update, alpha >= 0, beta >= 0, kappa >= 0, a0 >= 0, b0 >= 0 doesn't make sense, please don't guess the syntax, check the document of NDSolve carefully. $\endgroup$
    – xzczd ♦
    Commented May 26, 2019 at 13:17
  • 1
    $\begingroup$ @J.A You can (and arguably should) use pastebin.com. $\endgroup$
    – AccidentalFourierTransform
    Commented May 26, 2019 at 22:18

1 Answer 1

Reset to default
7
+50
$\begingroup$

The subject of parameter fitting comes up frequently on MSE. Parameter fitting is a difficult subject and will depend on your data quality, your model, and your intial guesses. I have been dabbling with StringTemplates as a potential way to encapsulate some of the basic parameter fitting work flow.

Approach

  • Use ParametricNDSolveValue to create the model.
  • Use StringTemplates to handle lists of parameters and variables.
  • Generate a Manipulate slider model to debug model and understand the effects of parameter changes.
  • Transfer initial guesses from manipulate to perform a fit.

Implementation

I commented the code so I hope it self explanatory. First assign the constants and prep the data.

(* Evaluate data first *)
(* Constants *)
l = 10^(-5);
k = 1/l;
chic = 0.5;
T = 550;
(* Get unique R0s *)
R0s = Union@data[[All, 2]];
(* Subset Matching R0 and Delete 2nd Column *)
rdat = (Cases[data, {_, #, _}][[All, {1, 3}]] & /@ R0s);

Now, set up the equations and the Manipulate slider to view how the model behaves and try to improve initial parameters estimates.

(* Generate System of Differential Equations *)
e1 = R'[t] == -a[t]*R[t] + b[t];
e3 = b'[t] == beta/2*(Tanh[(chi[t] - chic)*k] - 1);
e2 = a'[t] == -alpha/2*(Tanh[(chi[t] - chic)*k] - 1);
e4 = chi'[t] == -kappa*R[t]*(chi[t] - 2*chic);
ics = {a[0] == a0, b[0] == b0, R[0] == R0, chi[0] == 0};
eqns = {e1, e2, e3, e4}~Join~ics;
(*Variables*)
vbles = {R, a, b, chi};
(*Parameters with target and desired ranges*)
mat = {
   {alpha, 0.1, 0.00025, 0.5},
   {beta, 0.1, 0.00025, 0.5},
   {kappa, 0.05, 0.0125, 0.1},
   {a0, 0.01, 0.00005, 0.1},
   {b0, 3, 1, 6},
   {R0, 17, 17, 185}
   };
(* reduce the matrix because R0 does not participate in parameter \
fits *)
rmat = mat[[1 ;; -2]];
(* Build Manipulate sliders *)
sfun =  StringRiffle[(StringTemplate[
         "{{`1`,`2`},`3`,`4`,Appearance\[Rule]\"Labeled\"}"] @@ #) & \
/@ #, ","] &;
sliders = sfun[rmat];
(* Extract Parameters from mat *)
parms = mat[[All, 1]];
rparms = rmat[[All, 1]];
(* Create String Representations of parms *)
sparms = StringRiffle[ToString[#] & /@ parms, ","];
rsparms = StringRiffle[ToString[#] & /@ rparms, ","];
(* Create patterns and string reps of parameters *)
pats = Pattern @@@ (#*_ & /@ parms);
spats = StringRiffle[ToString[#] & /@ pats, ","];
(* List Plot of the data *)
lp = Graphics[{Hue[#2/185], PointSize[0.01], Point[{#1, #3}]} & @@@ 
    data, Axes -> True];
(* ParametricNDSolveValue *)
pfun = ParametricNDSolveValue[eqns, vbles, {t, 0, T}, parms];
(*Create an appropriate model function to fit*)
modelstring = "(#[[1]])&";
(* Create some PlotLegends *)
pl = ",PlotLegends\[Rule]{" <> 
   StringRiffle["\"R0=" <> ToString[#] <> "\"" & /@ R0s, ","] <> "}";
(* Build the model expression *)
ToExpression[
 StringTemplate[
   "model[`pats`][t_]:=`ms`@Through[pfun[`params`][t],List]\
/;And@@NumericQ/@{`params`};"][<|"pats" -> spats, "params" -> sparms, 
   "ms" -> modelstring|>]]
(* Create slider model *)
globalstring = 
  StringTemplate["global={`params`};"][<|"params" -> rsparms|>];
mantemp = 
  "Manipulate[`g`\[IndentingNewLine]Show[lp,Plot[Evaluate@({model[\
alpha,beta,kappa,a0,b0,#][t]}&/@R0s),{t,0,T},PlotRange\[Rule]{0,200}`\
pl`],ImageSize->Large],`sliders`]";
ToExpression@
 StringTemplate[mantemp][<|"sliders" -> sliders, "params" -> rsparms, 
   "pl" -> pl, "g" -> globalstring|>]
(*Display global variable*)
Dynamic@global

Manipulate Image

Now set up to fit the funtions for each R0 value.

(* Grab The initial parameter guesses *)
initguess = MapThread[List, {rparms, First@Dynamic@global}];
(* Create a fit function to operate on different R0s *)
fitfn = FindFit[rdat[[#]], 
    model[alpha, beta, kappa, a0, b0, R0s[[#]]][t], initguess, t, 
    Method -> "Gradient"] &;
(* Perform Fits on R0s *)
fits = fitfn[#][[All, 2]] & /@ Range@Length@R0s;
(* Display Results *)
fits // MatrixForm
Mean@fits

Parameter Estimates

The data is noisy leading to some dodgy results for the high R0. You can experiment with different fitting options, but you may need to improve your model and/or your data acquisition.

Update to Fit per $R_0$ Data set

As requested, here is a way to fit per data set. I also allowed $R_0$ to be fit using the column value as an initial guess. In this case, each fitted row is plotted. A word of caution, some fitting methods will run forever, so you may need to experiment.

(* Grab The initial parameter guesses from dynamic variable of slider \
*)
initguess = 
  MapThread[List, {parms, (First@Dynamic@global)~Join~{R0s[[#]]}}] &;
(* Create a fit function to operate on different R0s *)
fitfn = FindFit[rdat[[#]], model[alpha, beta, kappa, a0, b0, R0][t], 
    initguess[#], t, Method -> "Gradient", WorkingPrecision -> 10] &;
(* Perform Fits on R0s *)
(*fits = fitfn[#][[All,2]]&/@Range@Length@R0s;*)
fits = fitfn[#][[All, 2]] & /@ {1, 2, 3, 4, 5};
(* Display Results *)
fits // MatrixForm
mfit = Mean@fits
mat2 = rmat;
mat2[[All, 2]] = mfit[[1 ;; -2]];
Show[{lp, 
  Plot[Evaluate@((model @@ #)[t] & /@ fits), {t, 0, T}, 
   PlotRange -> {0, 200}, 
   PlotLegends -> {"R0=17.", "R0=22.", "R0=60.", "R0=119.", 
     "R0=185."}]}, ImageSize -> Large]

Individual Fits

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks this is a great job, and a very nice way to proceed. The only question I have is with regard to the fitting process. Is taking the mean the only way to proceed ? Meaning is there another way to gather the informations from the different data-sets ? $\endgroup$
    – J.A
    Commented May 27, 2019 at 15:43
  • $\begingroup$ I added a plot of the fitted parameters for each row. $\endgroup$
    – Tim Laska
    Commented May 27, 2019 at 21:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.