ClearAll["Global`*"]

Метод Монте-Карло

Нормальное распределение

Clear[μ,x]
σ=1
Table[PDF[NormalDistribution[μ,σ],x],{μ,-1,2}]
Plot[%,{x,-5,5},Filling->Axis]

Clear[λ,x]
μ=0
Table[PDF[NormalDistribution[μ,σ],x],{σ,1/2,2,1/2}]
Plot[%,{x,-5,5},Filling->Axis]

μ=0;
σ=2;
norm=RandomVariate[NormalDistribution[μ,σ],1000];
Histogram[norm]

Усеченное  распределение

norm=NormalDistribution[1,2];
TD=TruncatedDistribution[{0,+∞},norm];
Plot[{PDF[norm,x],PDF[TD,x]},{x,-3,7},Filling->Axis]

μ=1;
σ=2;
norm1=RandomVariate[NormalDistribution[μ,σ],1000];
norm2=RandomVariate[TruncatedDistribution[{2,5},NormalDistribution[μ,σ]],1000];
Histogram[{norm1,norm2},{-10,10,0.25}]

Треугольное распределение

Clear[c,x]
Table[PDF[TriangularDistribution[{0,4},c],x],{c,1,3}]
Plot[%,{x,-1,5},Filling->Axis]

norm=RandomVariate[TriangularDistribution[{0,4},2],1000];
Histogram[norm]

PERT распределение

Clear[c,x]
Table[PDF[PERTDistribution[{0,4},c],x],{c,1,3}]
Plot[{%,PDF[PERTDistribution[{2,6},4,1],x]},{x,-1,8},Filling->Axis]

norm=RandomVariate[PERTDistribution[{0,4},2],1000];
Histogram[norm]

Равномерное распределение

Clear[c,x]
Table[PDF[UniformDistribution[{0,c}],x],{c,1,3}]
Plot[%,{x,-1,8},Filling->Axis]

μ=0;
λ=2;
norm=RandomVariate[UniformDistribution[{0,2}],1000];
Histogram[norm,{-1,3,0.25}]

Логонормальное распределение

Clear[c,x]
μ=1;
Table[PDF[LogNormalDistribution[μ,σ],x],{σ,1,3}]
Plot[%,{x,-1,8},PlotRange->{0,1},Filling->Axis]

PDF[LogNormalDistribution[0,1],x]
Plot[%,{x,-1,8},PlotRange->{0,1},Filling->Axis]

lognorm1=PDF[LogNormalDistribution[0,1],x]
t=4;
lognorm2= PDF[LogNormalDistribution[0,1],x-t]
Plot[{lognorm1, lognorm2},{x,-1,8},PlotRange->{0,1},Filling->Axis]

t=4
log1=RandomVariate[LogNormalDistribution[0,1],10000];
log2=RandomVariate[LogNormalDistribution[0,1],10000]+t;
Histogram[{log1,log2},{-1,10,0.25}]

Моделирование вывода лекарства на рынок

Стандартная модель

Clear["Global`*"]
c={6,6.05,6.1};
n={802000,967000,1132000};
v=c*n;
s=0.55*v;
vp=v-s;
op=0.15 vp;
cd=vp-op;
ng=0.32 cd;
cf=cd-ng;
NPV=TimeValue[Cashflow[Prepend[cf,-3400000]],0.10,0]
IRR=irr/.FindRoot[
TimeValue[Cashflow[Prepend[cf,-3400000]],irr,0]==0,
{irr,0.15}]

Модель методом монте-карло

Clear["Global`*"]
ts=SessionTime[];
Ftr[min_,cnt_,max_]:=RandomVariate[TriangularDistribution[{min,max},cnt]]
Fnorm[cnt_,sig_]:=RandomVariate[NormalDistribution[cnt,sig]]

tn=10;

c=Transpose[{
Array [Ftr[5.9,6,6.1]&,tn],
Array [Ftr[5.95,6.05,6.15]&,tn],
Array [Ftr[6.0,6.1,6.2]&,tn]
}];
n=Transpose[{
Array [Fnorm[802000,25000]&,tn],
Array [Fnorm[967000,30000]&,tn],
Array [Fnorm[1132000,25000]&,tn]
}];
v=c*n;
s=Array [Ftr[0.5,0.55,0.65]&,{tn,3}]*v;
vp=v-s;
op=Array [Fnorm[0.15,0.02]&,{tn,3}]*vp;
cd=vp-op;
ng=0.32 cd;
cf=cd-ng;
cff=Prepend[#,-3400000]&/@ cf;
NPV=TimeValue[Cashflow[#],0.10,0]&/@ cff
IRR=irr/.FindRoot[TimeValue[Cashflow[#],irr,0]==0,{irr,0.15}]&/@ cff
SessionTime[]-ts

Clear["Global`*"]
cff=Quiet[ TimeValue[Cashflow[{-3400000,cf[[i,1]],cf[[i,2]],cf[[i,3]]}],p,0]]
Quiet[  Solve[-3400000+((cf[[i,1]])/(1+p))+((cf[[i,2]])/(((1+p))^2))+((cf[[i,3]])/(((1+p))^3))==0,p,Reals]];
irr=p/.First[%]

ts=SessionTime[];
Ftr[min_,cnt_,max_]:=RandomVariate[TriangularDistribution[{min,max},cnt]]
Fnorm[cnt_,sig_]:=RandomVariate[NormalDistribution[cnt,sig]]

tn=10;
p=0.1;
c=Transpose[{
Array [Ftr[5.9,6,6.1]&,tn],
Array [Ftr[5.95,6.05,6.15]&,tn],
Array [Ftr[6.0,6.1,6.2]&,tn]
}];
n=Transpose[{
Array [Fnorm[802000,25000]&,tn],
Array [Fnorm[967000,30000]&,tn],
Array [Fnorm[1132000,25000]&,tn]
}];
v=c*n;
s=Array [Ftr[0.5,0.55,0.65]&,{tn,3}]*v;
vp=v-s;
op=Array [Fnorm[0.15,0.02]&,{tn,3}]*vp;
cd=vp-op;
ng=0.32 cd;
cf=cd-ng;
NPV=Table[cff,{i,1,tn}]
IRR=Table[irr,{i,1,tn}]
SessionTime[]-ts

Итог

Clear["Global`*"]
ts=SessionTime[];
cff=Quiet[ TimeValue[Cashflow[{-3400000,cf[[i,1]],cf[[i,2]],cf[[i,3]]}],p,0]];
Quiet[  Solve[-3400000+((cf[[i,1]])/(1+p))+((cf[[i,2]])/(((1+p))^2))+((cf[[i,3]])/(((1+p))^3))==0,p,Reals]];
irr=p/.First[%];
Ftr[min_,cnt_,max_]:=RandomVariate[TriangularDistribution[{min,max},cnt]]
Fnorm[cnt_,sig_]:=RandomVariate[NormalDistribution[cnt,sig]]
tn=10000;
p=0.1;
c=Transpose[{
Array [Ftr[5.9,6,6.1]&,tn],
Array [Ftr[5.95,6.05,6.15]&,tn],
Array [Ftr[6.0,6.1,6.2]&,tn]
}];
n=Transpose[{
Array [Fnorm[802000,25000]&,tn],
Array [Fnorm[967000,30000]&,tn],
Array [Fnorm[1132000,25000]&,tn]
}];
v=c*n;
s=Array [Ftr[0.5,0.55,0.65]&,{tn,3}]*v;
vp=v-s;
op=Array [Fnorm[0.15,0.02]&,{tn,3}]*vp;
cd=vp-op;
ng=0.32 cd;
cf=cd-ng;
NPV=Table[cff,{i,1,tn}];
IRR=Table[irr,{i,1,tn}];
SessionTime[]-ts

Print["----------------------------------------"]
Print[" "]
Print["NPV ",Mean[NPV],"     IRR ",Mean[IRR]]
Print["Вероятность отриц. NPV ",Length[Select[NPV,#≤ 0&]]/tn*1.]

Histogram[NPV,{-500000,800000,50000}]
Histogram[IRR,{0,0.3,0.01}]
DistributionChart[NPV]
BoxWhiskerChart[NPV,"Mean"]


Grid[{
{" ","Min","Mean","Max"},
{"NPV",Min[NPV],Mean[NPV],Max[NPV]},
{"IRR",Min[IRR],Mean[IRR],Max[IRR]}
},Frame->All]

Grid[{
{" ","P10","P50","P90"},
Flatten[{"NPV",Quantile[NPV,{1/10,1/2,9/10},{{1/2, 0}, {0, 1}}]}],
Flatten[{"IRR",Quantile[IRR,{1/10,1/2,9/10},{{1/2, 0}, {0, 1}}]}]
},Frame->All]

Grid[{
{" ","Q 1/4","Q 1/2","Q 3/4"},
Flatten[{"NPV",Quartiles[NPV]}],
Flatten[{"IRR",Quartiles[IRR]}]
},Frame->All]