18. Probabilistic System Assessment

18.1 Use of Statistics Packages

Needs["Statistics`Master`"]
Needs["Graphics`Graphics`"]

LogUniformDistribution::usage =
"LogUniformDistribution[min, max] represents the 
continuous loguniform distribution with endpoints 
min and max with min < max. Log[10, x] of the 
random variable x is uniformly
distributed on [Log[10, min], Log[10, max]]";

LogUniformDistribution/: 
Domain[LogUniformDistribution[min_:1, max_:10]] :=
        {min, max}

LogUniformDistribution/: 
Random[LogUniformDistribution[min_:1, max_:10]] :=
        10.0^Random[Real, {N[Log[10, min]], N[Log[10, max]]}]

LogUniformDistribution/: 
PDF[LogUniformDistribution[min_:1, max_:10], 
x_] := With[{result = 1/(x*Log[max/min])},
	If[NumberQ[N[result]],N[result],result]]

?LogUniformDistribution

q = LogUniformDistribution[3, 300]

PDF[q, x]

Plot[PDF[q, x], {x, 3, 300}];

logunidata = Table[Random[q], {i, 1, 1000}];

binned = BinCounts[logunidata, {3, 300, 25}]

plotted = BarChart[binned, PlotRange -> All];

18.2 An Environmental Example: Migration

T[L_, P_, H_, K_] := (3.1688*10^-8*L*P)/(H*K)

T[100, 0.38, 0.01, 3*10^-5]

18.2.1 Mathematica Implementation

TruncateRandom[dist_, a_, b_] := 
Module[{q},
(Label[start];
q = Random[dist];
If[q < a || q > b, Goto[start], q])]

Randomizer := 
{L -> 100, 
 P -> Random[UniformDistribution[0.35, 0.40]],
 H -> Random[LogUniformDistribution[0.003, 0.03]],
 K -> TruncateRandom[NormalDistribution[3 10^-5, 10^-5], 
 0.5 10^-5, 5.5 10^-5]}

T[L, P, H, K]

T[L, P, H, K] /. Randomizer

Tdata = Table[T[L, P, H, K] /. Randomizer, {i, 1, 1000}];

LocationReport[Tdata]

DispersionReport[Tdata]

Quartiles[Tdata]

{Min[Tdata], Max[Tdata]}

BinCounts[Log[10, Tdata], {-0.5, 2, 0.125}]

BarChart[%];

18.3 A Financial Example: Aggregation

Emms[Eoss_, Ecss_, noss_, ncss_] := 
  (Eoss*noss + Ecss*ncss)/(noss + ncss)

Emms[26.3, 21.1, 1, 0.5]

18.3.1 Mathematica Implementation of Probabilistic Aggregation

Randomizer := 
{noss -> 1, 
 Eoss -> Random[UniformDistribution[24, 30]],
 Ecss -> Random[LogUniformDistribution[5, 50]],
 ncss -> 0.5}

Emms[Eoss, Ecss, noss, ncss]

Emms[Eoss, Ecss, noss, ncss] /. Randomizer

Edata = Table[Emms[Eoss, Ecss, noss, ncss] /. Randomizer, 
{i, 1, 1000}];

LocationReport[Edata]

bindata = BinCounts[Edata, {15, 40, 2.5}]

bins = Table[16.25 + 2.5 (n - 1), {n, 1, 10}]

BarChart[bindata, PlotRange -> All, BarValues -> True,
BarLabels -> Map[ToString, bins], 
DefaultFont -> {"Courier", 6}];

18.4 A More Structured Approach to PSA

CleanSlate[];
Needs["Prosysas`"]

$ContextPath

