Teorie spolehlivosti

Priklad porovnani serioveho a paralelniho systemu

(Subscript[X, 1]),...,(Subscript[X, n]) iid s rozdelenim exponencialnim s parametrem θ=1/λ, tj. s hustotou

Clear[pexp,θ];
d1=ExponentialDistribution[(1/θ)];
pexp[x_]:=PDF[d1,x]
pexp[x]

Stredni doba do poruchy kazdeho prvku je

Mean[d1]

Stredni doba do poruchy zalohovaneho systemu s nezatizenymi zalohami E (X^((n)))=n θ.

Doba do poruchy zalohovaneho systemu se zatizenymi zalohami ma distribucni funkci

P{(Subscript[X, max])<x}=P{∀i=1,...n (Subscript[X, i])<x}=(F^n)(x)=

Clear[dfFn];
dfFn[x_]:=(((1-Exp[-(x/θ)]))^n)

Stredni hodnota E (Subscript[X, max]) je tedy

eXmax=Assuming[θ>0&&n>0&&n∈Integers,(∫_0^∞)(1-dfFn[x])\[DifferentialD]x]

Expand[dfFn[x]]

(1-(((1-(E^(-(x/θ)))))^n))/.(1-(E^(-(x/θ))))->y

Solve[y==(1-(E^(-(x/θ)))),x]

D[x/.{{x->-θ Log[1-y]}},y]

EXmax=θ Integrate[(((1-(y^n)))/(1-y)),{y,0,1},Assumptions->{n∈Integers&&n>0}]

Series[HarmonicNumber[n],{n,∞,1}]//Normal//FullSimplify//PowerExpand
EulerGamma//N

Z toho vyplyva, ze Stredni hodnota E (Subscript[X, max]) je priblizne C+log n. Podil strednich hodnot je priblizne

ListPlot[Table[((Log[n]+ EulerGamma)/n),{n,30}],Joined->True,AxesOrigin->{0,0}]

Pravdepodobnostni rozdeleni

Exponencialni rozdeleni Exp(θ)

Clear[pexp,θ];
d1=ExponentialDistribution[(1/θ)];
pexp[x_]:=PDF[d1,x]
pexp[x]

distribucni funkci a funkci spolehlivosti

Clear[rexp];
CDF[d1,x]
rexp=1-CDF[d1,x]//Simplify
Plot[rexp/.θ->1,{x,0,5}]

Intensita poruch je

Clear[rexp];
rexp[x_]:=((pexp[x])/(E^(-(x/θ))))

rexp[y]//Simplify

Momenty

ExpectedValue[(x^k),d1,x,Assumptions->{k∈Integers&&k>0}]

{Mean[d1],Variance[d1]}

{Skewness[d1],Kurtosis[d1]}

Variacni koeficient (Sqrt[var X])/E X

((StandardDeviation[d1])/(Mean[d1]))

Rozdeleni (Subscript[N, t])

PDF[GammaDistribution[p,a],y]

Konvoluce rozdeleni gamma

Integrate[PDF[GammaDistribution[p1,a],y]PDF[GammaDistribution[p2,a],x-y],{y,0,x},Assumptions->{p1>0,p2>0,a>0,x>0}]//Simplify//PowerExpand

CharacteristicFunction[GammaDistribution[p,a],y]

P{(Subscript[N, t])=k}=

Assuming[k∈Integers&&k>0&&θ>0,CDF[GammaDistribution[k,θ],t]-CDF[GammaDistribution[k+1,θ],t]]//FunctionExpand//FullSimplify

Integrate[f[x] f[y],{x,-∞,∞},{y,x,∞}]

(∑_(j=1)^k)(1/(((n-j+1))^2))//FunctionExpand

Limit[(∑_(j=1)^k)(1/(((n-j+1))^2)),n->∞,Assumptions->{k>0,k∈Integers}]

Limit[(∑_(j=1)^n)(1/(((n-j+1))^2)),n->∞]

{Limit[(∑_(j=1)^n)(1/(j^2)),n->∞],(∑_(j=1)^∞)(1/(j^2))}

