Effects of CAS on understanding probabilistic concepts

Mihály Klincsik, University of Pécs

Pollack Mihály Faculty of Engineering

Define random variables and working with them

Create random variables and use of the main characteristic features

>    restart:with(Statistics):

Create random variables with known distribution

A random variable X can be creat ed by the Maple command "RandomVariable " from  Statistics package.    

Choose one of the known distribution from Maple list

>    X:=RandomVariable(Normal(mu,sigma));

X := _R

X is a normally distributed  continuous random variable with mean mu  and standard deviation sigma .

Calculations of main features of the random variables

 PDF procedure gives the formula of the Probability Density Function

>    assume(sigma>0):f:=unapply(PDF(X,u),u);

f := u -> 1/2*2^(1/2)/Pi^(1/2)/sigma*exp(-1/2*(u-mu)^2/sigma^2)

The expected value of X can get from the Mean  or ExpectedValue  procedures and by the definition E(X) = int(x*f(x),x = -infinity .. infinity)  using integration

>    Int(x*f(x),x=-infinity..infinity) = ExpectedValue(X);value(%);

Int(1/2*x*2^(1/2)/Pi^(1/2)/sigma*exp(-1/2*(x-mu)^2/sigma^2),x = -infinity .. infinity) = mu

mu = mu

The standard deviation  of X can get by the StandardDeviation procedure and from the definition: STD(X) = E((X-mu)^2)  deduce to the mean procedure.

>    sqrt(ExpectedValue((X-mu)^2))=StandardDeviation(X);

sigma = sigma

Maple can give the formulas of the

(i)   Cumulative Distribution Function of X with "CDF" command: CDF(x) = int(f(t),t = -infinity .. x)

(ii)  Characteristic function of X with  "CharacteristicFunction" command Phi(t) = E(exp(I*t*X))

>    int(f(t),t=-infinity..x)=CDF(X,u);

>    Mean(exp(I*t*X))=CharacteristicFunction(X,t);

1/2+1/2*erf(1/2*2^(1/2)*(x-mu)/sigma) = 1/2+1/2*erf(1/2*2^(1/2)*(u-mu)/sigma)

exp(1/2*I*t*(t*sigma^2*I+2*mu)) = exp(mu*t*I-1/2*sigma^2*t^2)

The list of all procedures which are available for a distribution  we can get by using the " Distribution " command.  

>    Z:=Distribution(Normal(mu,sigma));

Z := module () export Conditions, ParentName, Parameters, CharacteristicFunction, CGF, Mean, Median, MGF, Mode, PDF, StandardDeviation, Variance, CDFNumeric, QuantileNumeric, RandomSample, RandomSample...
Z := module () export Conditions, ParentName, Parameters, CharacteristicFunction, CGF, Mean, Median, MGF, Mode, PDF, StandardDeviation, Variance, CDFNumeric, QuantileNumeric, RandomSample, RandomSample...
Z := module () export Conditions, ParentName, Parameters, CharacteristicFunction, CGF, Mean, Median, MGF, Mode, PDF, StandardDeviation, Variance, CDFNumeric, QuantileNumeric, RandomSample, RandomSample...
Z := module () export Conditions, ParentName, Parameters, CharacteristicFunction, CGF, Mean, Median, MGF, Mode, PDF, StandardDeviation, Variance, CDFNumeric, QuantileNumeric, RandomSample, RandomSample...
Z := module () export Conditions, ParentName, Parameters, CharacteristicFunction, CGF, Mean, Median, MGF, Mode, PDF, StandardDeviation, Variance, CDFNumeric, QuantileNumeric, RandomSample, RandomSample...
Z := module () export Conditions, ParentName, Parameters, CharacteristicFunction, CGF, Mean, Median, MGF, Mode, PDF, StandardDeviation, Variance, CDFNumeric, QuantileNumeric, RandomSample, RandomSample...

These theoretically konwn features are implemented or realized in the practice.

The "Conditions" procedure gives assumptions about the parameters mu  and sigma .

