In [1]:

    

from sympy import *; init_session()


    

    




IPython console for SymPy 1.0 (Python 2.7.12-64-bit) (ground types: python)

These commands were executed:
>>> from __future__ import division
>>> from sympy import *
>>> x, y, z, t = symbols('x y z t')
>>> k, m, n = symbols('k m n', integer=True)
>>> f, g, h = symbols('f g h', cls=Function)
>>> init_printing()

Documentation can be found at http://docs.sympy.org/1.0/

In [2]:

    

import random as rand
import numpy as np


    

In [3]:

    

from applpy import *


    

In [4]:

    

# Example 1


    

In [5]:

    

rand.seed(3)
results = []
for i in range(5):
    nrep = 4000000
    count = 0
    for i in range(nrep):
        x = sum([rand.random() for j in range(8)])
        if x > 3 and x < 6:
            count += 1
    results.append(count/nrep)
print(results)


    

    




[0.88131075, 0.88124125, 0.88117425, 0.88130075, 0.88114075]

In [6]:

    

n = 8
X = UniformRV(Rational(0),Rational(1))
Y = ConvolutionIID(X,n)


    

In [7]:

    

CDF(Y,6) - CDF(Y,3)


    

    
Out[7]:





$$\frac{35531}{40320}$$

In [8]:

    

# Example 2


    

In [9]:

    

X = TriangularRV(1,2,4)
Y = TriangularRV(1,2,3)
V = X*Y
V.simplify()


    

    




continuous pdf
for 1 <= x <= 2
---------------------------
  4⋅x   2⋅(x + 1)⋅log(x)   4
- ─── + ──────────────── + ─
   3           3           3
---------------------------
 
 
for 2 <= x <= 3
---------------------------
          ⎛ 7⋅x   14    5⋅x    ⎞    
          ⎜ ─── + ──  - ─── - 4⎟    
10⋅x      ⎜  3    3      3     ⎟    
──── + log⎝2        ⋅x         ⎠ - 8
 3                                  
---------------------------
 
 
for 3 <= x <= 4
---------------------------
         ⎛    7⋅x   2  -2⋅x         ⎞    
         ⎜    ─── + ─  ─────        ⎟    
         ⎜     3    3    3    -x - 2⎟    
         ⎜16⋅2       ⋅3     ⋅x      ⎟    
2⋅x + log⎜──────────────────────────⎟ - 4
         ⎝            9             ⎠    
---------------------------
 
 
for 4 <= x <= 6
---------------------------
           ⎛ -7⋅x   -2⋅x   4⋅x   22⎞     
           ⎜ ─────  ─────  ─── + ──⎟     
           ⎜   3      3     3    3 ⎟     
  8⋅x      ⎜2     ⋅3     ⋅x        ⎟   44
- ─── + log⎜───────────────────────⎟ + ──
   3       ⎝         147456        ⎠   3 
---------------------------
 
 
for 6 <= x <= 8
---------------------------
           ⎛    -4⋅x   x  x   4⎞    
           ⎜    ─────  ─  ─ + ─⎟    
           ⎜      3    3  3   3⎟    
  2⋅x      ⎜81⋅2     ⋅3 ⋅x     ⎟   8
- ─── + log⎜───────────────────⎟ + ─
   3       ⎝        256        ⎠   3
---------------------------
 
 
for 8 <= x <= 12
---------------------------
         ⎛  x        x    ⎞    
         ⎜  ─ + 4  - ─ - 4⎟    
2⋅x      ⎜  3        3    ⎟    
─── + log⎝12     ⋅x       ⎠ - 8
 3                             
---------------------------

In [10]:

    

PlotDist(V)
plt.xlabel('v')
plt.ylabel('f(v)')


    

    
Out[10]:





<matplotlib.text.Text at 0x7f935ffe3410>

In [12]:

    

print V.latex()


    

    




