Quiet@Remove["Global`*"];
SetDirectory[NotebookDirectory[]];
(*jelikoz "Full" je hezci, ale pripravuji texty take na laptopu, ktery je stary a snese maximalne Mathematica 7.0 a ta moznost "Full" nema, udelam si globalni promennou gridLines abych pri prechodu mezi verzemi menil jen jednu promennou na jednom miste*)
gridLines=Full;

(*Data obvykle v technicke praxi alespon napoprve nedostaneme v nejake privetive podobe*)
(*Tady ukazeme jen, ze mame kumulativni soucty jednne veliciny a denni prirustky druhe veliciny.*)
(*Vetsinou je to v praxi mnohem horsi, formaty dat, nestejny krok dat, podivne znaky atd...*)

(*z webu "https://onemocneni-aktualne.mzcr.cz/covid-19?fbclid=IwAR2l7Yo-jDxnsKoCf-3O0X5BC9DZMNi8hU3qSUbWWjrJHlgFGIJycOy3adg" mame od 1.3. odectene denni prirustky pozitivne testovanych na covid*)
(*promenna positvniVeDni se mela jmenovat positivniVeDni, ale kdyz uz bylo vsechno napsane, je zbytecne to menit*)
(*Casovy krok mezi dvema hodnotami v listu je jeden den, muzeme to chapat jako posloupnost nebo funkci defiovanou tabulkou*)
positvniVeDni={0,3,0,3,11,7,6,6,25,31,22,25,48,109,85,67,110,205,124,158,114,126,185,291,259,373,262,160};

(*Promenna testuchKumulatiune (jmeno promenne inspirovae Zemeplochou) maji opet casovy odstup jeden den a jde o celkovy pocet provedenych testu, tedy kumulativni sumu ci integral jako funkce horni meze, distribucni funkci v pravdepodobnosti,jak je libo *)
testuchKumulatiune={201,262,340,407,483,594,787,928,1193,1358,1816,2356,3094,4065,5068,6302,7664,9402,11619,13704,15584,17377,19624,22601,26698,31127,36374,40700,43498};
(*mame-li vice souboru dat, podivame se na jejich delku, zde namapujeme fukci Length a list se dvema prvky, zminenymi soubory dat*)
(*otazka pro zvidave hlavy: proc potrebujeme mit akumulovanych o jeden prvek vice abychom popsali stejnou situaci? Nesouvisi to nahodou s pocatecni podminkou diferencialni rovnice prvniho radu? *)
{nap,nkt}=Length/@{positvniVeDni,testuchKumulatiune}