>    Z:-Conditions;

[mu::real, 0 < sigma]

Features of a discrete type random variable are different from the continuous one.

>    U:=Distribution(Poisson(lambda));

U := module () export Conditions, ParentName, Parameters, CDF, CharacteristicFunction, CGF, Kurtosis, Mean, MGF, ProbabilityFunction, Skewness, Variance, VariationCoefficient, CDFNumeric, QuantileNumer...
U := module () export Conditions, ParentName, Parameters, CDF, CharacteristicFunction, CGF, Kurtosis, Mean, MGF, ProbabilityFunction, Skewness, Variance, VariationCoefficient, CDFNumeric, QuantileNumer...
U := module () export Conditions, ParentName, Parameters, CDF, CharacteristicFunction, CGF, Kurtosis, Mean, MGF, ProbabilityFunction, Skewness, Variance, VariationCoefficient, CDFNumeric, QuantileNumer...
U := module () export Conditions, ParentName, Parameters, CDF, CharacteristicFunction, CGF, Kurtosis, Mean, MGF, ProbabilityFunction, Skewness, Variance, VariationCoefficient, CDFNumeric, QuantileNumer...
U := module () export Conditions, ParentName, Parameters, CDF, CharacteristicFunction, CGF, Kurtosis, Mean, MGF, ProbabilityFunction, Skewness, Variance, VariationCoefficient, CDFNumeric, QuantileNumer...
U := module () export Conditions, ParentName, Parameters, CDF, CharacteristicFunction, CGF, Kurtosis, Mean, MGF, ProbabilityFunction, Skewness, Variance, VariationCoefficient, CDFNumeric, QuantileNumer...

We know from the theory that the parameter lambda  of the Poisson distribution is non negative.

>    U:-Conditions;

[0 <= lambda]

>    assume(lambda>0):
Y:=RandomVariable(Poisson(lambda));

Y := _R0

>    p:=unapply(ProbabilityFunction(Y,k),k):p(k);

PIECEWISE([0, k < 0],[lambda^k*exp(-lambda)/k!, otherwise])

>    F:=unapply(CDF(Y,k),k);

F := k -> -1/(floor(k)+1)!*GAMMA(floor(k)+2)+1+(floor(k)+1)/(floor(k)+1)!*GAMMA(floor(k)+1,lambda)

>    plot(subs(lambda=2,p(floor(k))),k=-1..9,title=`Poisson distribution`);

[Maple Plot]

>    plot(subs(lambda=2,F(k)),k=-1..8,thickness=3,discont=true,title=`Cumulative distribution function` );

[Maple Plot]

Create random variables with new distribution

The following function f is appropriate to a probability density function.

>    f:=t -> piecewise(t < 0, 0, t < 1, 1-exp(-t), exp(1-t)-exp(-t)):f(t);

PIECEWISE([0, t < 0],[1-exp(-t), t < 1],[exp(1-t)-exp(-t), otherwise])

>    plot(f(t),t=0..4);

[Maple Plot]

>    Int(f(t),t=0..infinity)=int(f(t),t=0..infinity);simplify(%);

Int(PIECEWISE([0, t < 0],[1-exp(-t), t < 1],[exp(1-t)-exp(-t), otherwise]),t = 0 .. infinity) = exp(1)*exp(-1)

Int(PIECEWISE([1-exp(-t), t < 1],[exp(1-t)-exp(-t), 1 <= t]),t = 0 .. infinity) = 1

Create with "Distribution" command a new distribution object T for wich the PDF is the function f.

>    T := Distribution(PDF = f );

T := module () export PDF, Type; options Distribution, Continuous; end module

Let X random variable with distribution type T.

>    X:=RandomVariable(T);

X := _R1

Probability density function is the given f

>    PDF(X,u);

PIECEWISE([0, u < 0],[1-exp(-u), u < 1],[exp(1-u)-exp(-u), otherwise])

We can get the mean and standard deviation of X.

>    Mean(X),StandardDeviation(X);