Weibullovo rozdeleni W(θ, β)

CDF[d2=WeibullDistribution[β,θ],x]

r[x_]=((PDF[d2,x])/(1-CDF[d2,x]))

Plot[Table[r[x]/.{θ->1,β->β1},{β1,{0.5,1,1.5,2,2.5}}],{x,0,1}]

Assuming[k∈Integers&&k>0,ExpectedValue[(x^k),d2,x]]

{Mean[d2],Variance[d2]}

{θ Gamma[1+(1/β)],(θ^2) (-((Gamma[1+(1/β)])^2)+Gamma[1+(2/β)])}

Variacni koeficient

vk=((StandardDeviation[d2])/(Mean[d2]))//Simplify

nezavisi na θ !!!!

Plot[vk,{β,0.1,3},AxesOrigin->{0,0}]

Priklad (Odhad parametru β momentovou metodou)

SeedRandom[13];
{Mean[WeibullDistribution[1.5,1]],StandardDeviation[WeibullDistribution[1.5,1]]}
randomsample=RandomReal[WeibullDistribution[1.5,1],50];
{ex=Mean[randomsample],sd=StandardDeviation[randomsample]}
fr=FindRoot[((Sqrt[-((Gamma[1+(1/β)])^2)+Gamma[((2+β)/β)]])/(Gamma[1+(1/β)]))==(sd/ex),{β,1}]
β/.fr
odhadθ= (ex/(Gamma[1+(1/β)]))/.fr

Clear[simulace];
simulace[θ_,β_,n_,nN_]:=Module[{pomocna,randomsample,fr,ex,sd,odhadθ,odhadβ,b},
pomocna=Table[randomsample=RandomReal[WeibullDistribution[β,θ],n];{ex=Mean[randomsample],sd=StandardDeviation[randomsample]};
fr=FindRoot[((Sqrt[-((Gamma[1+(1/b)])^2)+Gamma[((2+b)/b)]])/(Gamma[1+(1/b)]))==(sd/ex),{b,1}];
{odhadβ=b/.fr,odhadθ= (ex/(Gamma[1+(1/b)]))/.fr},{nN}]
								]

SeedRandom[13];
sim50a1000=simulace[1,1.5,20,25000];
{odhadyβ,odhadyθ}=Transpose[sim50a1000];
{exβ=Mean[odhadyβ],exθ=Mean[odhadyθ]}
{sdβ=StandardDeviation[odhadyβ],sdθ=StandardDeviation[odhadyθ]}
{Histogram[odhadyβ,Automatic,"ProbabilityDensity",Epilog->First@Plot[PDF[NormalDistribution[exβ,sdβ],x],{x,Min[odhadyβ],Max[odhadyβ]},PlotStyle->Red]],
Histogram[odhadyθ,Automatic,"ProbabilityDensity",Epilog->First@Plot[PDF[NormalDistribution[exθ,sdθ],x],{x,Min[odhadyθ],Max[odhadyθ]},PlotStyle->Red]]}

SeedRandom[13];
sim50a1000=simulace[1,1.5,50,10000];
{odhadyβ,odhadyθ}=Transpose[sim50a1000];
{exβ=Mean[odhadyβ],exθ=Mean[odhadyθ]}
{sdβ=StandardDeviation[odhadyβ],sdθ=StandardDeviation[odhadyθ]}
{Histogram[odhadyβ,Automatic,"ProbabilityDensity",Epilog->First@Plot[PDF[NormalDistribution[exβ,sdβ],x],{x,Min[odhadyβ],Max[odhadyβ]},PlotStyle->Red]],
Histogram[odhadyθ,Automatic,"ProbabilityDensity",Epilog->First@Plot[PDF[NormalDistribution[exθ,sdθ],x],{x,Min[odhadyθ],Max[odhadyθ]},PlotStyle->Red]]}

