Archive

Posts Tagged ‘Mathematica’

Probability manipulations with Mathematica…

There comes a time that a statistician needs to do some ananytic calculations. There more than a bunch of tools to use but I usually prefer Mathematica or Maple. Today, I’m gonna use Mathematica to do a simple exhibition.

Let’s set this example upon the  U(2 \theta _1-\theta _2\leq x\leq 2 \theta _1+\theta _2)   distribution.

pfun = PDF[UniformDistribution[{2*Subscript[θ, 1] - Subscript[θ, 2],
2*Subscript[θ, 1] + Subscript[θ, 2]}], x]

\begin{cases}  \frac{1}{2 \theta _2} & 2 \theta _1-\theta _2\leq x\leq 2 \theta _1+\theta _2 \\  0 & \text{True}  \end{cases}

One of the most intensive calculations is the characteristic function (eq. the moment generating function). This is straightforward to derive.

cfun=CharacteristicFunction[UniformDistribution[
{2*Subscript[θ, 1]-Subscript[θ, 2],2*Subscript[θ, 1]+Subscript[θ, 2]}],x]

-\frac{i \left(-e^{i x \left(2 \theta _1-\theta _2\right)}+e^{i x \left(2 \theta _1+\theta _2\right)}\right)}{2 x \theta _2} .

The Table[] command calculates for us the raw moments for our distribution.

Table[Limit[D[cfun, {x, n}], x -> 0]/I^n, {n, 4}]

\left\{2 \theta _1,\frac{1}{3} \left(12 \theta _1^2+\theta _2^2\right),2 \theta _1 \left(4 \theta _1^2+\theta _2^2\right),16 \theta _1^4+8 \theta _1^2 \theta _2^2+\frac{\theta _2^4}{5}\right\} .

Calculate the sample statistics.

T=List[8.23,6.9,1.05,4.8,2.03,6.95];
{Mean[T],Variance[T]}

\{4.99333,8.46171\} .

Now, we can use a simple moment matching technique to get estimates for the parameters.

Solve[{Mean[T]-2*Subscript[θ, 1]==0,-(2*Subscript[θ, 1])^2+
1/3 (12 Subscript[θ, 1]^2+\!\*SubsuperscriptBox[\(θ\), \(2\), \(2\)])-
Variance[T]==0},{Subscript[θ, 2],Subscript[θ, 1]}]

\left\{\left\{\theta _1\to 2.49667,\theta _2\to -5.03836\right\},\left\{\theta _1\to 2.49667,\theta _2\to 5.03836\right\}\right\} .

Check the true value for the \theta _2.

Reduce[2 Subscript[θ, 1]-Subscript[θ, 2]<=2 Subscript[θ, 1]+Subscript[θ, 2],
Subscript[θ, 2]]

\theta _1\in \text{Reals}\&\&\theta _2\geq 0 .

Then, \left\{\left\{\theta _1\to 2.49667,\theta _2\to  5.03836\right\}\right\} .

May the rule be with you ;)

Actuaries are not famous for their sense of humour. Sometimes, they exhibit a remarkable sense of it though!

A second year undergrad student is familiar with the following formula,

EX=\int_{0}^{\infty }{\left( 1-F(x) \right)dx},

where X is a random variable satisfying X\ge 0 with probability 1. Most call 1-F(x)\, the right tail and others the survival function of X,denoted by S(x)\ . This is a really convenient way to calculate expected values instead of using integration by parts.

Actuaries find this quite handy especially when dealing with deductibles that shift the variable or stop limits that truncate the variable from above. The first is usually treated whether the latter makes no real difference as the survival function will be 0 over the stop limit.

Where’s the joke that actuaries tipped in? Well, they called this the Darth Vader Rule given that Darth Vader was a true survivor…

Below you can see the expected value of a Gamma(2,5) random variable.

To get the plot it’s easy with Mathematica for instance.

Plot[CDF[GammaDistribution[2, 5], x], {x, 0, 25}, PlotRange -> {{0, 29}, {0, 1}}, Filling -> Top]

Next calculate the expected value.

Integrate[1 – CDF[GammaDistribution[2, 5], z], {z, 0, 25}] // N
9.76417

This is alpha+…

It’s great. It could be the new google for us I can see this coming…

Wolfram’s Alpha goes NYSE : Quote NBG (National Bank of Greece)

Categories: infos Tags: , , , , ,

A merry wolfram xmas!

December 14, 2009 Leave a comment

Christmas coming & time for fun!  Wolfram’s demonstrations give you a sense of the holiday season along some nice demonstrations. I particulary like the “Ornamental Holiday Decoration