\begin{cases} - \frac{4 x}{3} + \frac{1}{3} \left(2 x + 2\right) \log{\left (x \right )} + 1.33333333333333 & \text{for}\: x \geq 1 \\\frac{10 x}{3} + \log{\left (2^{\frac{7 x}{3} + 4.66666666666667} x^{- \frac{5 x}{3} - 4} \right )} - 8 & \text{for}\: x \geq 2 \\2 x + \log{\left (\frac{16}{9} 2^{\frac{7 x}{3} + 0.666666666666667} \cdot 3^{- \frac{2 x}{3}} x^{- x - 2} \right )} - 4 & \text{for}\: x \geq 3 \\- \frac{8 x}{3} + \log{\left (\frac{1}{147456} 2^{- \frac{7 x}{3}} \cdot 3^{- \frac{2 x}{3}} x^{\frac{4 x}{3} + 7.33333333333333} \right )} + 14.6666666666667 & \text{for}\: x \geq 4 \\- \frac{2 x}{3} + \log{\left (\frac{81}{256} 2^{- \frac{4 x}{3}} \cdot 3^{\frac{x}{3}} x^{\frac{x}{3} + 1.33333333333333} \right )} + 2.66666666666667 & \text{for}\: x \geq 6 \\\frac{2 x}{3} + \log{\left (12^{\frac{x}{3} + 4} x^{- \frac{x}{3} - 4} \right )} - 8 & \text{for}\: x \geq 8 \\0 & \text{otherwise} \end{cases}

In [13]:

    

# Example 3


    

In [14]:

    

X = ExponentialRV(Rational(1,2))
Y = ExponentialRV(Rational(1,3))
Z = X+Y
Z.display()


    

    




continuous pdf
for 0 <= x <= oo
---------------------------
   -x     -x 
   ───    ───
    2      3 
- ℯ    + ℯ   
---------------------------

In [15]:

    

# Example 4


    

In [16]:

    

x = Symbol('x',positive=True)
X = RV([x,2*x],[0,1,oo],['continuous','hf'])
simplify(Mean(X))


    

    
Out[16]:





$$\frac{e^{x} - \frac{1}{2}}{x e^{x}}$$

In [17]:

    

X.display()


    

    




continuous hf
for 0 <= x <= 1
---------------------------
x
---------------------------
 
 
for 1 <= x <= oo
---------------------------
2⋅x
---------------------------

In [18]:

    

# Example 5


    

In [22]:

    

x = Symbol('x')
X = RV([x/21],[1,6],['Discrete','pdf'])
Y = ConvolutionIID(X,5)


    

In [23]:

    

[Mean(Y),Variance(Y),Skewness(Y),Kurtosis(Y)]


    

    
Out[23]:





$$\left [ \frac{65}{3}, \quad \frac{100}{9}, \quad - \frac{13}{50}, \quad \frac{1431}{500}\right ]$$

In [ ]:

    




    

In [24]:

    

# Example 3.1


    

In [72]:

    

theta = Symbol('theta',positive=True)
X = ExponentialRV(theta)
Y = ExponentialRV(Rational(1,100))
data = [97, 51, 11, 4, 14, 18, 142, 68, 77, 80, 1, 16, 106, 206,
82, 54, 31, 216, 46, 111, 39, 63, 18, 191, 18, 163, 24]
#P = Posterior(X,Y,data,theta)


    

In [ ]:

    

P


    

In [26]:

    

CredibleSet(P,0.10)


    

    
Out[26]:





$$\left [ 0.00972185579651469, \quad 0.018189624855718\right ]$$

In [27]:

    

Mean(P)


    

    
Out[27]:





$$\frac{28}{2047}$$

In [28]:

    

Mean(P).evalf()


    

    
Out[28]:





$$0.0136785539814362$$

In [34]:

    

X = BinomialRV(12,theta)
Y = TriangularRV(Rational(0),Rational(1,2),Rational(1))
data = [5,6,3,2,5]
P = Posterior(X,Y,data,theta)


    

In [35]:

    

P.display()


    

    




continuous pdf
for 0 <= x <= 1/2
---------------------------
                                       22        39 
-543985284099935107948831513667174400⋅x  ⋅(x - 1)   
────────────────────────────────────────────────────
                 395821499369962291                 
---------------------------
 
 
for 1/2 <= x <= 1
---------------------------
                                      21        40
543985284099935107948831513667174400⋅x  ⋅(x - 1)  
──────────────────────────────────────────────────
                395821499369962291                
---------------------------

In [36]:

    

# Example 3.2


    

In [59]:

    

rand.seed(3)
treatment = [16,23,38,94,99,141,197]
n = len(treatment)
B = 50
x = []
for i in range(B):
    samp = [rand.sample(treatment,1)[0] for i in range(n)]
    med = np.median(samp)
    x.append(med)
np.std(x)


    

    
Out[59]:





$$42.3912302251$$

In [60]:

    

treatment = [16,23,38,94,99,141,197]
X = BootstrapRV(treatment)
Y = OrderStat(X,7,4)
sqrt(Variance(Y)),sqrt(Variance(Y)).evalf()


    

    
Out[60]:





$$\left ( \frac{2 \sqrt{242712738519382}}{823543}, \quad 37.834674393866\right )$$

In [61]:

    

# Example 3.3


    

In [4]:

    

X = MarkovChain(P = [ [Rational(45,100), Rational(48,100),Rational(7,100)],
                    [Rational(5,100),Rational(70,100),Rational(25,100)],
                    [Rational(1,100),Rational(50,100),Rational(49,100)]],
                init = [Rational(5,100),Rational(80,100),Rational(15,100)],
                states = ['stormy','cloudy','sunny'])


    

In [5]:

    

X.display()


    

    




The transition probability matrix:
       stormy cloudy   sunny
stormy   9/20  12/25   7/100
cloudy   1/20   7/10     1/4
sunny   1/100    1/2  49/100
----------------------------------------
The initial system state:
        Prob
stormy  1/20
cloudy   4/5
sunny   3/20

In [6]:

    

X.display(option='steady state',method ='rational')


    

    




The steady state probabilities are:
             Prob
stormy     35/561
cloudy  1399/2244
sunny     235/748

In [13]:

    

X.probability([(3,'stormy')],given=[(1,'cloudy')],method='rational')


    

    
Out[13]:





$$\frac{3}{50}$$

In [8]:

    

X.probability([(5,'cloudy')],method='rational')


    

    
Out[8]:





$$\frac{974185141}{1562500000}$$

In [1]:

    

# Example 3.4


    

In [44]:

    

Y = UniformRV(0,1/2)


    

In [45]:

    

A = 1 - exp(-2*y)
B = exp(-2*y) - 1/2
C = Rational(1,2) - exp(-2* (1-y))


    

In [46]:

    

ys = solve(B-C,y)[0].evalf()
yss = solve(A-C,y)[3].evalf()


    

In [47]:

    

g = [[A,B,C],[0,ys,yss,0.5]]


    

In [ ]:

    




    

In [29]:

    

# Example 3.5


    

In [73]:

    

x=Symbol('x')
X1 = BinomialRV(23, Rational(21,23))
X1 = Transform(X1,[[x/23],[0,23]])
X2 = BinomialRV(28,Rational(27,28))
X2 = Transform(X2,[[x/28],[0,28]])
X3 = BinomialRV(84,Rational(82,84))
X3 = Transform(X3,[[x/84],[0,84]])


    

In [74]:

    

T = X1*X2*X3


    

In [75]:

    

IDF(T,0.05)


    

    
Out[75]:





$$\frac{6723}{9016}$$

In [47]:

    

# Example 3.6


    

In [76]:

    

theta = Symbol('theta',positive=True)
mu = Symbol('mu',positive=True)
X = ExponentialRV(theta)
Y = ExponentialRV(mu)


    

In [77]:

    

T = Queue(X,Y,4,0,1)


    

In [78]:

    

T.simplify()


    

    




continuous pdf
for 0 <= x <= oo
---------------------------
 -3.0        -40.0 ⎛ 4.0        30.0 ⎛       2.0  3        5.0        2       
μ    ⋅(μ + θ)     ⋅⎝μ   ⋅(μ + θ)    ⋅⎝- 4.0⋅μ   ⋅θ ⋅(μ + θ)    - 2.0⋅θ ⋅(μ + θ
──────────────────────────────────────────────────────────────────────────────
                                                                              

 6.0 ⎛ 1.0      2.0⎞                9.0              10.0⎞      5.0           
)   ⋅⎝μ   ⋅θ + μ   ⎠ - 2.0⋅θ⋅(μ + θ)    + 2.0⋅(μ + θ)    ⎠ + 2⋅μ   ⋅θ⋅x⋅(μ + θ
──────────────────────────────────────────────────────────────────────────────
                                                                              

 21.0 ⎛       1.0  3        14.0    2.0          15.0        2        13.0 ⎛ 2
)    ⋅⎝- 1.0⋅μ   ⋅θ ⋅(μ + θ)     + μ   ⋅θ⋅(μ + θ)     + 1.0⋅θ ⋅(μ + θ)    ⋅⎝μ 
──────────────────────────────────────────────────────────────────────────────
                                                                2             