(*mame-li data ruznych formatu, neni spatne si je ztransformovat tak, abychom meli k dispozici pro kazdy JEV stejny format*)
(*data v promenne positvniVeDni jsou funkcni hodnoty a v promenne testuchKumulatiune kumulativni hodnoty, vyrobime si tedy kumulativni hodnoty pozitivne testovanych*)
(*chcete-li videt, jak funguje Accumulate, zkuste Help nebo prikaz Accumulate[(Subscript[a, #])&/@Range[5]] a jeho modifikace...*)

akuPositivni=Accumulate[positvniVeDni]

(*bylo to jednoduche... pro kumulativni data jsme pouzili Accumulate*)
(*ted naopak ze souctu {(Subscript[a, 1]),(Subscript[a, 1])+(Subscript[a, 2]),(Subscript[a, 1])+(Subscript[a, 2])+(Subscript[a, 3]),(Subscript[a, 1])+(Subscript[a, 2])+(Subscript[a, 3])+(Subscript[a, 4]),(Subscript[a, 1])+(Subscript[a, 2])+(Subscript[a, 3])+(Subscript[a, 4])+(Subscript[a, 5])} potrebujeme udelat {(Subscript[a, 1]),(Subscript[a, 2]),(Subscript[a, 3]),(Subscript[a, 4]),(Subscript[a, 5])} *)
(*ruzne rychlych ci elegantnich moznosti je ve Wolframu mnoho, nam se libi postupy zacinajici tim, ze si udelame takovehle dvojice*)
(*pp=Accumulate[(Subscript[a, #])&/@Range[5]];
dvojice=Partition[pp,2,1]
{{(Subscript[a, 1]),(Subscript[a, 1])+(Subscript[a, 2])},{(Subscript[a, 1])+(Subscript[a, 2]),(Subscript[a, 1])+(Subscript[a, 2])+(Subscript[a, 3])},{(Subscript[a, 1])+(Subscript[a, 2])+(Subscript[a, 3]),(Subscript[a, 1])+(Subscript[a, 2])+(Subscript[a, 3])+(Subscript[a, 4])},{(Subscript[a, 1])+(Subscript[a, 2])+(Subscript[a, 3])+(Subscript[a, 4]),(Subscript[a, 1])+(Subscript[a, 2])+(Subscript[a, 3])+(Subscript[a, 4])+(Subscript[a, 5])}}
*)
(*vidime, ze odecteme-li od druheho prvku kazde dvojice prvek prvni, ziskame {(Subscript[a, 1]),(Subscript[a, 2]),(Subscript[a, 3]),(Subscript[a, 4]),(Subscript[a, 5])}*)
(*tyto rozdily muzeme vyrobit ruzne, ukazeme nektere *)
(*tenhle prvni vyuziva skalarniho soucinu, funkce {-1,1}.#& vraci rozdil druheho a prvniho prvku parametru (ktery musi byt List)*)
(*zkuste si:
{-1,1}.#&[{x,y}]
vrati to 
-x+y
*)
veDnechTesty={-1,1}.#&/@Partition[testuchKumulatiune,2,1]

(*nebo jinak, kazdy mame jine preference, jinak myslime... tak treba pomoci patternu a pravidel zameny*)
veDnechTesty=Partition[testuchKumulatiune,2,1]/.{x_,y_}:>y-x

(*jasne, v ramci exhibicionismu a fukcni libosti muzeme treba takhle :-)*)
veDnechTesty=First/@ListConvolve[{{1,-1}},Partition[testuchKumulatiune,2,1]]

(*jasne, milovnici proceduralnejsich zpusobu programovani to mohou napsat treba takhle*)
veDnechTesty=Table[testuchKumulatiune[[i+1]]-testuchKumulatiune[[i]],{i,Length[testuchKumulatiune]-1}]

(*a pokud vas neco pekneho napadne, jak to napsat jinak, poslete nema to a my vam radi dame body! :-0*)

(*better safe than sorry... kdyz jsme vyrobili listy stejnych vyznamu (kumulativni/denniPrirustek), zkontrolujeme, jestli maji stejnou delku, stejny pocet prvku, stejnou dimenzi, chapeme-li je jako vektor *)
(*jasne, ted to sedi... ale napoprve a napodruhe, kdyz jsme tenhle program psali, tak to nesedelo... chybicky se vloudi a to, ze vysledek vypada logicky, neznamena, ze to tak bylo na zacatku psani kodu*)
Length/@{positvniVeDni,veDnechTesty}

(*ted troska programovani...  zatim jsme meli "den" jen mysleny odstup mezi dvema sousednimi prvky listuu *)
(*pro dalsi ale bude vyhodne mit misto dat typu {a,b,c,d...} data {{1,a},{2,b},{3,c},{4,d},..}*)
(*priprava na to... sice jsme se vizualne presvedcili, ze delky listu jsou stejne, presto si radeji vezmu minimum jejich delek -zde totozne s jejich delkami/delkou- a mam jistotu, ze neprelezu mimo data...*)
n=Min[Length/@{positvniVeDni,veDnechTesty}];

(*vytvorim si tedy promennou "dny" *)
dny=Range[n]

(*a ted jeden postup jak udelat z listu {1,2,3,4,..} a listu {a,b,c,d...} list {{1,a},{2,b},{3,c},{4,d},..}*)
(*to "T" nahore vpravo znamena Transpose a napise se jako "esc tr esc"*)
(*klasicka transpozice matice, z matice
MatrixForm[Take[#,5]&/@{dny,positvniVeDni}]
tedy
(({{1, 2, 3, 4, 5}, {0, 3, 0, 3, 11}}))
vrati
MatrixForm[(Take[#,5]&/@{dny,positvniVeDni})\[Transpose]]
tedy
(({{1, 0}, {2, 3}, {3, 0}, {4, 3}, {5, 11}}))
*)
(*Take[#,5]& jen vybira prvnih pet prvku z listu, ktery je jejim parametrem...a bychom neukazovali zbytecne cela data, tahle ukazka by mela k pochopeni stacit*)
(*{dny,positvniVeDni})\[Transpose] je ekvivalentni Transpose[{dny,positvniVeDni}]*)

dataAku={dny,positvniVeDni}\[Transpose];

(*totez pro akumulovana data pozitivnich pripadu*)
(*a vzdycky je fajn se na data podivat!*)
dataAkuPos={dny,akuPositivni}\[Transpose];
pl1=ListPlot[dataAkuPos,AxesLabel->{"den","nPositivnich"},PlotStyle->{Red,PointSize[0.01]},GridLines->gridLines,PlotRange->All]

(*no a ted zacneme samanit, hrat si... podivejte, neni vubec zadny duvod proc by mel byt celkovy pocet pozitivne testovanych vyjadritelny vzorcem, navic jednoduchym vzorcem*)
(*pokud jsme deterministe a pozitiviste (jako vetsina techniku), pripustime, ze by mohla byt funkci zavislost mezi vsemi pozitivnimi (nezavisle na tom, jestli byli testovani, nebo ne) a casem*)
(*takze tenhle model je velika PITOMOST!*)
(*ale prosim, uvdomte si, ze pitomost jsou i modely, o kterych se doctete na netu, opravdu*)
(*takze si klidne zkuste jine modely! Budou se treba pekne motat kolem bodu... ale o skutecnosti (pokud jsme pozitivisti a predpokladame existenci, ci lepe byti, skutecnosti) nevypovidaji nic*)
(*mozna proto jsme zvolili docela hloupiy vzorecek*)

model[den_]:=a*(Exp[((b*den+f*(den^2))/(1+d*den+e*(Sqrt[den])))]-1);
(*opet: bereme eukleidovskou normu, ale proc? Jen ze zvyku. Neni pro ni zadny duvod: nejde o Fit (linearni regresi, linearni problem), kde mela misto, nebot po derivovani se problem hledani minima redukoval na reseni soustavy linearnich rovnic. Ted je to jedno... pomuze, kdyz ma norma vsude derivaci podle sveho parametru... ale to uz je spis magisterska latka z predmetu IAP :-)*)

norm:=((#)^2)&;
(*naprosto analogicky minulemu textiku... *)
err[{den_,nPositivnich_}]:=norm[model[den]-nPositivnich];
errCelk=Total[err/@dataAkuPos];
(*najdeme si promenne, ktere hledame, to, co vidite modre v definici model *)
(*no jo... jenze kdyz to takhle delame, ze hledame promenne typu Symbol (pattern _Symbol), tak najdeme i E, protoze tohle iracionalni cislo je ve Wolframu typu Symbol... stejne jako treba π. Tak ho zrusime! :-) *)
(*zrusime ho pomoci DeleteCases[kdeZrusit,pattern] *)
(*a E je zcela konkretni pattern, konkretnost neubira patternu pou patternovitosti :-)*)
(*jelikoz Cases najde vsechny vyskyty patternu v cele struktire, poskrtame duplicity pomoci Union*)
promsMin=DeleteCases[Union@Cases[errCelk,_Symbol,{0,∞}],E]
(*nebo lepe, kdybychom chteli odstranit vse, co neni cislo... podivejte, ve Wolframu jsou Pi, E a podobne reprezentovane jako Symbol, ale mohou mit ciselnou hodnotu*)
(*z boolske funkce bool udelame pattern ?bool, && je klasicke And jako v jinych jazycich, vykricnik znamena negaci, !True==False *)
promsMin=Union@Cases[errCelk,_?((Head[#]==Symbol)&&!NumericQ[#]&),{0,∞}]
(*tady mame totez zapsane "matemematickologicteji"*)
promsMin=Union@Cases[errCelk,_?(And[Head[#]==Symbol,Not[NumericQ[#]]]&),{0,∞}]

(*tak si najdeme optimalni hodnoty parametru modelu*)
(*treti parametr NMinimize, tedy Method->"DifferentialEvolution", rika, ze se ma hledat minimum pomoci genetickych algoritmu... snad bude cas a v lepsich casech vam povim, co to znamena, zatim zkuste google*)
(*bereme NMinimize[errCelk,promsMin,Method->"DifferentialEvolution"][[2]]  protoze NMIinimize vraci List {hodnotaMinima,{a->aMin,...}}*)
(*a my potrebujeme jen ten druhy prvek toho listu, tedy {a->aMin,...}*)
min=NMinimize[errCelk,promsMin,Method->"DifferentialEvolution"][[2]];
(*ted uz je to jednoduche: vys je vyraz, kde jsme do vyrazu model[den] dosadili za "modre" promenne *)
vys=model[den]/.min;
pl2=Plot[vys,{den,1,n}];
pl5=Show[pl1,pl2,PlotRange->All]
(*jeste si pro "prosty lid" vypocteme limitu, jak by to skoncilo, kdyby platil uvazovany model*)
Print["nakonec bude  ",Limit[vys,den->∞],"  nakazenych"]
(*a vidite, ze pitomy model dal pitomou limitu... CR bude tezko nekomu dluzit dvanact setin nakazeneho :-)*)

(*a ted prejdeme od vlastniho obsahu VYA k vlastnimu obsahu VYA... zacneme premyslet:-)*)
(*na tyhle grafy se podivame a... vykasleme se na ne*)
(*lepe se premysli v prirustcich, nez v akumulovanych datech*)
(*LINEARNI difoska ma reseni nelinearni vuci parametru*)

Show[ListPlot[testuchKumulatiune,Joined->True,GridLines->gridLines],ListPlot[testuchKumulatiune],pl5,PlotRange->All,AxesLabel->{"den","n"}]

(*spocteme si koeficient korelace... a vysvetlime si, co na  rika*)

corr=N@Correlation[positvniVeDni,veDnechTesty];
Print["koeficient korelace=  ",corr," -" ]

(*koeficient korelace je cislo mezi munus 1 a 1 a dava nam miru linearni zavislosti mezi dvema vektory (zde Listy) dat stejne delaky*)
(*zkusme treba*)
nDat=10;
SeedRandom[1];
data=Table[Random[],{nDat}];
funkcniZavislost1[x_]:=0.1+5*x;
data1=funkcniZavislost1/@data;
(*vidime cistou jednicku: data jsou svazana linearni funkci s kladnou derivaci (to je ta petka v definici funkcniZavislost1[x_]:=0.1+5*x)*)
Correlation[data,data1]

(*podobne*)
funkcniZavislost2[x_]:=0.1-5*x;
data2=funkcniZavislost2/@data;
(*vidime cistou minus jednicku: data jsou svazana linearni funkci se zapornou derivaci (to je ta minus petka v definici funkcniZavislost1[x_]:=0.1+5*x)*)
Correlation[data,data2]

(*pridejme nahodna cisla *)
funkcniZavislost3[x_]:=0.1-5*x+3*RandomReal[{-1,1}];
data3=funkcniZavislost3/@data;
(*vidime ze data prilis linearne svazana nejsou, zaporna hodnota ale znamena, ze prirustek prvnich dat spise provazi pokles druhych*)
Correlation[data,data3]

(*pridejme  nelinearitu*)
funkcniZavislost4[x_]:=Exp[10x];
data4=funkcniZavislost4/@data;
(*vidime jakesi kladne cislo, narust prvni veliciny spise vede k narustu druhe, moc to ale linearni zavislost neni *)
Correlation[data,data4]

(*vytvorme si list pomeru pocetPositivneTestovanych/pocetVsechTestu*)
(*podil listu vraci podily "po prvcich", tedy (({a,b,c})/({d,e,f})) vrati {(a/d),(b/e),(c/f)}*)
pomer=(positvniVeDni/veDnechTesty)//N

(*pro srovnani za cele obdobi vezmeme celkovy pocet pozitivnich vysledku/celkovy pocet testu*)
(*v kumulativnich datech jsou to posledni prvky prislusnych listu*)
(*to dalsi jsou jen grafy, to uz byste meli umet*)
mean=((Last[akuPositivni])/(Last[testuchKumulatiune]))//N;
pl3=ListPlot[pomer,PlotRange->All,AxesLabel->{"den","(nPositivnich/nTestu)"},PlotStyle->{Red, PointSize[0.01]},Joined->True,PlotLabel->"mean= "<>ToString[mean]<>"   corr= "<>ToString[0.001*Round[1000corr]],GridLines->gridLines];
pl4=ListPlot[mean&/@positvniVeDni,Joined->True,PlotStyle->{Black,Thickness[0.004]}];
plvys=Show[pl3,pl4]

(*tak a ted se pustime do ... do klasicke prace inzenyra, ktery nema dost dat, eni odbornik na zadnou cast problemu (epidemiologie, matematicka statistika a podobne), ale ma prehled aspon trochu a chce se v problemu orientovat*)

(*na netu najdete moudra o "R faktoru", coz neni nic nez kvocient gemotricke posloupnost pozitivne testovanych*)
(*opravdu neni v normalnich podminkach duvod, proc by melo jit o geometrickou posloupnost (natoz radu, jak nejaci "inteligenti" pisou)*)
(*podivejme se na to matematicteji, logicteji, technicteji...*)

(*1/ pocet pozitivne testovanych je vlastne nezajimava velicina... pokud se ve state Dolni Ampera (omlouvam se Horni Volte, podobnost ciste nahodna) nebude testovat vubec, bude to vypadat, ze to nikdo nema
2/dulezity je totiz počet skutečně pozitivnich a nikoli pozitivne testovanych:
Neumira se na pozitivni test, ale na nemoc samotnou. Predstavte si,Ze posuzujete mnozství hub v lesich podle toho,kolik houbari nasbiraji hub v danem dni a neberete v uvahu pocet houbarohodin v lesich. Dojdete k zvlastnimu objevu,ze houby se v lesich vyskytuji zejmena o weekendech. Proto vyhodnocujeme korelaci pozitivních vysledku testu s poctem provedenych testu... kdyby byla mala,informaci o promorenosti populace neziskam. Jenze ona vychazi velmi blizka jednicce a to uz nám informaci dava, resp. dava ji dohromady s informaci, ze pozitivnich je pod 7%,tedy se testuje vice z hlediska covid19 zdravych,než pozitivnich. Dokud to tak bude,lze usuzovat,že promorenost populace nestoupa, tak jako nestoupa pocet hub v lesich o weekendech,ackoli je jich o weekendech vice nalezeno.*)
(*samozrejme mame malo dat, nevime, kolik bylo opakovanych testu u jedne osoby, nakoli se testuji hypochondri bez priznaku a podobne*)

(*a na konec?*)
(*studenty obvykle zajimaji body a znamka...*)
(*muzete doplnovat data, hledat rozumnejsi modely a poloci limity tipovat vyvoj situace*)
(*za dobra reseni radi prodelime body*)

(*a nebojte se moc, podle koeficientu korelace a pomeru to vypada docela nadejne*)

(*zaverem jen, ze ty pomery uz oficialni stranka zverejnuje, tu korelaci ne... ale puvodne tam nebyly ani ty pomery a tak jsem nam to naprogramoval*)