?Prosysas`*

T[L_, P_, H_, K_] := (3.1688*10^-8*L*P)/(H*K)

ranset := 
MakeRandom[T[L, P, H, K], 
{{L, 100}, {P, Random[UniformDistribution[0.35, 0.40]]},
{H, Random[LogUniformDistribution[0.003,
 0.03]]}, {K, TruncateRandom[NormalDistribution[3 10^-5, 
 10^-5], 0.5 10^-5, 5.5 10^-5]}}]

ranset

Tdata = Table[ranset, {i, 1, 1000}];

18.5 Expectations and Volume Integration

prod[x_, y_, z_] := x y z

Randomizera := 
{x -> Random[UniformDistribution[0.0, 1.0]],
 y -> Random[UniformDistribution[0.0, 1.0]],
 z -> Random[UniformDistribution[0.0, 1.0]] 
}

prod[x, y, z] /. Randomizera

average[n_] := Mean[
 Table[prod[x, y, z] /. Randomizera, {i, 1, n}]]; 

nvalues = { 50, 50, 50, 250, 250, 250, 
	500, 500, 500, 750, 750, 750};

averages = Table[average[nvalues[[b]]], {b, 1, 12}];

PsaIntegralData = ListPlot[Transpose[{nvalues, averages}], 
PlotStyle -> {PointSize[0.01]},
PlotRange -> {0.09, 0.14},
AxesOrigin -> {0, 0.125}];

18.6 Wozniakowski's Method and Hammersley Points

a[k_Integer, j_Integer] := IntegerDigits[k, Prime[j]];

a[38, 2]

phi[k_, j_] := Module[{b, len, pj},
	pj = Prime[j]; 
	b = IntegerDigits[k, pj];
	len = Length[b];
	Sum[b[[len - i]]/(pj^(i+1)), {i, 0, len-1}]
	]

phi[38, 2]

1/3^4 + 1/3^3 + 0/3^2 + 2/3

Anzero[e_, d_] = Floor[(1/e)*(Log[1/e])^((d - 1)/2)];
Bnzero[e_, d_] = Floor[(1/e)*(Log[1/e])^(d - 1)];
Cnzero[e_, d_] = Floor[(1/e)*(Log[1/e])^d];

M[d_, n_] := (Product[Prime[i], 
		{i, 1, d-1}])^Ceiling[N[Log[n]]]

M[3, n]

X[d_, n_] := Table[
	Which[m < 2, 1 - k/n, 
	m > 1, 1 - phi[k, m-1]],
	{k, 1, n}, {m, 1, d}]

n = Bnzero[0.1, 3]

Short[S = X[3, 53],4]

Show[Graphics3D[Map[Point, S]]];

18.7 Efficient Integration

g[r_List] := r[[1]]*r[[2]]*r[[3]]

WozIntegral = N[Apply[Plus, Map[g, S]]/n]

WozInt[f_, e_, d_] := 

Module[{n, S},
	n = Bnzero[e, d];
	S = X[d, n];
	{n, N[Apply[Plus, Map[f, S]]/n]}];

WozInt[g, 0.06, 3]

wozdata = Table[WozInt[g, u, 3], {u, 0.02, 0.2, 0.02}]

We can compare these values with the sampled estimates established earlier:

ListPlot[wozdata, PlotJoined -> True, 
PlotRange -> {0.09, 0.15},
AxesOrigin -> {0, 0.125},
PlotStyle -> {Thickness[0.001]},
Evaluate[Epilog -> {PointSize[0.01], 
Map[Point, Transpose[{nvalues, averages}]]}]   ];

18.9 Appendix: Prosysas Package

BeginPackage["Prosysas`","Statistics`Master`",
             "Graphics`Graphics`"]

Prosysas::usage = "Prosysas.m is a prototype package for
probabilistic system assessment."

(* Usage Messages for exported  objects *)

LogUniformDistribution::usage = 
 "LogUniformDistribution[min, max] represents the 
  continuous loguniform distribution with endpoints 
  min and max with min < max. Log[10, x] of the 
  random variable x is uniformly 
  distributed on [Log[10, min], Log[10, max]]"

TruncateRandom::usage = "TruncateRandom[dist, a, b]
produces a random sampling of dist truncated between a and b"

MakeRandom::usage = "MakeRandom[fn, {{x1, dist1}, {x2, dist2} ...}]
 produces a random sampling f with each variable xi replaced by
 a random sampling from disti"

Begin["`Private`"]

(* Function Definitions *)

LogUniformDistribution/:
Domain[LogUniformDistribution[min_:1, max_:10]] := 
       {min, max}

LogUniformDistribution/: 
Random[LogUniformDistribution[min_:1, max_:10]] :=
10.0^Random[Real, {N[Log[10, min]], N[Log[10, max]]}]

LogUniformDistribution/: 
PDF[LogUniformDistribution[min_:1, max_:10], x_] := 
With[{result = 1/(x*Log[max/min])}, 
If[NumberQ[N[result]],N[result],result]]

TruncateRandom[dist_, a_, b_] := 
Module[{q},
(Label[start];
q = Random[dist];
If[q < a || q > b, Goto[start], q])]

 

MakeRandom[fn_, vbles_List] := 
fn /. Map[((#[[1]] -> #[[2]])&), vbles] 

End[]

EndPackage[]