Manipulate[
 Module[{level0, level1, level2},
 level0 = C[spikey, 0];
 level1 = 
 Flatten[daughterPolyhedra[C[spikey, 0], {d1, ω1, s1}, ρ s1]];
 level2 = 
 Flatten[daughterPolyhedra[#, {d2, ω2, s2}, ρ  s1 s2] & /@ 
 Cases[level1, _C]];
 Graphics3D[{EdgeForm[],
 {Red, egc[level0]}, gg1 = {col1, egc[level1]}, 
 gg2 = {col2, ControlActive[{}, egc[level2]]},
 Directive[GrayLevel[0.2], Specularity[colc, 12]], 
 ecc[{level1, level2}]}, Boxed -> False, 
 ImageSize -> {400, 400}]],
 "layer 1:",
 {{d1, 1.5, "distance"}, -3, 3, ImageSize -> Tiny},
 {{ω1, 0, "rotation"}, -Pi, Pi, ImageSize -> Tiny},
 {{s1, 1/2, "size"}, 0, 1, ImageSize -> Tiny},
 {{col1, Yellow, "color"}, Blue, ControlType -> None},
 Delimiter, 
 "layer 2:",
 {{d2, 1, "distance"}, -3, 3, ImageSize -> Tiny},
 {{ω2, 0, "rotation"}, -Pi, Pi, ImageSize -> Tiny},
 {{s2, 1/2, "size"}, 0, 1, ImageSize -> Tiny},
 {{col2, Green, "color"}, Green, ControlType -> None},
 Delimiter,
 "connectors:",
 {{ρ, 0.3, "radius"}, 0, 1, ImageSize -> Tiny},
 {{colc, Brown, "color"}, Yellow, ControlType -> None}, 
 AutorunSequencing -> {1, 3, 5, 7},
 Initialization :> {
 spikey = 
 MapAt[Developer`ToPackedArray, 
 MapAt[Developer`ToPackedArray, N[PolyhedronData["Spikey"]][[1]], 
 1], {2, 1}]; 
 mirrorRotateAndShift[gc_GraphicsComplex, 
 n_, {distance_, angle_, size_}, ρ_] := 
 With[{aux = 
 mirrorRotateAndShiftCF[gc[[1]], gc[[1, n]], distance, angle, 
 size]}, {C[GraphicsComplex[aux, gc[[2]]], n], 
 Cylinder[{gc[[1, n]], aux[[n]]}, ρ]}]; 
 mirrorRotateAndShiftCF = 
 Compile[{{vertices, _Real, 2}, {rPoint, _Real, 1}, distance, 
 angle, size}, Read more...

Web analytics revisited : The Mathematica way!

December 12, 2009 Leave a comment

I found out something really interesting. Men do it with keyboards (not clicks!) and so does Patrick Collison using Mathematica to analyse web data he got from Google Analytics.

See : Hacking for fun and profit with Mathematica and the Google Analytics API [link]

Simulating the maximum…

Studying the distributional properties of the maximum is gaining more of attention nowadays given the relative failure of the conventional tools the quants used before the recession.

Let’s suppose that an insurance company faces the challenge to cover claims incurring in random time and size. The intensity of the process is \lambda (t)=100t(1-t), \lambda (0)=25 and the process itself is (surprise!!!) Poisson. The claims are iid from a lognormal distribution with parameters \mu=0.75t and \sigma=1, where t is the intermediate time between occurence i and i-1.

Some Mathematica code from my Graduate Simulation class is the following.

First we set up the intersity function and the generic vectors recoring Claims (Claimsv) and maximums (Maxv)

λ0 = 25; t0 = 1;
λ[t_] := 100 (t – t^2); Claimsv = {}; Maxv = {};

The main program is the following

Do[S = {}; t = 0; n = 0;
X = -Log[Random[]]/λ0; t = t + X;
While[t < t0, If[Random[] < λ[t]/λ0, n = n + 1; S = Append[S, t]];
X = -Log[Random[]]/λ0; t = t + X];
Claims = 0; S[[0]] = 0; Maxx = 0; Do[rr = Exp[
Random[NormalDistribution[0.75*(S[[i]] – S[[i – 1]]), 1]]];
If[rr > Maxx, Maxx = rr];
Claims = Claims + rr, {i, 1, n}] AppendTo[Claimsv, Claims];
AppendTo[Maxv, Maxx];, {10000}]

And plotting the histograms of the simulated distribution of the Claims and the maximum we nearly finished…

<< Histograms`

Histogram[Claims.v]
claimsv Histogram[Max.v]

maxv

to be continued…