(*Parametry Skellamova rozdělení*)
lambda1True=7;
lambda2True=7;

(*Funkce pro odhad pomocí metody momentů*)
methodOfMoments[data_]:=Module[{mean,variance},mean=Mean[data];
variance=Variance[data];
{(variance+mean)/2//N,(variance-mean)/2}//N]

(*Funkce pro simulace a odhady pomocí MM*)
simulateAndEstimateMM[sampleSize_]:=Module[{data,methodMomentsEstimate},
SeedRandom[7];
data=RandomVariate[SkellamDistribution[lambda1True,lambda2True],{500,sampleSize}];
methodMomentsEstimate=Table[methodOfMoments[data[[i]]],{i,1,500}];
methodMomentsEstimate]

(*Log-likelihood funkce*)
logL[lambda1_,lambda2_,data_]:=Module[{n,sumX},n=Length[data];
sumX=Total[data];
-n (lambda1+lambda2)+(Log[lambda1]-Log[lambda2]) (sumX/2)+Sum[Log[BesselI[Abs[data[[i]]],2 Sqrt[lambda1 lambda2]]],{i,1,n}]]

(*Globální maximalizace Log-likelihood funkce*)
optimizeForSample[data_]:=Module[{result,lambda1,lambda2},result=NMaximize[{logL[lambda1,lambda2,data],lambda1>0,lambda2>0},{lambda1,lambda2},Method->"NelderMead"];
{lambda1,lambda2}/. result[[2]]]

(*Funkce pro simulace a odhady pomocí MLE*)
simulateAndEstimateMLE[sampleSize_]:=Module[{data,MLEEstimateResult},
SeedRandom[7];
data=RandomVariate[SkellamDistribution[lambda1True,lambda2True],{500,sampleSize}];
(*Odhady metodou maximální věrohodnosti*)MLEEstimateResult=Table[optimizeForSample[data[[i]]],{i,1,500}];
MLEEstimateResult]

(*Rozsah 50*)
dvojiceMM10= simulateAndEstimateMM[50];
dvojiceMLE10 = simulateAndEstimateMLE[50];
lambdaMLE110 = Table[dvojiceMLE10[[i,1]],{i,1,500}];
lambdaMLE210 = Table[dvojiceMLE10[[i,2]],{i,1,500}];
lambdaMM110 = Table[dvojiceMM10[[i,1]],{i,1,500}];
lambdaMM210 = Table[dvojiceMM10[[i,2]],{i,1,500}];

biasMM10=Mean[dvojiceMM10-ConstantArray[{lambda1True,lambda2True},500]];
biasMLE10=Mean[dvojiceMLE10-ConstantArray[{lambda1True,lambda2True},500]];

varMM10=Variance[dvojiceMM10];
varMLE10=Variance[dvojiceMLE10];

mseMM10=varMM10+biasMM10^2;
mseMM10//Mean;
mseMLE10=varMLE10+biasMLE10^2;
mseMLE10//Mean;

(*Rozsah 100*)
dvojiceMM50= simulateAndEstimateMM[100];
dvojiceMLE50 = simulateAndEstimateMLE[100];
lambdaMLE150 = Table[dvojiceMLE50[[i,1]],{i,1,500}];
lambdaMLE250 = Table[dvojiceMLE50[[i,2]],{i,1,500}];
lambdaMM150 = Table[dvojiceMM50[[i,1]],{i,1,500}];
lambdaMM250 = Table[dvojiceMM50[[i,2]],{i,1,500}];

biasMM50=Mean[dvojiceMM50-ConstantArray[{lambda1True,lambda2True},500]];
biasMLE50=Mean[dvojiceMLE50-ConstantArray[{lambda1True,lambda2True},500]];

varMM50=Variance[dvojiceMM50];
varMLE50=Variance[dvojiceMLE50];

mseMM50=varMM50+biasMM50^2;
mseMM50//Mean;
mseMLE50=varMLE50+biasMLE50^2;
mseMM50=varMM50+biasMM50^2;
mseMLE50//Mean;

(*Rozsah 250*)
dvojiceMM100= simulateAndEstimateMM[250];
dvojiceMLE100 = simulateAndEstimateMLE[250];
lambdaMLE1100 = Table[dvojiceMLE100[[i,1]],{i,1,500}];
lambdaMLE2100 = Table[dvojiceMLE100[[i,2]],{i,1,500}];
lambdaMM1100 = Table[dvojiceMM100[[i,1]],{i,1,500}];
lambdaMM2100 = Table[dvojiceMM100[[i,2]],{i,1,500}];

biasMM100=Mean[dvojiceMM100-ConstantArray[{lambda1True,lambda2True},500]];
biasMLE100=Mean[dvojiceMLE100-ConstantArray[{lambda1True,lambda2True},500]];

varMM100=Variance[dvojiceMM100];
varMLE100=Variance[dvojiceMLE100];

mseMM100=varMM100+biasMM100^2;
mseMLE100=varMLE100+biasMLE100^2;
mseMLE100//Mean;

(*Rozsah 500*)
dvojiceMM1000= simulateAndEstimateMM[500];
dvojiceMLE1000 = simulateAndEstimateMLE[500];
lambdaMLE11000 = Table[dvojiceMLE1000[[i,1]],{i,1,500}];
lambdaMLE21000 = Table[dvojiceMLE1000[[i,2]],{i,1,500}];
lambdaMM11000 = Table[dvojiceMM1000[[i,1]],{i,1,500}];
lambdaMM21000 = Table[dvojiceMM1000[[i,2]],{i,1,500}];

biasMM1000=Mean[dvojiceMM1000-ConstantArray[{lambda1True,lambda2True},500]];
biasMLE1000=Mean[dvojiceMLE1000-ConstantArray[{lambda1True,lambda2True},500]];

varMM1000=Variance[dvojiceMM1000];
varMLE1000=Variance[dvojiceMLE1000];

mseMM1000=varMM1000+biasMM1000^2;
mseMM1000//Mean;
mseMLE1000=varMLE1000+biasMLE1000^2;
mseMLE1000//Mean;

(*Rozsah 1000*)
dvojiceMM5000= simulateAndEstimateMM[1000];
dvojiceMLE5000 = simulateAndEstimateMLE[1000];
lambdaMLE15000 = Table[dvojiceMLE5000[[i,1]],{i,1,500}];
lambdaMLE25000 = Table[dvojiceMLE5000[[i,2]],{i,1,500}];
lambdaMM15000 = Table[dvojiceMM5000[[i,1]],{i,1,500}];
lambdaMM25000 = Table[dvojiceMM5000[[i,2]],{i,1,500}];

biasMM5000=Mean[dvojiceMM5000-ConstantArray[{lambda1True,lambda2True},500]];
biasMLE5000=Mean[dvojiceMLE5000-ConstantArray[{lambda1True,lambda2True},500]];

varMM5000=Variance[dvojiceMM5000];
varMLE5000=Variance[dvojiceMLE5000];

mseMM5000=varMM5000+biasMM5000^2;
mseMM5000//Mean;
mseMLE5000=varMLE5000+biasMLE5000^2;
mseMLE5000//Mean;

(*Střední hodnota odhadů MM*)
estimatesMM1 = {lambdaMM110 //Mean, lambdaMM150//Mean,lambdaMM1100//Mean,lambdaMM11000//Mean , lambdaMM15000//Mean  };
estimatesMM2 = {lambdaMM210 //Mean, lambdaMM250//Mean,lambdaMM2100//Mean,lambdaMM21000//Mean , lambdaMM25000//Mean  };

(*Střední hodnota odhadů MLE*)
estimatesMLE1 = {lambdaMLE110 //Mean, lambdaMLE150//Mean,lambdaMLE1100//Mean,lambdaMLE11000//Mean , lambdaMLE15000//Mean };
estimatesMLE2 = {lambdaMLE210 //Mean, lambdaMLE250//Mean,lambdaMLE2100//Mean,lambdaMLE21000//Mean , lambdaMLE25000//Mean  };

(*Konsistence*)
xValues={50,100,250,500,1000};
ListLinePlot[{Transpose[{xValues,estimatesMM1}],Transpose[{xValues,estimatesMLE1}],Transpose[{xValues,ConstantArray[lambda1True,Length[xValues]]}]},PlotRange->All,ScalingFunctions->{"Log",None},Frame->True,FrameLabel->{"Sample Size","Estimate"},FrameTicks->{Automatic,Automatic}]
ListLinePlot[{Transpose[{xValues,estimatesMM2}],Transpose[{xValues,estimatesMLE2}],Transpose[{xValues,ConstantArray[lambda2True,Length[xValues]]}]},PlotRange->All,ScalingFunctions->{"Log",None},Frame->True,FrameLabel->{"Sample Size","Estimate"},FrameTicks->{Automatic,Automatic}]

(*Boxploty*)
GraphicsRow[{Show[BoxWhiskerChart[{lambdaMM110,lambdaMLE110}],Graphics[{Red,Dashed,Line[{{0.5,lambda1True},{2.5,lambda1True}}]}]],
Show[BoxWhiskerChart[{lambdaMM150,lambdaMLE150}],Graphics[{Red,Dashed,Line[{{0.5,lambda1True},{2.5,lambda1True}}]}]]}]
GraphicsRow[{Show[BoxWhiskerChart[{lambdaMM1100,lambdaMLE1100}],Graphics[{Red,Dashed,Line[{{0.5,lambda1True},{2.5,lambda1True}}]}]],
Show[BoxWhiskerChart[{lambdaMM11000,lambdaMLE11000}],Graphics[{Red,Dashed,Line[{{0.5,lambda1True},{2.5,lambda1True}}]}]]}]
GraphicsRow[{Show[BoxWhiskerChart[{lambdaMM15000,lambdaMLE15000}],Graphics[{Red,Dashed,Line[{{0.5,lambda1True},{2.5,lambda1True}}]}]]}]

GraphicsRow[{Show[BoxWhiskerChart[{lambdaMM210,lambdaMLE210}],Graphics[{Red,Dashed,Line[{{0.5,lambda2True},{2.5,lambda2True}}]}]],
Show[BoxWhiskerChart[{lambdaMM250,lambdaMLE250}],Graphics[{Red,Dashed,Line[{{0.5,lambda2True},{2.5,lambda2True}}]}]]}]
GraphicsRow[{Show[BoxWhiskerChart[{lambdaMM2100,lambdaMLE2100}],Graphics[{Red,Dashed,Line[{{0.5,lambda2True},{2.5,lambda2True}}]}]],
Show[BoxWhiskerChart[{lambdaMM21000,lambdaMLE21000}],Graphics[{Red,Dashed,Line[{{0.5,lambda2True},{2.5,lambda2True}}]}]]}]
GraphicsRow[{Show[BoxWhiskerChart[{lambdaMM25000,lambdaMLE25000}],Graphics[{Red,Dashed,Line[{{0.5,lambda2True},{2.5,lambda2True}}]}]]}]

(*Příklad konkrétních hodnot pro lambda1MM,lambda2MM,lambda1MLE,lambda2MLE*)
rozsahy = {50,100,250,500,1000};
lambda1MMValues=estimatesMM1; (*Např.hodnoty pro lambda1MM*)
lambda2MMValues=estimatesMM2; (*Např.hodnoty pro lambda2MM*)
lambda1MLEValues=estimatesMLE1; (*Např.hodnoty pro lambda1MLE*)
lambda2MLEValues=estimatesMLE2; (*Např.hodnoty pro lambda2MLE*)

(*Vytvoření tabulky*)
table=Transpose[{rozsahy, lambda1MMValues,lambda2MMValues,lambda1MLEValues,lambda2MLEValues}];

(*Tabulka výběrových průměrů odhadů pro daný rozsah*)
Grid[Prepend[table,{"rozsah","lambda1MM","lambda2MM","lambda1MLE","lambda2MLE"}],Spacings->{2,2}]

(*Příklad konkrétních hodnot pro lambda1MM,lambda2MM,lambda1MLE,lambda2MLE*)
lambda1MMBias={biasMM10[[1]],biasMM50[[1]],biasMM100[[1]],biasMM1000[[1]],biasMM5000[[1]]}; (*Např.hodnoty pro lambda1MM*)
lambda2MMBias={biasMM10[[2]],biasMM50[[2]],biasMM100[[2]],biasMM1000[[2]],biasMM5000[[2]]}; (*Např.hodnoty pro lambda2MM*)
lambda1MLEBias={biasMLE10[[1]],biasMLE50[[1]],biasMLE100[[1]],biasMLE1000[[1]],biasMLE5000[[1]]}; (*Např.hodnoty pro lambda1MLE*)
lambda2MLEBias={biasMLE10[[2]],biasMLE50[[2]],biasMLE100[[2]],biasMLE1000[[2]],biasMLE5000[[2]]}; (*Např.hodnoty pro lambda2MLE*)

(*Vytvoření tabulky*)
table=Transpose[{rozsahy, lambda1MMBias,lambda2MMBias,lambda1MLEBias,lambda2MLEBias}];

(*Tabulka BIAS odhadů pro daný rozsah*)
Grid[Prepend[table,{"rozsah","Bias lambda1MM","Bias lambda2MM","Bias lambda1MLE","Bias lambda2MLE"}],Spacings->{2,2}]

(*Příklad konkrétních hodnot pro lambda1MM,lambda2MM,lambda1MLE,lambda2MLE*)
lambda1MMVar={varMM10[[1]],varMM50[[1]],varMM100[[1]],varMM1000[[1]],varMM5000[[1]]}; (*Např.hodnoty pro lambda1MM*)
lambda2MMVar={varMM10[[2]],varMM50[[2]],varMM100[[2]],varMM1000[[2]],varMM5000[[2]]}; (*Např.hodnoty pro lambda2MM*)
lambda1MLEVar={varMLE10[[1]],varMLE50[[1]],varMLE100[[1]],varMLE1000[[1]],varMLE5000[[1]]}; (*Např.hodnoty pro lambda1MLE*)
lambda2MLEVar={varMLE10[[2]],varMLE50[[2]],varMLE100[[2]],varMLE1000[[2]],varMLE5000[[2]]}; (*Např.hodnoty pro lambda2MLE*)

(*Vytvoření tabulky*)
table=Transpose[{rozsahy, lambda1MMVar,lambda2MMVar,lambda1MLEVar,lambda2MLEVar}];

(*Tabulka výběrových rozptylů odhadů pro daný rozsah*)
Grid[Prepend[table,{"rozsah","Var lambda1MM","Var lambda2MM","Var lambda1MLE","Var lambda2MLE"}],Spacings->{2,2}]

(*Příklad konkrétních hodnot pro lambda1MM,lambda2MM,lambda1MLE,lambda2MLE*)

lambda1MMmse={mseMM10[[1]],mseMM50[[1]],mseMM100[[1]],mseMM1000[[1]],mseMM5000[[1]]}; (*Např.hodnoty pro lambda1MM*)
lambda2MMmse={mseMM10[[2]],mseMM50[[2]],mseMM100[[2]],mseMM1000[[2]],mseMM5000[[2]]}; (*Např.hodnoty pro lambda2MM*)
lambda1MLEmse={mseMLE10[[1]],mseMLE50[[1]],mseMLE100[[1]],mseMLE1000[[1]],mseMLE5000[[1]]}; (*Např.hodnoty pro lambda1MLE*)
lambda2MLEmse={mseMLE10[[2]],mseMLE50[[2]],mseMLE100[[2]],mseMLE1000[[2]],mseMLE5000[[2]]}; (*Např.hodnoty pro lambda2MLE*)

(*Vytvoření tabulky*)
table=Transpose[{rozsahy, lambda1MMmse,lambda2MMmse,lambda1MLEmse,lambda2MLEmse}];

(*Tabulka MSE odhadů pro daný rozsah*)
Grid[Prepend[table,{"rozsah","MSE lambda1MM","MSE lambda2MM","MSE lambda1MLE","MSE lambda2MLE"}],Spacings->{2,2}]