3/2, 1/6*39^(1/2)

>   

Manipulating with continuous random variables

>    restart:with(Statistics):

Linear operations  on a continuous random variable

Consider Z as a standard normal random variable, i.e. the mean 0 and standard deviation 1.

>    Z:=RandomVariable(Normal(0,1));
fZ:=PDF(Z,u);

Z := _R

fZ := 1/2*2^(1/2)/Pi^(1/2)*exp(-1/2*u^2)

Define a new random variable R by the operations 2*Z+3 = R , where Z~N(0,1).

>    R:=2*Z+3;

R := 2*_R+3

Maple can calculate the probability density function for the new random variable R such as the original Z.  

>    fR:=PDF(R,u);

>    plot([fZ,fR],u=-3..9,color=[blue,red],thickness=[3,3]);

fR := 1/4*2^(1/2)/Pi^(1/2)*exp(-1/2*(1/2*u-3/2)^2)

[Maple Plot]

If we compare this density function with the general form of  the density function of the normal distribution then we can recognize 2*Z+3 ~ Normal(3,2) where Z~Normal(0,1) .

>    PDF(Normal(3,2),u);

1/4*2^(1/2)/Pi^(1/2)*exp(-1/8*(u-3)^2)

Maple can identify that the repeated sum is a multiplication

>    Z+Z+Z;

3*_R

Operations on discrete version are not working

Similar simple operations are not working  in discrete  cases.

>    U:=RandomVariable(Poisson(1));

>    ProbabilityFunction(U,k);

U := _R1

PIECEWISE([0, k < 0],[exp(-1)/k!, otherwise])

>    U2:=U+U;

U2 := 2*_R1

>    ProbabilityFunction(U2,k);

FAIL

Nonlinear operations  on continuous random variable

Multiplication or power

>    R2:=Z*Z;

R2 := _R^2

>    PDF(R2,u)=PDF(ChiSquare(1),u);

PIECEWISE([0, u <= 0],[1/2*1/u^(1/2)*2^(1/2)/Pi^(1/2)*exp(-1/2*u), 0 < u]) = PIECEWISE([0, u < 0],[1/2*1/u^(1/2)*2^(1/2)/Pi^(1/2)*exp(-1/2*u), otherwise])

We compared  the probablity density function of  the square of a Z~N(0,1) and a Chi-Square with 1 degree of freedom.

We get the theoretical result that Z^2 = (chi^2)[1]

Operations between two random variables

We can continue these investigations with two independent random variables having N(0,1) distribution.

>    Z1:=RandomVariable(Normal(0,1));
Z2:=RandomVariable(Normal(0,1));

Z1 := _R3

Z2 := _R4

In this case the random variables Z1 and Z2 are independent. Consider the sum of their squares!

>    C:=Z1^2+Z2^2;

>    PDF(C,u) = PDF(RandomVariable(ChiSquare(2)),u);

C := _R3^2+_R4^2

PIECEWISE([0, u <= 0],[1/2*exp(-1/2*u), 0 < u]) = PIECEWISE([0, u < 0],[1/2*exp(-1/2*u), otherwise])

The result is well-know that   Z1^2+Z2^2  ~ (chi^2)[2]   , where Z1, Z2  ~ N(0,1)

Operations between two random variables with different distributions

Consider a random variable X having exponential distribution with parameter 1 and Y having uniform distribution on interval [0,1].

>    X:=RandomVariable(Exponential(1));
Y:=RandomVariable(Uniform(0,1));

X := _R6

Y := _R7

Form the sum of X and Y

>    Z:=X+Y;

Z := _R6+_R7

We can get the probability density function of the sum.

>    PDF(Z,u);

>    plot(%,u=-4..4);

PIECEWISE([0, u <= 0],[1-exp(-u), u <= 1],[exp(-u+1)-exp(-u), 1 < u])

[Maple Plot]

>   

Theoretical results without making a great effort

Sum of independent exponentially distributed random variables

Experimenting and make conjecture

Create two independent random variables both having exponential distribution with parameter lambda

>    restart:with(Statistics):

>    assume(lambda>0):
V1:=RandomVariable(Exponential(lambda));
V2:=RandomVariable(Exponential(lambda));;

>    p(x)=PDF(V1,x);

V1 := _R

V2 := _R0

p(x) = PIECEWISE([0, x < 0],[1/lambda*exp(-x/lambda), otherwise])

Create the sum of the two independent exponentially distributed random variable

>    V:=V1+V2;

V := _R+_R0

We obtaine the probability density function of the sum

>    fV:=PDF(V,x);

fV := PIECEWISE([0, x <= 0],[x/lambda^2*exp(-x/lambda), 0 < x])

Maple can't search from the built-in distribution that this is a known distribution. However we can compare it with gamma distribution with parameters lambda  and 2.

>    PDF(Gamma(lambda,2),x);

PIECEWISE([0, x < 0],[x/lambda^2*exp(-x/lambda), otherwise])

  Conjecture

If X[1], X[2], `...`, X[n]  are independent random variables having Exponential(lambda)  distribution, then the distribution of the sum S = X[1]+X[2]+`...`+X[n]  is Gamma with parameters lambda  and n.    

Sum of normally distributed random variables

Experimenting and make conjecture

Creating two indepenedent random variables V1 and V2 with normal distributions and let the parameters are different.

>    restart:with(Statistics):

>    assume(sigma1>0,sigma2>0):
V1:=RandomVariable(Normal(mu1,sigma1));
V2:=RandomVariable(Normal(mu2,sigma2));

V1 := _R

V2 := _R0

Let us form the sum of V1 and v2.

>    S:=V1+V2;

S := _R+_R0

Maple can give the probability density function of the sum S !

>    fS:=PDF(S,x);

fS := 1/2*exp(-1/2*(mu1^2+x^2-2*x*mu2+mu2^2-2*mu1*x+2*mu1*mu2)/(sigma2^2+sigma1^2))/Pi^(1/2)*2^(1/2)/(sigma2^2+sigma1^2)^(1/2)

Now we can compare this function with probability distribtion function having normal distribution with mean mu1+mu2  and standard deviation sqrt(sigma1^2+sigma2^2) .

>    f:=PDF(Normal(mu1+mu2,sqrt(sigma1^2+sigma2^2)),x);

f := 1/2*2^(1/2)/Pi^(1/2)/(sigma2^2+sigma1^2)^(1/2)*exp(-1/2*(x-mu1-mu2)^2/(sigma2^2+sigma1^2))

>    fS-f=simplify(fS-f);