.0      3.0⎞     ⎛ 1.0        ⎞        16.0          18.0⎞                    
  ⋅θ + μ   ⎠ - θ⋅⎝μ    + 1.0⋅θ⎠⋅(μ + θ)     + (μ + θ)    ⎠ + 0.333333333333333
──────────────────────────────────────────────────────────────────────────────
                                                                              

  7.0  3  3        37.0    2  2        32.0 ⎛ 5.0        3.0 ⎛ 2        2.0   
⋅μ   ⋅θ ⋅x ⋅(μ + θ)     + θ ⋅x ⋅(μ + θ)    ⋅⎝μ   ⋅(μ + θ)   ⋅⎝μ ⋅(μ + θ)    + 
──────────────────────────────────────────────────────────────────────────────
                                                                              

     2.0          ⎞        7.0          4.0⎞⎞  -μ⋅x
1.0⋅μ   ⋅θ⋅(μ + θ)⎠ + 1.0⋅μ   ⋅θ⋅(μ + θ)   ⎠⎠⋅ℯ    
───────────────────────────────────────────────────
                                                   
---------------------------

In [79]:

    

Mean(T)


    

    
Out[79]:





$$\frac{1}{\mu \left(\mu + \theta\right)^{22.0}} \left(1.0 \mu \theta^{3} \left(\mu^{1.0} + \theta\right) \left(\mu + \theta\right)^{17.0} + 2.0 \mu \theta^{3} \left(\mu + \theta\right)^{16.0} \left(\mu^{1.0} \theta + \mu^{2.0}\right) + \mu \theta^{2} \left(\mu + \theta\right)^{18.0} \left(1.0 \mu^{1.0} + 2.0 \theta\right) + \mu \theta^{2} \left(\mu + \theta\right)^{19.0} + 2.0 \theta^{3} \left(\mu + \theta\right)^{19.0} + 1.0 \theta \left(\mu + \theta\right)^{21.0} + 1.0 \left(\mu + \theta\right)^{22.0}\right)$$

In [80]:

    

Variance(T)


    

    
Out[80]:





$$\frac{1}{\mu^{11.0} \left(\mu + \theta\right)^{75.0}} \left(14.0 \mu^{9.0} \theta^{3} \left(\mu + \theta\right)^{72.0} + 4.0 \mu^{9.0} \theta \left(\mu + \theta\right)^{74.0} - \mu^{9.0} \left(\mu + \theta\right)^{31.0} \left(1.0 \mu \theta^{3} \left(\mu^{1.0} + \theta\right) \left(\mu + \theta\right)^{17.0} + 2.0 \mu \theta^{3} \left(\mu + \theta\right)^{16.0} \left(\mu^{1.0} \theta + \mu^{2.0}\right) + \mu \theta^{2} \left(\mu + \theta\right)^{18.0} \left(1.0 \mu^{1.0} + 2.0 \theta\right) + \mu \theta^{2} \left(\mu + \theta\right)^{19.0} + 2.0 \theta^{3} \left(\mu + \theta\right)^{19.0} + 1.0 \theta \left(\mu + \theta\right)^{21.0} + 1.0 \left(\mu + \theta\right)^{22.0}\right)^{2} + 2.0 \mu^{9.0} \left(\mu + \theta\right)^{75.0} + 6.0 \mu^{10.0} \theta^{4} \left(\mu + \theta\right)^{70.0} + 10.0 \mu^{10.0} \theta^{3} \left(\mu + \theta\right)^{71.0} + 6 \mu^{10.0} \theta^{2} \left(\mu + \theta\right)^{72.0} + 6.0 \mu^{11.0} \theta^{3} \left(\mu^{1.0} + \theta\right) \left(\mu + \theta\right)^{69.0} + 8.0 \mu^{11.0} \theta^{3} \left(\mu + \theta\right)^{70.0} + 4.0 \mu^{11.0} \theta^{2} \left(\mu + \theta\right)^{71.0}\right)$$

In [81]:

    

# Example 3.8


    

In [82]:

    

Mean(ChiRV(3))/sqrt(3)


    

    
Out[82]:





$$\frac{2 \sqrt{6}}{3 \sqrt{\pi}}$$

In [83]:

    

Mean(ChiRV(99))/sqrt(99)


    

    
Out[83]:





$$\frac{39614081257132168796771975168 \sqrt{22}}{105095150568296034723763017975 \sqrt{\pi}}$$

In [60]:

    

# Example 3.9


    

In [71]:

    

U1 = UniformRV(0,Rational(1))
U2 = UniformRV(0,Rational(1))
V1 = U1 - U2
V1.simplify()


    

    




continuous pdf
for -1 <= x <= 0
---------------------------
x + 1
---------------------------
 
 
for 0 <= x <= 1
---------------------------
-x + 1
---------------------------

In [72]:

    

g1 = [[x*x,x*x],[-oo,0,oo]]
V2 = Transform(V1,g1)


    

In [73]:

    

V2.display()


    

    




continuous pdf
for 0 <= x <= 1
---------------------------
     1 
-1 + ──
     √x
---------------------------

In [104]:

    

#V3 = V2+V2


    

In [70]:

    

V3.display()


    

    




continuous pdf
for 0 <= x <= 1
---------------------------
√x⋅(-log(x) + 4) - log(x) - 4
─────────────────────────────
              √x             
---------------------------

In [73]:

    

g2 = [[sqrt(x)],[0,2]]
V4 = Transform(V3,g2)


    

In [75]:

    

V4.simplify()


    

    




continuous pdf
for 0 <= x <= 1
---------------------------
-4⋅x⋅log(x) + 8⋅x - 4⋅log(x) - 8
---------------------------

In [78]:

    

# Example 3.13


    

In [84]:

    

X = WeibullRV(2,2)


    

In [60]:

    

alist = X.variate(n=10)
blist = X.variate(n=12)
clist = X.variate(n=17)


    

In [4]:

    

alist = [0.2135,0.3153,0.3841,0.3946,0.4707,0.5107,0.5783,0.5960,0.6404,0.7601]
blist = [0.0482,0.0945,0.1149,0.2441,0.287,0.311,0.3362,0.3924,0.4166,0.4194,0.6217,0.698]
clist = [0.0474,0.2732,0.2828,0.2952,0.337,0.3627,0.3698,0.4348,0.4968,0.5061,0.5093,
         0.5211,0.5266,0.6654,0.6951,0.6975,0.7694]


    

In [5]:

    

alist = [0.2135, 0.3153, 0.3841, 0.3946, 0.4707, 0.5107, 0.5783, 0.5960,
           0.6404,0.7601]
blist = [0.0482, 0.0945, 0.1149, 0.2441, 0.287, 0.311, 0.3362, 0.3924,
           0.4166, 0.4194, 0.6217, 0.698]
clist = [0.0474, 0.2732, 0.2828, 0.2952, 0.337, 0.3627, 0.3698, 0.4348,
           0.4968, 0.5061, 0.5093, 0.5211, 0.5266, 0.6654, 0.6951, 0.6975,
           0.7694]


    

In [6]:

    

A = BootstrapRV(alist)
B = BootstrapRV(blist)
C = BootstrapRV(clist)


    

In [7]:

    

Sstar = Minimum( MaximumIID(A,2),
               Minimum(MaximumIID(B,2),
                      MaximumIID(C,2)))
Sstar1 = Minimum(MaximumIID(A,3),
                 MaximumIID(B,3),
                 MaximumIID(C,3))
S = Minimum(MaximumIID(A,2),
           MaximumIID(B,2),
           MaximumIID(C,2))


    

In [8]:

    

Mean(Sstar),Mean(Sstar1), Mean(S)


    

    
Out[8]:





$$\left ( 0.375442473087274, \quad 0.438389764206852, \quad 0.375442473087274\right )$$

In [ ]:

    

n = 2
target = 1.05
m = Mean(Sstar1)
mux = m
while (mux < target * m):
    n+=1
    System = Minimum(MaximumIID(A,n),
                 MaximumIID(B,3),
                 MaximumIID(C,3))
    mux = Mean(System)
    print n,mux
print n


    

In [ ]: