py6. The function choose and percentile in Python

For Python programming, it is best to use a package like Anaconda with almost all the modules. It is easy to update modules using conda command. Lastly if the module is not in conda, you can use the pip program to install it. Remember to use the pip inside Anaconda.

The function np.random.choice selects random values from a given array, with or without replacement. With replacement, we can generate arrays of length greater than that of the the input array.

The np.unique() function finds the unique entries in an array. The percentile function orders the array, and finds the quantiles.

# ex6.py
from __future__ import division, print_function
from numpy.random import choice
import numpy as np
x = choice(100, size = 300, replace = True)
y = len(x)-len(np.unique(x))
print('The number of repeats is %d' % y)
stat = x.mean(), x.std(ddof = 1)
print('The mean is %.2f and standard deviation is %.2f' % stat)
print('The quantiles are:')
v = np.arange(25,85,10)
p = np.percentile(x, v)
print(np.vstack((v,p)))
## The number of repeats is 206
## The mean is 48.65 and standard deviation is 29.26
## The quantiles are:
## [[ 25.    35.    45.    55.    65.    75.  ]
##  [ 24.    34.    42.55  51.    63.7   74.  ]]

R5. Cumulative Binomial Distribution in R

The pbinom function generates cumulative probabilities based on the Binomial Distribution.

A uniform distribution is tested over 1000 trials to get similar curve.

The number of successful outcomes in each trial is 0 (never) to 100 (always) with 50 maximum probability, since the success probability is 0.5.

# ex5.R
num_trials <- 1000
each_trial <- 100
p <- 0.5
x <- 0:each_trial
sig <- numeric(num_trials)
sig1 <- pbinom(x, size = each_trial, prob = p)
for (t in 1:num_trials) {
  tr <- runif(each_trial)
  sig[t] <- sum(tr<=p)
}
sigc <- numeric(length(x))
for (i in x){
  sigc[i+1] <- sum(sig==i)
}
plot(x, cumsum(sigc)/num_trials, type = 'l', col = 'red',
     xlab ='Success', ylab = 'Cumulative Probability',
     main = 'Cumulative Binomial Distribution')
lines(x, sig1, col = rgb(0, 1, 0, 0.5))

Output:

py5. Cumulative Binomial Distribution in Python

We can use scipy.stats.binom.cdf to find cumulative binomial distribution probabilities.

This is contrasted to a uniform distribution generated from 1000 trials, each of size of 100.

The probability of success is 0.5 so the number of successful outcomes ranges from 0 to 100 with 50 most likely.

# ex5.py
from __future__ import division, print_function
from numpy.random import rand
from scipy.stats import binom
import numpy as np
import matplotlib.pyplot as plt
num_trials = 1000
each_trial = 100
p = 0.5
x = np.arange(each_trial+1)
sig = np.zeros(num_trials)
sig1 = binom.cdf(x, n = each_trial, p = p)
for t in range(num_trials):
    tr = rand(each_trial)
    sig[t] = sum(tr<=p)
sigc = np.zeros(len(x))
for i in x:
  sigc[i] = sum(sig==i)
plt.title('Cumulative Binomial Distribution')
plt.plot(x, sigc.cumsum()/num_trials, 'r', 
         x, sig1, 'g')
plt.xlabel('Successes in %d' % each_trial)
plt.ylabel('Cumulative Probability')
plt.show()

Output:

R4. Binomial Distribution in R

The function rbinom gives the binomial distribution

The test distribution is based on 1000 trials with 100 in each trial and success probability of 0.5.

A distribution based on the uniform distribution is tested with that based on the rbinom function.

# ex4.R
num_trials <- 1000
each_trial <- 100
p <- 0.5
sig <- numeric(num_trials)
sig1 <- rbinom(size = each_trial, prob = p, n = num_trials)
for (t in 1:num_trials) {
  tr <- runif(each_trial)
  sig[t] <- sum(tr<=p)
}
hist(sig, breaks = 20, xlab = 'Number of Successes', 
     ylab = 'Frequency', main = 'Binomial Distribution', col='red')
hist(sig1, breaks = 20, , add = T, col=rgb(0, 1, 0, 0.5))

Output:

py4. Binomial Distribution in Python

We can use numpy.random.binomial to create distributions based on binomial distribution.

The probability of success, here, is p = 0.5. We conduct 1000 trials, each of a size of 100. Thus the number of successes for each of 1000 trials is between 0 and 100, but near 50, because of 0.5 probability.

# ex4.py
from __future__ import division, print_function
from numpy.random import binomial,  rand
import numpy as np
import matplotlib.pyplot as plt
num_trials = 1000
each_trial = 100
p = 0.5
sig = np.zeros(num_trials)
sig1 = binomial(each_trial,p,num_trials)
for t in range(num_trials):
    tr = rand(each_trial)
    sig[t] = sum(tr<=p)
plt.subplot(211)
plt.title('Trials')
plt.ylabel('Frequency')
plt.hist(sig, range=(30,70), bins = 20, color = 'r')
plt.subplot(212)
plt.hist(sig1, range=(30,70), bins = 20, color = 'g')
plt.xlabel('Successes in %d' % each_trial)
plt.ylabel('Frequency')
plt.show()

Output:

R3. Cumulative Normal Distribution in R

A for loop is used here to find the fraction of points upto a particular x.

The pnorm function is used to get the actual cumulative normal distribution.

Since the calculated, from the random sample set, is very close to the actual, the resulting plot is mostly magenta, rather than pure red or blue.

# ex3.R
num <- 10000
sig <- rnorm(10000)
x <- seq(from = min(sig), to = max(sig), length.out = 50)
y <- numeric(length(x))
for (i in 1:50) {
  y[i] <- sum(sig<=x[i])/num
}
plot(x, y, type = 'l', col = 'red',
     main = 'Cumulative Normal Distribution')
abline(x, pnorm(x), col = 'blue')

Output:

 

py3. Cumulative Normal Distribution in Python

A for loop is used here to find the fraction of points upto a particular x.

The cdf function from scipy.stats.norm is used to get the actual cumulative normal distribution.

Since the calculated, from the random sample set, is very close to the actual, the resulting plot is mostly magenta, rather than pure red or blue.

# ex3.py
from __future__ import division, print_function
from numpy.random import normal
from scipy.stats import norm
import numpy as np
import matplotlib.pyplot as plt
num = 10000
sig = normal(size=num)
x = np.linspace(min(sig),max(sig))
y = np.zeros(len(x))
for i,v in enumerate(x):
    y[i] = sum(sig<=v)/num
plt.plot(x,y,'b',x,norm.cdf(x),'r')
plt.title('Cumulative Normal Distribution')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

Output: