Off[General::spell]
Off[General::spell1]
Off[FindRoot::cvnwt]
Off[FindRoot::frmp]
<< Statistics`

\(\(\(Clear[k, \ P, \ ode1, ode2, a, b, c, \ CountSS, \ LSS, \ j, \ 
      i\ , \ q, \ SolutionTensor, \ ParamTensor, \ 
      PlotResolution]\)\(
\)\(
\)
    \) (*\ Parameter\ Values\ *) \)
a = 0.5;
b =  .25;
c = 0.75;

(*Define Parameters  -  
        replace these with whatever parameter names you want*)

(*Define parameters and distributions .   List is given by parameter name,  
      PDFTYPE,  {min, 
            max} .   The sampling is for uniform distributions .   You can add more parameters by just adding more elements to this list.*) 
ParamTensor = 
{
{a,  
            1,  { .3,  .3},  {}},  
{b,  
            1,  { .01, 10.0},  {}},  
{c,  
            1,  {0, 
              8},  {}},  
};

NN = 
        5;
paramnumber = 
        3;   
solvetime = 
        20;

(*Generate the lists of sampled parameter values *)

\(
 (*Divide\ Parameter\ Range\ into\ Equiprobable\ Intervals
Assume\ UNIFORM\ distribution*) 
\(For[j = 1, \ 
        j ≤ 
          paramnumber, \ \(j++\), 
IntervalLength = \((ParamTensor[\([j, 3, 2]\)] - ParamTensor[\([j, 3, 1]\)])\)/
            NN; 
IntervalPoints\  = \ {ParamTensor[\([j, 3, 
                1]\)]}; 
For[k = 1, \ 
          k ≤ 
            NN, \(k++\), \ \(IntervalPoints = 
              Append[IntervalPoints, \ 
                ParamTensor[\([j, 3, 1]\)] + 
                  k*IntervalLength];\)
\ \ ]; 

 (*Perform\ Sampling\ on\ the\ parameter\ intervals*) 

 (*Generate\ List\ of\ random\ numbers\ to\ choose\ intervals*) 
RandomIntervals = 
          Random[]; 
For\ [k = 1, \ 
            k < NN, \ \(k++\), 
\(RandomIntervals = {RandomIntervals, \ 
                  Random[Real, \ 
                    1 - k/NN]};\)
]
\(RandomIntervals = 
            Flatten[RandomIntervals]\); 

 (*Come\ up\ with\ an\ ordering\ of\ parameter\ intervals\ based\ on\ the\ list\ of\ random\ numbers\ generated . \ \ The\ ordering\ is\ generated\ without\ replacement*) 
SortIntervals = {}; 

        For[k = 1, \ 
          k <= NN, \ \(k++\), 
\(SortIntervals = {SortIntervals, \ 
                k};\)
]; 
SortIntervals = 
          Flatten[SortIntervals]; 


        For[k = 0, \ 
          k < NN, \ \(k++\), 
BinPicked = 
            Ceiling[\((NN)\)*
                RandomIntervals[\([k + 1]\)]\ ]; 
Swap = 
            SortIntervals[\([BinPicked + 
                  k]\)]; 
SortIntervals[\([BinPicked + 
                  k]\)] = 
            SortIntervals[\([1 + 
                  k]\)]; 
SortIntervals[\([1 + k]\)] = 
            Swap;
]; 

 (*Randomly\ select\ values\ from\ the\ intervals\ chosen*) 
ParamValue = {}; 

        For[k = 1, \ 
          k ≤ \ 
            NN, \ \(k++\), 
\(ParamValue = {ParamValue, \ 
                Random[Real, \ {IntervalPoints[\([SortIntervals[\([k]\)]]\)], \ IntervalPoints[\([SortIntervals[\([k]\)] + 
                          1]\)]}]};\)
]; 

ParamValue = 
          Flatten[ParamValue]; 
ParamTensor[\([j, 4]\)] = 
          ParamValue;
];\)\)
ParamTensor;

q = 
        1;
\(Do[
 (*Pick\ out\ a\ set\ of\ parameters\ from\ the\ ParamTensor*) 
a = 
        Rationalize[ParamTensor[\([1, 4, q]\)]]; 

      b\  = \ Rationalize[ParamTensor[\([2, 4, q]\)]]; 

      c\  = Rationalize[
          ParamTensor[\([3, 4, q]\)]]; 

ode1 = \ \(P'\)[
            t] == P[t] \((1.  - P[t] - Q[t])\); 

      ode2 = \ \(Q'\)[t]\  == 
          Q[t] \((a - b\ Q[t] - 
                c\ P[t])\); 

sol = 
        NDSolve[{ode1, \ ode2, \ P[0. ] == 1.25, \ Q[0. ] == 1.25}, \ {P[t], 
            Q[t]}, \ {t, 0, solvetime}]; 
SolutionTensor[
          q]\  = sol;

, \ {q, \ 
        NN}]\ \)

SolNum  = 5;
\(\(p1 = 
        Plot[Evaluate[{P[t]} /. SolutionTensor[SolNum]], \ {t, 0, 
            solvetime}, 
		\ PlotStyle -> {RGBColor[1, 0, 
                0]}, 
		PlotLabel -> {"P1 Population size", a, 
              Rationalize[ParamTensor[\([2, 4, SolNum]\)]], 
              Rationalize[
                ParamTensor[\([3, 4, 
                    SolNum]\)]]}];\)
\)
\(\(p2 = 
        Plot[Evaluate[{Q[t]} /. SolutionTensor[SolNum]], \ {t, 0, 
            solvetime}, 
		\ PlotStyle -> {RGBColor[1, 0, 0]}, 
          PlotRange -> {0, 
              1}, 
		PlotLabel -> {"P2 Population size", a, 
              Rationalize[ParamTensor[\([2, 4, SolNum]\)]], 
              Rationalize[ParamTensor[\([3, 4, SolNum]\)]]}];\)\)