SeedRandom[13];
sim50a1000=simulace[1,1.5,100,5000];
{odhadyβ,odhadyθ}=Transpose[sim50a1000];
{exβ=Mean[odhadyβ],exθ=Mean[odhadyθ]}
{sdβ=StandardDeviation[odhadyβ],sdθ=StandardDeviation[odhadyθ]}
{Histogram[odhadyβ,Automatic,"ProbabilityDensity",Epilog->First@Plot[PDF[NormalDistribution[exβ,sdβ],x],{x,Min[odhadyβ],Max[odhadyβ]},PlotStyle->Red]],
Histogram[odhadyθ,Automatic,"ProbabilityDensity",Epilog->First@Plot[PDF[NormalDistribution[exθ,sdθ],x],{x,Min[odhadyθ],Max[odhadyθ]},PlotStyle->Red]]}

Intensita poruch

Manipulate[Plot[(β/θ)((((x/θ)))^(β-1)),{x,0.01,10}],{{θ,1},0.5,10},{{β,1.5},0.3,3}]

Plot[Table[(β/θ)((((x/θ)))^(β-1))/.θ->1,{β,{0.5,1,1.5,2,2.5}}],{x,0.01,2}]

Rayleighovo rozdeleni

{mr=Mean[d2/.{β->2,θ->σ(Sqrt[2])}],sdr=StandardDeviation[d2/.{β->2,θ->σ(Sqrt[2])}],(sdr/mr)//Simplify,(sdr/mr)//N}

Normalni rozdeleni N(μ,(σ^2))

d3=NormalDistribution[μ,σ]

Distribucni funkce N(μ,(σ^2)) je Φ(((x-μ)/σ))

CDF[NormalDistribution[μ,σ],x]

Pst, ze nah. vel. dude mensi nez μ-3σ?

{CDF[d3,μ-3σ],CDF[d3,μ-3σ]//N}

Pravidlo 2σ: P(|N(μ,σ) - μ| < 2σ)

CDF[d3,μ+2σ]-CDF[d3,μ-2σ]//N

Pravidlo 3σ: P(|N(μ,σ) - μ| < 3σ)

CDF[d3,μ+3σ]-CDF[d3,μ-3σ]//N

Hustota useknuteho normalniho rozdeleni

Clear[fr];
ftruncated[x_,μ_,σ_]:=(1/(σ CDF[NormalDistribution[0,1],μ/σ]))PDF[NormalDistribution[0,1],((x-μ)/σ)]

(Integrate[ftruncated[x,μ,σ],{x,0,∞}]/.If[a_,b_,c___]->b)//FullSimplify

Integrate[ftruncated[x,μ,σ],{x,0,∞},Assumptions->σ>0]

Plot[Table[ftruncated[x,μ,1],{μ,{0,1,3,5}}],{x,0,5}]

Stredni hodnota

meantruncated=Integrate[x ftruncated[x,μ,σ],{x,0,∞},Assumptions->σ>0]//Simplify

Table[{mm,"    "meantruncated/.{μ->mm,σ->1}}//N,{mm,0.5,5,0.5}]//TableForm

Intensita poruch N(μ,σ)

PDF[NormalDistribution[μ,σ],x]

Manipulate[Plot[{((x-μ)/(σ^2)),PDF[NormalDistribution[μ,σ],x]/(1-CDF[NormalDistribution[μ,σ],x])},{x,0,5}],{μ,0,5},{σ,0.1,3}]

PDF[NormalDistribution[μ,σ],x]/.{μ->1,σ->0.5,x->5}

(1-CDF[NormalDistribution[μ,σ],x])/.{μ->1,σ->0.5,x->5}

ser1=Assuming[σ>0,Series[PDF[NormalDistribution[μ,σ],x]/(1-CDF[NormalDistribution[μ,σ],x]),{x,∞,3}]]//Simplify

Priklad (prvek, poruchy nahle i postupne -> seriovy system) se spolehlivosti

Clear[Φ];
Φ[x_]=CDF[NormalDistribution[0,1],x]

Odpovidajici spolehlivost systemu je

rR[x_,θ_,μ_,σ_]:=Exp[((-x)/θ)]Φ[((μ-x)/σ)]

expected=Assuming[θ>0&&σ>0,(∫_0^∞)rR[x,θ,μ,σ]\[DifferentialD]x]

Konec 3. 11. 2009

Log-normalni rozdeleni LN(μ,(σ^2))

PDF[d4=LogNormalDistribution[μ,σ],x]

Momenty

{Mean[d4],Variance[d4],Skewness[d4],Kurtosis[d4]}

Plot[-3+3 (E^(2 (σ^2)))+2 (E^(3 (σ^2)))+(E^(4 (σ^2))),{σ,0,2}]

Manipulate[Plot[PDF[d4,x]/.{μ->μ1,σ->σ1},{x,0.5,4},PlotRange->All,AxesOrigin->{0,0}],{μ1,0,2},{σ1,0.1,2}]

Intensita poruch

Manipulate[Plot[PDF[d4,x]/(1-CDF[d4,x])/.{μ->μ1,σ->σ1},{x,0.5,4},PlotRange->All,AxesOrigin->{0,0}],{μ1,0,2},{σ1,0.1,2}]

Gamma rozdeleni Gamma(a,p)

PDF[d5=GammaDistribution[p,1/a],x]

{Mean[d5],Variance[d5],Skewness[d5],Kurtosis[d5]}

Prevracena hodnota intensity poruch:

Assuming[p>1&&a>0,(((∫_x^∞)(t^(p-1))(E^(-a t))\[DifferentialD]t)/((x^(p-1))(E^(-a x))))]

?Gamma

expr=Assuming[x>0&&p>0&&a>0,(((∫_x^∞)(t^(p-1))(E^(-a t))\[DifferentialD]t)/((x^(p-1))(E^(-a x))))]

dexpr=D[expr,x]//FullSimplify

Reduce[dexpr<0&&a>0&&p>0,{x,a,p},Reals]

Grafy intensity poruch

Table[(((x^(p-1))(E^(-x)))/(Assuming[x>0,(∫_x^∞)(t^(p-1))(E^(- t))\[DifferentialD]t])),{p,{0.5,1,2,3,5}}]//Evaluate

Plot[Table[(((x^(p-1))(E^(-x)))/(Assuming[x>0,(∫_x^∞)(t^(p-1))(E^(- t))\[DifferentialD]t])),{p,{0.5,1,2,3,5}}]//Evaluate,{x,0.001,5}]

Plot[(((x^(p-1))(E^(-x)))/((∫_x^∞)(t^(p-1))(E^(- t))\[DifferentialD]t))/.p->0.5,{x,0.001,5}]

CharacteristicFunction[d5,z]

CharacteristicFunction[GammaDistribution[p1,1/a],z]*CharacteristicFunction[GammaDistribution[p2,1/a],z]

Souvislost rozdeleni gamma s rozdelenim (χ^2)

(χ^2)(1)=(N^2)(0,1)

ma hustotu

PDF[ChiSquareDistribution[1],x]

PDF[ChiSquareDistribution[n],x]

To je hustota

PDF[GammaDistribution[(n/2),2],x]

⟹ (χ^2)(n)=Gamma(1/2,n/2)

Paretovo rozdeleni

{CDF[ParetoDistribution[k,α],x],PDF[ParetoDistribution[k,α],x]}

{Mean[ParetoDistribution[k,α]],Variance[ParetoDistribution[k,α]]}

Zobecnene Paretovo rozdeleni

Definition
The probability density function for the generalized Pareto distribution with shape parameter k ≠ 0, scale parameter ó, and threshold parameter è, is

for è < x, when k > 0, or for è < x < -ó/k when k < 0.
In the limit for k = 0, the density is

for è < x.
If k = 0 and è = 0, the generalized Pareto distribution is equivalent to the exponential distribution. If k > 0 and è = ó/k, the generalized Pareto distribution is equivalent to the Pareto distribution.

The density

Clear[pdfGPD,k,σ,θ];
pdfGPD[x_,k_,σ_,θ_]:=(1/σ)(((1+k((x-θ)/σ)))^(-1-1/k))

Manipulate[Plot[pdfGPD[x,k,σ,θ],{x,θ,10}],{θ,-1,3},{k,0.1,0.5},{σ,1,3}]

Case k>0

Indefinite integral

ix=Integrate[pdfGPD[x,k,σ,θ],x,Assumptions->σ>0&&k>0&&θ∈Reals]//PowerExpand//Simplify

Really, it is a density:

Density

Limit[ix,x->∞,Assumptions->σ>0&&k>0&&θ∈Reals]-(ix/.x->θ)

This is a mystery, however; why the additional assumption?

Integrate[pdfGPD[x,k,σ,θ],{x,θ,∞},Assumptions->σ>0&&k>0&&θ∈Reals]

Distribution function

ix-(ix/.x->θ)

Clear[cdfGPD];
cdfGPD[x_,k_,σ_,θ_]:=1-(((1+k((x-θ)/σ)))^(-1/k))

and the corresponding reliability function

Quantile function

Solve[1-(((1+k((x-θ)/σ)))^(-1/k))==y,x]//Simplify

Clear[quantileGPD];
quantileGPD[y_,k_,σ_,θ_]=((k θ+(-1+(((1-y))^(-k))) σ)/k)

Random numbers from GPD

tab1=Map[quantileGPD[#,1/2,1,1]&,RandomReal[{0,1},10000]];
Histogram[tab1]
{Min[tab1],Max[tab1]}

Moments

momentGPD[n_,k_,σ_,θ_]=Module[{indefiniteintegral},indefiniteintegral=Assuming[σ>0&&k>0&&θ∈Reals&&n<1/k,∫(x^n)pdfGPD[x,k,σ,θ]\[DifferentialD]x]//FullSimplify//PowerExpand;
Limit[indefiniteintegral,x->∞,Assumptions->σ>0&&k>0&&θ∈Reals&&n<1/k]-(indefiniteintegral/.x->θ)//FunctionExpand//Simplify
		]

m123=Table[momentGPD[n,k,σ,θ]/.{k->1/4},{n,3}]//Simplify//PowerExpand

(* variance *)
(θ^2)+((8 θ σ)/3)+((16 (σ^2))/3)-(((θ+((4 σ)/3)))^2)//Simplify

MLE

The likelihood

Log[pdfGPD[x,k,σ,θ]]//PowerExpand//Simplify

lik[x_,k_,σ_,θ_]=-(((1+k) Log[1+((k (x-θ))/σ)])/k)-Log[σ]

derlikk=D[lik[x,k,σ,θ],k]//FullSimplify

derlikσ=D[lik[x,k,σ,θ],σ]//FullSimplify

observations=Map[quantileGPD[#,1/2,1,1]&,RandomReal[{0,1},50]];
estimateθ=Min[observations]
sumk=Map[(derlikk/.θ->estimateθ)/.x->#&,observations]//Simplify;
sumσ=Map[(derlikσ/.θ->estimateθ)/.x->#&,observations]//Simplify;
totalk=Total[sumk];
totalσ=Total[sumσ];
FindRoot[{totalk==0,totalσ==0},{k,1/2},{σ,1}]

xx=Array[x,{3}];

sumk=Map[derlikk/.x->#&,xx]//Simplify

sumσ=Map[derlikσ/.x->#&,xx]//Simplify

Exponential distribution with location and scale parameters

Introduction

The support of the distribution is (θ,∞).

(∫_θ^x)(1/σ)Exp[-((y-θ)/σ)]\[DifferentialD]y

Clear[σ,θ,n];
pdfExp[x_,σ_,θ_]:=(1/σ)Exp[-((x-θ)/σ)]
cdfExp[x_,σ_,θ_]:=1-Exp[-((x-θ)/σ)]
momentExp[n_,σ_,θ_]=Assuming[θ>0&&σ>0&&θ∈Reals,(∫_θ^∞)(x^n)(1/σ)Exp[-((x-θ)/σ)]\[DifferentialD]x]

First three moments

Clear[n,σ,θ];
Table[momentExp[n,σ,θ],{n,3}]//FunctionExpand//Simplify
(* rozptyl *)(θ^2)+2 θ σ+2 (σ^2)-(((θ+σ))^2)//Simplify

Quantile function

Solve[1-Exp[-((x-θ)/σ)]==y,x]

quantileExp[y_,σ_,θ_]:=θ -σ Log[1-y]

tab1=Map[quantileExp[#,1,1]&,RandomReal[{0,1},1000]];

Histogram[tab1]

MLE

Parametr θ se odhadne maximem ze vsech pozorovani.

lik[x_,σ_,θ_]=Log[pdfExp[x,σ,θ]]//PowerExpand
derlikσ=D[lik[x,σ,θ],σ]

xx=Array[x,{3}];

sumσ=Map[derlikσ/.x->#&,xx]//Simplify
totalσ=Total[sumσ]
Solve[totalσ==0,σ]

Z toho je videt, ze MLE parametru σ je

σ=(1/n)(∑_(i==1)^n)(Subscript[X, i])-θ

Example

quantileExp[y_,σ_,θ_]:=θ -σ Log[1-y]
tab1=Map[quantileExp[#,1,1]&,RandomReal[{0,1},100]];

Print["Odhad θ = ",odhadθ=Min[tab1]," Odhad σ = ",odhadσ=Mean[tab1]-odhadθ]

quantileExp[y_,σ_,θ_]:=θ -σ Log[1-y]
nN=100;n=1000;σ=1;θ=1;
(* N vyberu o rozsahu n v nasledujici tabulce *)tabulkavelka=Table[Map[quantileExp[#,1,1]&,RandomReal[{0,1},n]],{nN}];
Dimensions[tabulkavelka]
odhadyθ=Map[Min,tabulkavelka];
nestranneodhadyθ=odhadyθ-(1/n);
odhadyσ=Map[Mean,tabulkavelka]-odhadyθ;
Histogram[odhadyθ]
Histogram[odhadyσ]

Ciselne charakteristiky

{Mean[odhadyθ],Mean[odhadyσ]}
{StandardDeviation[odhadyθ],StandardDeviation[odhadyσ]}

Interval spolehlivosti pro σ

Needs["HypothesisTesting`"]
MeanCI[odhadyσ]

Pomocne vypocty

Module[{n}, Assuming[n>0,(∫_0^∞)x Exp[-n x]\[DifferentialD]x]]

Least Squares for Location an Scale Parameters

Order statistics

Distribution function of the k-th order statistics

dfc[x_,k_,fF_]:=Sum[Binomial[n,j] ((fF[x])^j)(((1-fF[x]))^(n-j)),{j,k,n}]

dfcexp=Assuming[k>0&&n>0&&k∈Integers&&n∈Integers,Simplify[dfc[x,k,CDF[ExponentialDistribution[1],#]&]]]//FunctionExpand

relexp=1-dfcexp//FullSimplify

expexp=Integrate[relexp,{x,0,∞},Assumptions->k>0&&n>0&&k∈Integers&&n∈Integers]

relexp/.{n-> 100,k->50}

?MaxRecursion

NIntegrate[relexp/.{n-> 100,k->50},{x,0,∞},MaxRecursion->20]

Integrate[relexp/.{n-> 100,k->1},{x,0,∞}]

NIntegrate[relexp/.{n-> 100,k->100},{x,0,∞}]

Useful referencess

http://www.mathwave.com/articles/generalized_pareto_distribution.html

Hurt, J.: Asymptotic expansions ... (Truncated exponential distribution)

http://library.wolfram.com/infocenter/Conferences/7454/

Estimate of the failure rate of an Exponential distribution under random censorship

PDF[GammaDistribution[α,β],x]

(∑_(j==1)^n)(Subscript[W, j]) possesses

Clear[n,θ,γ];
distrdg=GammaDistribution[ n,(θ/(1+γ))];

{Mean[distrdg],Variance[distrdg]}

expected1z=Assuming[n>1,ExpectedValue[1/z,distrdg,z]]

It follows that E(OverHat[λ ])= E,Bi(n,1/(1+γ))·expectedz, i. .e.

Mean[BinomialDistribution[n,1/(1+γ)]] expected1z

Thus ((n-1)/n)(OverHat[λ]) is an unbiased estimate of λ.

Probit and Logit Analysis

ProbitModelFit

Heart Transplant Data

Odkaz zde