1/2*exp(-1/2*(mu1^2+x^2-2*x*mu2+mu2^2-2*mu1*x+2*mu1*mu2)/(sigma2^2+sigma1^2))/Pi^(1/2)*2^(1/2)/(sigma2^2+sigma1^2)^(1/2)-1/2*2^(1/2)/Pi^(1/2)/(sigma2^2+sigma1^2)^(1/2)*exp(-1/2*(x-mu1-mu2)^2/(sigma2^2+...

Conjecture

If X[1], X[2], `...`, X[n]  are (not necessarily independent) random variables with Normal  distribution, then the distribution of the sum S = X[1]+X[2]+`...`+X[n]  is also normal.

>   

Sum of two independent Poisson distribution.

A little bit harder to get similar results on discrete type random variables because of Maple can't compute the quantities. Have to give more inventions to get the results.

Define two independent random variables X and Y both having Poisson distribution with parameter lambda1  and lambda2 , respectively.

>    restart:with(Statistics):

>    assume(lambda1>0,lambda2>0):
X:=RandomVariable(Poisson(lambda1));
Y:=RandomVariable(Poisson(lambda2));

X := _R

Y := _R0

>    f1:=unapply(ProbabilityFunction(X,k),k);
f2:=unapply(ProbabilityFunction(Y,k),k);

f1 := k -> piecewise(k < 0,0,lambda1^k*exp(-lambda1)/k!)

f2 := k -> piecewise(k < 0,0,lambda2^k*exp(-lambda2)/k!)

Let S is the sum of  X and Y.

>    S:=X+Y;

S := _R+_R0

Maple can't give the distribution of the random variable S.

>    ProbabilityFunction(S,k);

FAIL

We have to determine the probabilities P(S=n)  for n = 0, 1, 2, ....

 Choose an integer n and since S = X+Y and the values of X, Y may also be 0, 1,2,.. therefore

 S = n = 0+n=1+(n-1)=...=k+(n-k) =...=n+0

From the condition of independence we get P(X=k, Y=n-k) =P(X=k)*P(Y=n-k).

Maple can easily give the results for specific n=0, 1, 2, ...., but we can write the approprate sums!

>    seq(factor(sum(f1(k)*f2(n-k),k=0..n)),n=0..5);

exp(-lambda1)*exp(-lambda2), exp(-lambda1)*exp(-lambda2)*(lambda2+lambda1), 1/2*exp(-lambda1)*exp(-lambda2)*(lambda2+lambda1)^2, 1/6*exp(-lambda1)*exp(-lambda2)*(lambda2+lambda1)^3, 1/24*exp(-lambda1)*...
exp(-lambda1)*exp(-lambda2), exp(-lambda1)*exp(-lambda2)*(lambda2+lambda1), 1/2*exp(-lambda1)*exp(-lambda2)*(lambda2+lambda1)^2, 1/6*exp(-lambda1)*exp(-lambda2)*(lambda2+lambda1)^3, 1/24*exp(-lambda1)*...

We can make a conjection on the basis of this values.

However  for the proof of the desired result for all non-negative integer n we have to make some more assumptions and symbolic simplifications.

>    assume(n>0):assume(k>=0):assume(n-k>=0):

>    simplify(sum(f1(k)*f2(n-k),k=0..n),symbolic);

exp(-lambda1-lambda2)*(lambda2+lambda1)^n/GAMMA(n+1)

>    P('S'=n)=convert(%,factorial);

P(S = n) = exp(-lambda1-lambda2)*(lambda2+lambda1)^n/n!

Theorem

The sum S=X+Y of two independent random variables both having Poisson distribution with parameters lambda1  and lambda2  is also Poisson distributed with paramater ( lambda1+lambda2 ).  

>   

A new result

We can make interesting explorations between the built-in functions easily!

Define two independent random variables having standard normal distributions.

>    restart:with(Statistics):

>    Z1:=RandomVariable(Normal(0,1)):
Z2:=RandomVariable(Normal(0,1)):

Investigate the multiplication of the two standard normal distributions.

>    A:=Z1*Z2;

A := _R*_R0

It is amazing that the probability density function can be given by the second kind BesselK-function.

>    fA:=unapply(PDF(A,u),u);

fA := u -> 1/Pi*BesselK(0,(u^2)^(1/2))

We get a new distribution having the BesselK - function as probabiblity density function.

>    plot(fA(u),u=-3..3);
Int(fA(u),u=-infinity..infinity)=int(fA(u),u=-infinity..infinity);

[Maple Plot]

Int(1/Pi*BesselK(0,(u^2)^(1/2)),u = -infinity .. infinity) = 1

>   

Example: Breaking - load of a steel cable which is twisted from wires   

Problem   [Maple Bitmap]

A thick steel cable are twisted from 100 identical thin steel thread. We want to load the cable with M=1960 kg mass. Assume that the breaking stress of each thin threads having normal distribution with the mean m=20 kg and the standard deviation sigma = 4  kg. The strength of the whole cable is the sum of the component forces.

What is the probability that the cable is not breaking under this loading?

(We assume that neither of the thin threads break at an earlier date as the whole cable is breaking!)

Solution

Denote by X the breaking stress of a thin steel thread measuring it in kg. It is a random value and on the basis of the assumptions  X ~ Normal (20 , 4).

Plot the probability density function of X!

>    restart:with(Statistics):

>    X:=RandomVariable(Normal(20,4));

>    plot(PDF(X,u),u=0..40,title=`Probability density function of Normal(20,4)`);

X := _R

[Maple Plot]

>   

Theoretical background

The cable consists of  n=100 thin steel wires. The strength of the cable is the sum of the strength of the wires:

 S = X[1] + X[2] +... + X[100]  ,

where   X[1], X[2], `...,`, X[n]  are independent and distributed by the same Normal(20;4). We know from the theory of random variables that

(i) S is normally distributed;

(ii) the expected value E(S) = sum(E(X[k]),k = 1 .. 100)  = 100*20 = 2000;

(iii) the variance Var(S)= sum(Var(X[k]),k = 1 .. 100)  = 100* 4^2  = 1600, because of the variables are independent.

>   

Practical considerations

We can create 100 random variables and add them by using the following Maple commands:
(a)
"sum"  Maple command;

(b) " add " Maple command;

(c) " for  ...  do  .." loop statement.

However use carefully with these cases!

The "sum" statement is not appropriate creating the sum because of it creates always the same variable. So the sum is the 100 times of one variable.

>    n:=100;

n := 100

>    U:=sum(RandomVariable(Normal(20,4)),k=1..n);

U := 100*_R0

>    Variance(U);

160000

The result obtained by the " add " command is convenient.

>    Z:=add(RandomVariable(Normal(20,4)),k=1..n);

Z := _R10+_R1+_R2+_R3+_R4+_R5+_R6+_R7+_R8+_R9+_R11+_R12+_R13+_R14+_R15+_R16+_R17+_R18+_R19+_R20+_R21+_R22+_R23+_R24+_R25+_R26+_R27+_R28+_R29+_R30+_R31+_R32+_R33+_R34+_R35+_R36+_R37+_R38+_R39+_R40+_R41+...
Z := _R10+_R1+_R2+_R3+_R4+_R5+_R6+_R7+_R8+_R9+_R11+_R12+_R13+_R14+_R15+_R16+_R17+_R18+_R19+_R20+_R21+_R22+_R23+_R24+_R25+_R26+_R27+_R28+_R29+_R30+_R31+_R32+_R33+_R34+_R35+_R36+_R37+_R38+_R39+_R40+_R41+...
Z := _R10+_R1+_R2+_R3+_R4+_R5+_R6+_R7+_R8+_R9+_R11+_R12+_R13+_R14+_R15+_R16+_R17+_R18+_R19+_R20+_R21+_R22+_R23+_R24+_R25+_R26+_R27+_R28+_R29+_R30+_R31+_R32+_R33+_R34+_R35+_R36+_R37+_R38+_R39+_R40+_R41+...
Z := _R10+_R1+_R2+_R3+_R4+_R5+_R6+_R7+_R8+_R9+_R11+_R12+_R13+_R14+_R15+_R16+_R17+_R18+_R19+_R20+_R21+_R22+_R23+_R24+_R25+_R26+_R27+_R28+_R29+_R30+_R31+_R32+_R33+_R34+_R35+_R36+_R37+_R38+_R39+_R40+_R41+...
Z := _R10+_R1+_R2+_R3+_R4+_R5+_R6+_R7+_R8+_R9+_R11+_R12+_R13+_R14+_R15+_R16+_R17+_R18+_R19+_R20+_R21+_R22+_R23+_R24+_R25+_R26+_R27+_R28+_R29+_R30+_R31+_R32+_R33+_R34+_R35+_R36+_R37+_R38+_R39+_R40+_R41+...
Z := _R10+_R1+_R2+_R3+_R4+_R5+_R6+_R7+_R8+_R9+_R11+_R12+_R13+_R14+_R15+_R16+_R17+_R18+_R19+_R20+_R21+_R22+_R23+_R24+_R25+_R26+_R27+_R28+_R29+_R30+_R31+_R32+_R33+_R34+_R35+_R36+_R37+_R38+_R39+_R40+_R41+...

The sum is also acceptable for us using the "for.. do " loop statement.

>    S:=RandomVariable(Normal(20,4));
for i from 2 to 100 do
S:=S+RandomVariable(Normal(20,4));
end do:S;

S := _R101

_R101+_R102+_R103+_R104+_R105+_R106+_R107+_R108+_R109+_R110+_R111+_R112+_R113+_R114+_R115+_R116+_R117+_R118+_R119+_R120+_R121+_R122+_R123+_R124+_R125+_R126+_R127+_R128+_R129+_R130+_R131+_R132+_R133+_R1...
_R101+_R102+_R103+_R104+_R105+_R106+_R107+_R108+_R109+_R110+_R111+_R112+_R113+_R114+_R115+_R116+_R117+_R118+_R119+_R120+_R121+_R122+_R123+_R124+_R125+_R126+_R127+_R128+_R129+_R130+_R131+_R132+_R133+_R1...
_R101+_R102+_R103+_R104+_R105+_R106+_R107+_R108+_R109+_R110+_R111+_R112+_R113+_R114+_R115+_R116+_R117+_R118+_R119+_R120+_R121+_R122+_R123+_R124+_R125+_R126+_R127+_R128+_R129+_R130+_R131+_R132+_R133+_R1...
_R101+_R102+_R103+_R104+_R105+_R106+_R107+_R108+_R109+_R110+_R111+_R112+_R113+_R114+_R115+_R116+_R117+_R118+_R119+_R120+_R121+_R122+_R123+_R124+_R125+_R126+_R127+_R128+_R129+_R130+_R131+_R132+_R133+_R1...
_R101+_R102+_R103+_R104+_R105+_R106+_R107+_R108+_R109+_R110+_R111+_R112+_R113+_R114+_R115+_R116+_R117+_R118+_R119+_R120+_R121+_R122+_R123+_R124+_R125+_R126+_R127+_R128+_R129+_R130+_R131+_R132+_R133+_R1...
_R101+_R102+_R103+_R104+_R105+_R106+_R107+_R108+_R109+_R110+_R111+_R112+_R113+_R114+_R115+_R116+_R117+_R118+_R119+_R120+_R121+_R122+_R123+_R124+_R125+_R126+_R127+_R128+_R129+_R130+_R131+_R132+_R133+_R1...
_R101+_R102+_R103+_R104+_R105+_R106+_R107+_R108+_R109+_R110+_R111+_R112+_R113+_R114+_R115+_R116+_R117+_R118+_R119+_R120+_R121+_R122+_R123+_R124+_R125+_R126+_R127+_R128+_R129+_R130+_R131+_R132+_R133+_R1...

The expected value and the variance of S is coinciding with the theoretically expected

>    expected_value=Mean(S);

>    variance_of_S=Variance(S);

expected_value = 2000

variance_of_S = 1600

Plot the density function of the sum variable S!

>    fS:=PDF(S,u);

>    density:=plot(fS,u=1850..2150,title=`Density function of the breaking-load of the cable`):density;

fS := 1/80*exp(-1250-1/3200*u^2+5/4*u)/Pi^(1/2)*2^(1/2)

[Maple Plot]

The cable is not breaking under the load M=1960 kg, when the cable breaking - load S is greater than the actual loading mass M, i.e. S > 1960. The probability P(S>1960) can visualize by painting the adequate area under the density function.

>    painted:=plot(PDF(S,u),u=1960..2150,filled=true,color=yellow):
plots[display]([density,painted]);

[Maple Plot]

Calculate the probability from the cumulative density function CDF(u)=P(S<u)!
P('S'>1960)=1-CDF(S,1960);
P('S'>1960)=evalf(rhs(%));

P(1960 < S) = 1/2+1/2*erf(1/2*2^(1/2))

P(1960 < S) = .8413447460

So nearly  84% is the probability that the cable does't break under the load M=1960 kg. Only 16% is the chance that the cable is breaking away.

>