Saturday, 25 January 2014

Networks 1: Scraping + Data visualization + Graph stats

These last weeks I have been reading about networks and optimization algorithms, I think is an interesting field with many applications, so my idea was write a new article (or series of articles) showing roughly how use some interesting python libraries like Networkx, for instance. In this article I also like to make a brief introduction to Scraping (a very good way to get data for your graphs).

And only to whet your appetite:

Figure 1. Undirected graph.


Web Scraping:

The first step to construct graphs like the previous one is get the data. As example, I wanted to visualize the whole foreign exchange market and the connections between each currency (Foreign exchange cross rates), but there are another good examples like user connections in social networks (Twitter, Facebook etc).
For this case study(Forex), there are several available websites where you can find financial data for free, nonetheless this data will be restricted to certain time periods, so do not expect 10 or 30 seconds candles for free.

List of Forex websites:

  1. Google Finance
  2. Yahoo Finance
  3. Dukascopy
  4. Investing

Once we have selected our source, we can think about how extract the info. Some websites like Yahoo Finance allow downloading end-of-day equities or forex data, normally in csv format. But let's image that our intention is to obtain the most possible current data. Here Scraping is the appropriate technique. Web scraping is the process of automatically collecting information from websites, usually simulating human interaction through a full-fledged web browser. Note: web scraping may be against the terms of use of several websites.

For this purpose, I used the next python modules:

  1. urllib2
  2. urlparse
  3. re
  4. BeautifulSoup

The goal then is extract all the foreign exchange cross rates. This data will be collected in three numpy arrays and finally we will create a dataframe with pandas to export a csv file. The tricky part of all this is parse the html file, however BeautifulSoup allows using regular expressions, which combined with urlparse creates a powerful and simple tool. So let's take a look at the script.

Script:

import numpy as np
import pandas as pd
import urllib2
import urlparse
import re
from bs4 import BeautifulSoup

#------------------------------------------------------------------------------
#Functions:
#------------------------------------------------------------------------------


def add_exchange(Data):
    Value = []
    for link in Data[0, :]:
        soup = connect(link)
        for span in soup.find_all('span', {'id': 'last_last'}):
            #remove thousands (",")
            value = (span.text).replace(',', '')
            Value.append(value)
    return np.vstack((Data, Value))


def connect(site):
    hdr = {'User-Agent': 'Mozilla/5.0 (X11; Linux x86_64)',
       'Accept': 'text/html,application/xhtml+xml,application/xml;q=0.9,*/*;q=0.8',
       'Accept-Charset': 'ISO-8859-1,utf-8;q=0.7,*;q=0.3',
       'Accept-Encoding': 'none',
       'Accept-Language': 'en-US,en;q=0.8',
       'Connection': 'keep-alive'}
    req = urllib2.Request(site, headers=hdr)
    page = urllib2.urlopen(req).read()
    soup = BeautifulSoup(page)
    return soup

#------------------------------------------------------------------------------
#Extract list of currencies + Matrix
#------------------------------------------------------------------------------
mainSites = {'http://www.investing.com/currencies/middle-east',
    'http://www.investing.com/currencies/north-america',
    'http://www.investing.com/currencies/central-america',
    'http://www.investing.com/currencies/south-america',
    'http://www.investing.com/currencies/asia',
    'http://www.investing.com/currencies/pacific',
    'http://www.investing.com/currencies/europe',
    'http://www.investing.com/currencies/middle-east',
    'http://www.investing.com/currencies/africa',
    'http://www.investing.com/currencies/caribbean'}

FxUrls = []
C1 = []
C2 = []
for mainSite in mainSites:
    soup = connect(mainSite)
    for link in soup.find_all(href=re.compile("www.investing.com/currencies/")):
        #extract currencies:
        url = link.get('href')
        FxUrls.append(url)
        currencies = urlparse.urlparse(url).path.split("/")[2].split("-")
        C1.append(currencies[0])
        C2.append(currencies[1])

#Create Matrix
Forex = np.vstack((FxUrls, C1, C2))
#------------------------------------------------------------------------------

#Extract all foreign exchange cross rates:
Import = add_exchange(Forex)

#Pandas dataframe
data = {'currency_1': Import[1, :],
    'currency_2': Import[2, :],
    'value': Import[3, :]}

df = pd.DataFrame(data)
df.to_csv('Forex.csv', index=False, cols=['currency_1', 'currency_2', 'value'])

Summarizing what we have done in this script:

I wrote two functions, one in order to connect with website and capture the html page, and the other one to parse this html and select the required information. There are different ways to extract the data, select the one that best suits the structure of your html file. Run the script, just take into account that this script may take a while depending on your internet connection.
With this tool, you might consider a dynamic system to obtain frequently data. Well, that would be complicated and unlikely, your IP address would be blocked if the number of requests is considerable.


Networkx - Graphs:

Now we have collected the data, it's time to show what we got. For the creation and study of the graphs we will use the python package called Networkx. This is a pretty helpful package which provides many functions and algorithm that facilitate the study of complex networks. However, this may be not the best package to draw your graph, although we can use it as an interface to use Graphviz.

The first thing I did was write a new module, called "graphs", with the Networkx's functionalities I needed. This includes two classes, "create_graph"and "graphStats". The first one reads the produced csv file and creates the two elements that form a graph, nodes and edges. In this case we also have weights, the exchange rates (third column), so we need a list in order to create the labels, which looks like this:

example: {('D', 'E'):9.0, ('F', 'N'):6.0}...

The class "graphStats" includes two functions in order to calculate the degree centrality and betweenness centrality, we will talk about this in more detail later on.

Graph class:
import networkx as nx
from operator import itemgetters
import pandas as pd
import numpy as np


class create_graph:

    def __init__(self):
        data = 'csv/Forex.csv'
        headers = ['start', 'end', 'weights']
        self.Network = pd.read_csv(data, skiprows=1, names=headers)
        self.graph = np.asarray(self.Network[['start', 'end']])
        self.nodes = np.unique(self.graph)
        self.weights = list(map(float, self.Network['weights']))

    def networkList(self):
        G = nx.DiGraph()
        for node in self.nodes:
            G.add_node(node)

        self.graph_l = []
        for edge in self.graph:
            G.add_edge(edge[0], edge[1])
            self.graph_l.append((edge[0], edge[1]))

        labels = dict(list(zip(self.graph_l, self.weights)))
        return G, labels


class graphStats:

    def calculate_degree_centrality(self, graph):
        dgc_key = []
        dgc_value = []
        g = graph
        dc = nx.degree_centrality(g)
        nx.set_node_attributes(g, 'degree_cent', dc)
        degcent_sorted = sorted(dc.items(), key=itemgetter(1), reverse=True)
        for key, value in degcent_sorted:
            dgc_key.append(str(key))
            dgc_value.append(value)
        return dgc_key, dgc_value

    def calculate_betweenness_centrality(self, graph):
        btc_key = []
        btc_value = []
        g = graph
        bc = nx.betweenness_centrality(g)
        betcent_sorted = sorted(bc.items(), key=itemgetter(1), reverse=True)
        for key, value in betcent_sorted:
            btc_key.append(str(key))
            btc_value.append(value)
        return btc_key, btc_value


Now we have prepared the module, we will use matplotlib to draw the graph. Actually, this is not a huge network, it has only 133 nodes and 1918 edges, but it is difficult get a good resolution if your output file is .jpg or png, so I recommend a pdf file and some adjustments in the figure size.

Graph - plot:
import networkx as nx
import matplotlib.pyplot as plt
import graphs

#------------------------------------------------------------------------------
# Read data + create graph
#------------------------------------------------------------------------------
gr = graphs.create_graph()
graph = gr.graph
G, labels = gr.networkList()

#------------------------------------------------------------------------------
# Graph plot:
#------------------------------------------------------------------------------
fig = plt.figure()
#fig = plt.figure(figsize=(8,8), dpi=1000) - pdf file
pos = nx.graphviz_layout(G, prog='neato')
node_color = [float(G.degree(v)) for v in G]
nx.draw_networkx_nodes(G, pos, node_size=[float(G.degree(v)) * 2.5 for v in G],
             alpha=0.6, node_color=node_color)
nx.draw_networkx_edges(G, pos, alpha=0.25, edge_color='blue', width=0.5,
arrows=False)
nx.draw_networkx_labels(G, pos, font_size=2.5)

plt.axis('off')
plt.show()
#plt.savefig("directedGraph.pdf")

Figure 1.1. Undirected graph (detail).
Figure 2. Directed graph.

The second graph is which shows the real Forex network. The Forex network is a directed graph, a graph where its edges have a direction associated with them. This digraph is strongly connected due to it contains several pairs of vertices with both paths, (euro -> dolar and dolar -> euro, for instance).

Graph - stats:

Networkx is a great python package to study graphs, its main purpose. In this point we are going to look into some statistics which will help us to extract real information from our graphs with the aim of interpret the results properly.

1) Betweenness centrality:

Betweenness centrality quantifies the number of times a node acts as bridge along the shortest path between two other nodes. It is basically a measure of connectedness between components of the graph.

2) Degree centrality:

The centrality for a node measures its relative importance within the graph. The degree centrality for a node (n) is the number of nodes it is connected. The degree centrality values are normalized by the maximum degree in the graph. For a graph with self loops (a digraph) the maximum degree might be higher than 1 and the degree centrality as well.

In this article I will not expand the theory behind these measures.
import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np
import graphs


def topTable(field1, field2, n_top):
    topM = max(field2) * 0.9
    right = len(field1) * 0.75
    plt.text(right, topM * 1.08, 'Top %s' % n_top, fontsize=12)
    for i in range(n_top):
        curr = field1[i]
        val = field2[i]
        plt.text(right, topM - i * topM / 20, '{}) {} = {}'.format(i + 1,
        curr.upper(), round(val, 3)), fontsize=8)

#------------------------------------------------------------------------------
# Read data + create graph
#------------------------------------------------------------------------------
gr = graphs.create_graph()
G, labels = gr.networkList()

#------------------------------------------------------------------------------
# Stats + plots:
#------------------------------------------------------------------------------
stats = graphs.graphStats()
betc_key, betc_value = stats.calculate_betweenness_centrality(G)
degc_key, degc_value = stats.calculate_degree_centrality(G)

#Average degree
N = G.order()
K = G.size()
avg_d = float(N) / K
avg_degree = 'Average degree: %.4f ' % (avg_d)


# Plot: Degree_centrality
plt.figure()

ax1 = plt.subplot(211)
plt.title('Degree centrality for nodes', fontsize=12)
a_lenght = np.arange(len(degc_value))
plt.bar(a_lenght, degc_value, color=cm.jet(degc_value), align='center')
plt.xticks(a_lenght, degc_key, size='small', rotation='vertical')
plt.tick_params(axis='x', labelsize=4)
plt.tick_params(axis='y', labelsize=8)
plt.autoscale(enable=True, axis='both', tight=None)

#Top degree centrality:
topTable(degc_key, degc_value, 10)
plt.text(len(degc_value) * 0.75, max(degc_value) * 0.4, avg_degree,
bbox={'facecolor': 'blue', 'alpha': 0.25, 'pad': 10}, fontsize=7)

# Plot: Betweenness_centrality
plt.subplot(212)
plt.title('Betweenness centrality for nodes', fontsize=12)
a_lenght = np.arange(len(betc_value))
plt.bar(a_lenght, betc_value, color=cm.jet(betc_value), align='center')
plt.xticks(a_lenght, betc_key, size='small', rotation='vertical')
plt.tick_params(axis='x', labelsize=4)
plt.tick_params(axis='y', labelsize=8)
plt.autoscale(enable=True, axis='both', tight=None)
plt.ylim(0, max(betc_value) * 1.1)
plt.plot(betc_value, '--b')

#Top degree betweenness:
topTable(betc_key, betc_value, 10)

plt.show()


The above plot shows the currencies which have more relevance into the network, and as we suspected these are: USD dolar, Euro and the pound sterling. We can make at least two significant considerations; The average degree centrality is quite low, what shows that there are currencies with a notable grade of popularity(Degree centrality : in-degree) and these same currencies play a broker role in the network due to cumulate most of trajectories for the shortest paths. For this particular example,this might mean that we can create a more simplistic network, which would represent pretty faithfully the Forex market, excluding those currencies with fewer number of connections.


Adjacency Matrix:

The adjacency matrix is a N x N matrix where the cells constitute edges. In this example the fields are coloured according the degree centrality and as the graph is directed, the matrix is not symmetric.
import numpy as np
import matplotlib.pyplot as plt
import graphs


gr = graphs.create_graph()
G = gr.networkList()[0]

n = len(gr.nodes)
oMapp = dict(list(zip(gr.nodes, np.arange(n))))

Mtx_w = (np.ones((n, n))) * np.inf
for i in range(len(gr.graph)):
    start = oMapp.get(gr.Network['start'][i])
    end = oMapp.get(gr.Network['end'][i])
    w = gr.Network['weights'][i]
    Mtx_w[start, end] = w

Mtx_nodes = []
Mtx_nodes = np.array([np.concatenate([Mtx_nodes, gr.nodes]) for _ in range(n)])
#-----------------------------------------------------------------------------

# Show matrix by colour:
Mtx_Image = Mtx_w.copy()
for i in range(n):
    count = ((Mtx_Image[i, :] != np.inf) & (Mtx_Image[i, :] != np.nan)).sum()
    Mtx_Image[i, :][Mtx_Image[i, :] != np.inf] = count

# Gradient based on degrees
plt.figure()
plt.title('Foreign Exchange cross rates', fontsize=16)
plt.imshow(Mtx_Image, interpolation='nearest')
plt.xticks(np.arange(len(gr.nodes)), gr.nodes, rotation=90)
plt.yticks(np.arange(len(gr.nodes)), gr.nodes)
plt.tick_params(axis='both', labelsize=4)

plt.tight_layout()
plt.show()
#plt.savefig("Forex_Matrix.pdf")




This first article related to networks has been just a brief introduction on this field. We will continue exploring graphs and graph algorithms in next posts. Finally, I would like to make some notes regarding the network images, these images have been modified reversing the RGB output (I preferred a black background), I used GIMP for that task.

Sunday, 5 January 2014

Numerical methods: ODE and Finite difference method

The goal of this post is show how solve an ordinary differential equation, numerically and using the finite difference method and compare the result with the analytic solution. The problem statement is as follow:

Numerically solve the ODE: $m\frac{\partial^{2}y}{\partial t^{2}} = -b\frac{\partial y}{\partial t} $ that represents an object whose velocity at $t=0$ is equal to $v_0$ and it is slowing due to a retardant force, which is proportional to its velocity.

The integration interval is $0 \leq t \leq 5$ seconds. For simplicity, we take $m = b = 1$.
Boundary conditions:

  • $y(0) = 0 m$
  • $y(5) = 1 m$

Analytic solution:
$Y(t) = \frac{m v_0}{b} (1 - e^{-\frac{b}{m}t})$ where $v_0 = 1$

Ordinary differential equation:
We will start solving the ordinary differential equation of second grade . I will explain the resolution step by step:

  • Step 1: Characteristic equation
  • Step 2: Finding constants of integration
$m\frac{\partial^{2}y}{\partial t^{2}} = -b\frac{\partial y}{\partial t} \Rightarrow mD^{2} = -bD \Rightarrow D(mD + b) = 0 $
Thus, we have two possible solutions: $s1 = 0; s2 = -\frac{b}{m}$. As these two roots are real and distinct ($s1 \ne s2$) we can use the general solution for this specific case:

$Y(x) = C1 e^{s_1 x} + C2 e^{s_2 x} \Rightarrow Y(t) = C1 e^{0 t} + C2 e^{-\frac{b}{m} t} = C1+C2 e^{-\frac{b}{m} t}$

Using the boundary conditions we can obtain the constants of integration C1 and C2:
$Y(t=0) \Rightarrow C1 = -C2$
$Y(t=5) \Rightarrow -C2 + C2e^{-\frac{b}{m} 5 } = 1 \Rightarrow (-1 + e^{-\frac{b}{m} 5 } ) C2 = 1 \Rightarrow C2 = \frac{1}{(-1 + e^{-\frac{b}{m} 5 } ) } $
Solution:
$Y(t) = -\frac{1}{(-1 + e^{-\frac{b}{m} 5 } ) } + \frac{1}{(-1 + e^{-\frac{b}{m} 5 } ) } e^{-\frac{b}{m} t}$

$m = b = 1 \Rightarrow Y(t) = 1.00678 - 1.00678 e^{-t}$

Finite difference method:
We use the centered difference at time t to obtain the central difference approximation in both equation sides.
\[y" = \frac{\partial^{2}y}{\partial t^{2}} = \frac{U_{i+1} - 2U_i + U_ {i-1}}{h^{2}} \enspace y' = \frac{\partial y}{\partial t} = \frac{U_{i+1}-U_{i-1}}{2h} \]
Now we can sustitude these above approximations in the initial expression:
\[\frac{m}{h^{2}}(U_{i+1} - 2U_i + U_ {i-1}) = -\frac{b}{2h}(U_{i+1}-U_{i-1})\] Isolate $U_i$:
$U_{i+1} - 2U_i + U_ {i-1} = -\frac{b}{2h}(U_{i+1}-U_{i-1})\frac{h^{2}}{m} \Rightarrow [-\frac{bh}{2m}(U_{i+1}-U_{i-1})-(U_{i+1}+U_{i-1})] \frac{1}{-2} = U_i$

Implementation:
For this simple example, we create a Python class called "ODE", which includes the analytic solution and the basic methodology based on the iterative scheme.

Iterations:

def solve(self, n, steps, m, b, h):
    v = np.zeros(n)
    #Boundary conditions
    #y(t = 0) = 0
    v[0] = 0
    #y(t = 5) = 1
    v[n - 1] = 1

    #Iterations
    for i in range(steps):
        for i in range(1, n - 1):
            a1 = (v[i - 1] - v[i + 1]) * (b * h / (2 * m))
            a2 = v[i + 1] + v[i - 1]
            v[i] = (a1 - a2) / (-2)
    return v
As shown, the first value in the array $= 0$, first boundary condition (we start the loop in the position 1) and the last value = 1 (position n - 1). Finally, the function returns the complete array.

The attached chart shows how the result is approximated by finite difference as the number of iterations increases.



Code:
import numpy as np
import matplotlib.pyplot as plt


class ODE(object):

    def __init__(self, m, b, v0, t_min, t_max, n):
        self.t = np.linspace(t_min, t_max, n)
        #Analytic solution
        self.s = m * v0 / b * (1 - np.exp(- b / m * self.t))
        #ODE
        C2 = 1 / (- 1 + np.exp(- b / m * t_max))
        C1 = -C2
        self.f = C1 + C2 * np.exp(- b / m * self.t)

    def solve(self, n, steps, m, b, h):
        v = np.zeros(n)
        #Boundary conditions
        #y(t = 0) = 0
        v[0] = 0
        #y(t = 5) = 1
        v[n - 1] = 1

        #Iterations
        for i in range(steps):
            for i in range(1, n - 1):
                a1 = (v[i - 1] - v[i + 1]) * (b * h / (2 * m))
                a2 = v[i + 1] + v[i - 1]
                v[i] = (a1 - a2) / (-2)
        return v

#inputs and solution:
m = 1
b = 1
v0 = b / m
t_min = 0.0
t_max = 5.0
n = 100
h = t_max / n

vel = ODE(m, b, v0, t_min, t_max, n)
plt.plot(vel.t, vel.f, 'b.-', label='ODE')
plt.plot(vel.t, vel.s, 'g.-', label='Analytic solution')

for i in range(250, 3000, 500):
    sol_i = vel.solve(n, i, m, b, h)
    plt.plot(vel.t, sol_i, '-.', label=r'Iteration = %i' % i)

plt.ylabel('Distance (m)')
plt.xlabel('Time (s)')

plt.legend(loc=2, prop={'size': 9})
plt.show()

Saturday, 30 November 2013

Weibull analysis

This first post in ComputSimu will explain the basic methodology to calculate Weibull distribution given a data set, in this case a data set provided by a three point bending test for specimens with circular cross section. The Weibull distribution is one of the most used distributions to model life data, due to its extreme flexibility to fit a wide range of data (Normal distribution or logarithmic distribution, for instance) and its applicability modelling different problems (weather forecasting, failure analysis, delivery times, etc...)

For creating the plot and calculate the Weibull distribution parameters we will use a python script (www.python.org). Python is widely extended as scripting language in computational science with an impressive support, which guarantees a constant number of new modules and improvements.
Schema:
  1. Three point bending test
  2. Weibull distribution
  3. Characteristic life
  4. Weibull distribution characteristics
  5. Extra: Gamma function
  6. Code

1) Three point bending test
This test consists in an applied force (F) acting at the center of a beam to assess the fracture resistance. Required formulas:

Beam deflection (double pinned contact):

$Moment(max) =\dfrac{F}{2}\dfrac{L}{2} = \dfrac{FL}{4}$ where \(L = span\) (mm) and \(F = load\) (N)
$Stress(max) = \dfrac{M_t}{\dfrac{I_p}{r}} = \dfrac{M_t}{W_p}$ where \(r = radius\) (mm); \(W_p = W_y = W_z = \dfrac{I_y}{r} = \dfrac{\dfrac{\pi r^4}{4}}{r} = \dfrac{\pi r^3}{4}\)
$Stress(max) = \dfrac{Moment(max)}{Resisting Bending Moment} = \dfrac{\dfrac{FL}{4}}{\dfrac{\pi r^3}{4}} = \dfrac{FL}{\pi r^3}$ (MPa)
We will use Numpy to import the data (there are many options, but this works pretty well). The completed code is on the bottom of the post, but I will explain quickly the most important parts. The first step is import the data into an array. These values are ordered ascendant, which is a must in order to calculate the failure probability, but you can also order your array once the values are imported. Example:

F(N)
333
440
440
...

 path = '../Data_Test.csv'
 aSamples = np.genfromtxt(path, delimiter=',', skip_header=1, dtype='f')
 
Now, we need two functions to create two new arrays with failure probabilities and the failure stresses. To calculate the probability we can choose one of the two methods below depending on the number of samples.

def FailureProb(n, N):
    if(N < 100):
        #Bernard's approximation
        return (n - 0.3) / (N + 0.4)
    else:
        #Mean Ranks
        return n / (N + 1)


def MaxStress(F, L, r):
    # Maximum Stress (cylinder 3 points)
    return (F * L) / (math.pi * r ** 3)

aSamples = np.genfromtxt(path, delimiter=',', skip_header=1, dtype='f')

aFailureProb = aSamples.copy()
aFailureStress = aSamples.copy()

for i in range(len(aFailureProb)):
    aFailureProb[i] = FailureProb((i + 1), float(len(aSamples) + 1))
    aFailureStress[i] = MaxStress(aSamples[i], L, r) 

2) Weibull distribution
The probability density function of a random variable is: \[f(x; \alpha, \beta) = \begin{cases} \frac{\alpha}{\beta} (\frac{x}{\alpha})^{(\alpha - 1)} e^{-(\frac{x}{\beta})^{\alpha}} & \text{if } x \geq 0 \\ 0 & \text{if } x < 0 \end{cases}\]
and the cumulative distribution function is: \[F(x; \alpha, \beta) = \begin{cases} 1-e^{-(\frac{x}{\beta})^{\alpha}} & \text{if } x \geq 0 \\ 0 & \text{if } x < 0 \end{cases}\]
Shape parameter ($\alpha$):
$\alpha$ shows how the failure rate develops in time. We can describe three different situations for the parameter $\alpha$:

  • $\alpha$ < 1: This first situation is characterized by a decreasing failure rate or hazard function.
  • $\alpha$ = 1: The second situation is characterized by a constant failure rate (independent of time), this is called lack of memory. In this situation Weibull distribution is identical to the exponential distribution.
  • $\alpha$ > 1: The third situation is characterized by an increasing failure rate. This represents that the more older is the "element" the more likely it is to fail.
Scale parameter or characteristic life ($\beta$):
$\beta$ is the time/cycle/sample at which 63,2% are expected to fail. This is a constant for all Weibull distribution, regardless the shape parameter.

Plot:

Script:
"""
Weibull:
    1) Probability density function
    2) Cumulative distribution function
"""

import numpy as np
import matplotlib.pyplot as plt
from weibull import weibProbDist, weibCumDist

x = np.linspace(0, 5, 100)
a_alpha = np.array([0.5, 1, 2, 5])
label = r'${\alpha} = %s, {\beta} = %s$'

fig = plt.figure('Weibull')
ax1 = fig.add_subplot(211)
ax2 = fig.add_subplot(212)
for i in a_alpha:
    ax1.plot(x, weibProbDist(x, i, 1), label=label % (i, 1))
    ax2.plot(x, weibCumDist(x, i, 1), label=label % (i, 1))
ax1.set_title('Probability density function', size=12)
ax2.set_title('Cumulative distribution function', size=12)
ax1.legend(loc=1, prop={'size': 10})
ax2.legend(loc=7, prop={'size': 10})

fig.subplots_adjust(hspace=0.5)
plt.show()

3) Characteristic life
To obtain the Weibull parameters the first step is linearise the cumulative distribution functions: $F(x) = 1-e^{-(\frac{x}{\beta})^{\alpha}} \Rightarrow P_f = 1-e^{-(\frac{\sigma}{\sigma_0})^{m}}$ where m = slope and \(\sigma\) = scale (characteristic stress)

$P_s = 1-P_f$
$ln(P_s^{-1}) = ln(e^{(\frac{\sigma}{\sigma_0})^{m}}) \Rightarrow ln[ln(P_s^{-1})] = m \enspace ln(\frac{\sigma}{\sigma_0}) \Rightarrow ln[ln(P_s^{-1})] = m \enspace ln(\sigma) - m \enspace ln(\sigma_0)$

Linear expression:
$\underbrace{ln[ln(P_s^{-1})]}_\text{f(x)} = \underbrace{ m \enspace ln(\sigma)}_\text{mx} - \underbrace{ m \enspace ln(\sigma_0)}_\text{b}$
Maximum likelihood estimation (MLE) - Linear regression:
$m (slope)= \dfrac{Cov(ln(\sigma), ln[ln(P_s^{-1})])}{Var(ln(\sigma))} = \dfrac{\dfrac{1}{n} \sum_{i=1}^{n}(ln(\sigma) - \overline{ln(\sigma)}) (ln[ln(P_s^{-1})] - \overline{ln[ln(P_s^{-1})])}}{\dfrac{1}{n} \sum_{i=1}^{n}(ln(\sigma) - \overline{ln(\sigma)})^2}$
$b (intercept)= \overline{ln[ln(P_s^{-1})])} - m \enspace \overline{ln(\sigma)}$
In order to solve this simple equation we use Scipy and its module called stats which apply the described methodology above.
 #Survival probability = 1 - Failure Probability
 aSurvivalProb = 1 - aFailureProb

 #Linear regression: parameters (m, b)
 aLnPs = np.log(- np.log(aSurvivalProb))
 aLnFS = np.log(aFailureStress / 1000)

 slope, intercept, r_square, s, std_error = stats.linregress(aLnFS, aLnPs)
 regressLine = slope * aLnPs + intercept

 #Characteristic life : Stress
 Stress = np.exp(- (intercept / slope)) * 1000
 Pr_Stress = weibCumDist(Stress, slope, Stress)
 
Finally, we can calculate the parameter $\beta$ with the next expression:
$\beta, \sigma = e^{-(\frac{intercept(b)}{slope(m)})}$

Plotting analysis results:
As described above, the final plot shows how a stress = 125,08 MPa is the characteristic stress whereby the failure probability is equal to 63,2%. The plots below try to be self-explanatory.




4) Weibull distribution characteristics
You will find all characteristic parameters in the Weibull module. These parameters are exported into a .csv file.
$Mean = \beta \enspace \Gamma(1+\dfrac{1}{\alpha})$
$Median = \beta \enspace \ln(2)^{(\frac{1}{\alpha})}$
$Mode = \begin{cases} \beta \enspace (\frac{\alpha - 1}{\alpha})^{\frac{1}{\alpha}} & \text{if } \alpha > 1 \\ 0 & \text{if } \alpha = 1 \end{cases}$
$Variance = \beta^2 \enspace \Gamma(1+\dfrac{2}{\alpha}) - \mu^2$
$Skewness = \beta^2[\Gamma(1+\dfrac{3}{\alpha})+2(\Gamma(1+\dfrac{1}{\alpha}))^3 - 3[\Gamma(1+\dfrac{2}{\beta})][\Gamma(1+\dfrac{1}{\beta})]]$
Skewness
  • $\text{if } \alpha < 3.6: Skewness > 0 $
  • $\text{if }\alpha = 3.6: Skewness = 0 $
  • $\text{if }\alpha > 3.6: Skewness < 0 $

5) Extra: Gamma function

The Weibull module also include one of the possible approximations for the Gamma function, Lanczos approximation. You can find more information about different implementations on the links below:
C++ library + pseudocode
"Calculators and the Gamma Function", Viktor T. Toth (Recommended)

6) Code
import numpy as np
import matplotlib.pyplot as plt
import math
from scipy import stats
from weibull import *
import csv

#inputs:
L, r = float(40), float(4)
path = '../Data_Test.csv'
pathReport = '../WeibullReport.csv'


def FailureProb(n, N):
    if(N < 100):
        #Bernard's approximation
        return (n - 0.3) / (N + 0.4)
    else:
        #Mean Ranks
        return n / (N + 1)


def MaxStress(F, L, r):
    # Maximum Stress (cylinder 3 points)
    return (F * L) / (math.pi * r ** 3)

aSamples = np.genfromtxt(path, delimiter=',', skip_header=1, dtype='f')

aFailureProb = aSamples.copy()
aFailureStress = aSamples.copy()

for i in range(len(aFailureProb)):
    aFailureProb[i] = FailureProb((i + 1), float(len(aSamples) + 1))
    aFailureStress[i] = MaxStress(aSamples[i], L, r)

MeanForces = np.mean(aSamples)
MeanStress = np.mean(aFailureStress)

#Survival probability = 1 - Failure Probability
aSurvivalProb = 1 - aFailureProb

#Linear regression: parameters (m, b)
aLnPs = np.log(- np.log(aSurvivalProb))
aLnFS = np.log(aFailureStress / 1000)

slope, intercept, r_square, s, std_error = stats.linregress(aLnFS, aLnPs)
regressLine = slope * aLnPs + intercept

m = slope
b = intercept

#Characteristic life : Stress
Stress = np.exp(- (intercept / slope)) * 1000
Pr_Stress = weibCumDist(Stress, slope, Stress)

# export results to .csv
Results = [
            ['Parameters', 'Value', 'Units'],
            ['shape (alpha)', slope],
            ['scale (beta)', Stress, 'MPa'],
            ['Mean', weibMean(slope, Stress), 'MPa'],
            ['Median', weibMean(slope, Stress), 'MPa'],
            ['Mode', weibMode(slope, Stress), 'MPa'],
            ['Variance', weibVar(slope, Stress), 'MPa^2'],
            ['Standard Deviation', np.sqrt(weibVar(slope, Stress)), 'MPa'],
            ['Skewness', weibSkew(slope, Stress)]
          ]

with open(pathReport, 'wb') as f:
    writer = csv.writer(f)
    writer.writerows(Results)

#Plots:
fig = plt.figure('Weibull')
ax1 = fig.add_subplot(211)
scatt = ax1.scatter(aLnFS, aLnPs, label='samples', c='DarkTurquoise', s=15)
regres = ax1.plot(aLnPs, regressLine, 'b-', label='Linear regression (MLE)')
ax1.set_ylabel('$Ln[-Ln(Ps)]$')
ax1.set_xlabel('$Ln(\sigma) : Ln(GPa)$')
ax1.set_title('Weibull analysis')
ax1.set_ylim(min(aLnPs) * 1.1, max(aLnPs) * 1.25)
ax1.set_xlim(min(aLnFS) * 1.1, max(aLnFS) * 0.9)
ax1.legend(loc=2, prop={'size': 10})
ax1.text(min(aLnFS) * 1.05, -1, r'$y=%.4f m + %.4f$' % (m, b))
ax1.text(min(aLnFS) * 1.05, -1.5, r'$R^2=%.4f$' % r_square)

x = np.linspace(0, round(Stress * 2, 0), 100)
ax2 = fig.add_subplot(212)
ax2.set_title('Weibull distribution')
ax2.plot(x, weibProbDist(x, slope, Stress), label='Prob. Density function')
ax2.fill_between(x, weibProbDist(x, slope, Stress), 0, color='0.90')
ax2.set_ylabel(r'$f(x, {\alpha}, {\beta})$')
ax2.set_xlabel(r'${\sigma} (MPa)$')
ax3 = ax2.twinx()
ax3.plot(x, weibCumDist(x, slope, Stress), 'r', label='Cum. distr. function')
ax3.fill_between(x, weibCumDist(x, slope, Stress), alpha=0.15, color='r')
ax3.set_ylim(0, 1)
ax3.set_ylabel(r'$F(x) = 1-exp(-(x/{\beta})^{\alpha})$')
ax3.legend(loc=7, prop={'size': 10})
ax2.legend(loc=2, prop={'size': 10})

#Text:
ax3.text(5, 0.7, r'${\alpha} = %.4f, {\beta} = %.4f$' % (m, Stress))
ax3.text(5, 0.6, r'$F(%s) = %s$' % (round(Stress, 3), round(Pr_Stress, 3)),
color='m')

#Arrows:
ax3.arrow(Stress, 0, 0, Pr_Stress - 0.1, head_width=2, head_length=0.1,
ec='m', fc='m')
ax3.arrow(Stress, Pr_Stress, round(Stress * 2, 0) - Stress - 10, 0,
head_width=0.02, head_length=10, ec='m', fc='m')

fig.subplots_adjust(hspace=0.5)
plt.show()
Weibull module:
"""
Weibull distribution:
    1) Probability density function
    2) Cumulative distribution function
    3) Mean
    4) Median
    5) Mode
    6) Variance - Standard deviation
    7) Skewness
    8) Gamma function (Lanczos approximation)
"""

import numpy as np
import cmath as cm


def weibProbDist(x, a, b):
    return (a / b) * (x / b) ** (a - 1) * np.exp(-(x / b) ** a)


def weibCumDist(x, a, b):
    return 1 - np.exp(-(x / b) ** a)


def weibMean(a, b):
    return b * gamma(1 + 1 / a)


def weibMedian(a, b):
    return b * np.log(2) ** (1 / a)


def weibMode(a, b):
    if a > 1:
        return b * ((a - 1) / a) ** (1 / a)
    else:
        return 0


def weibVar(a, b):
    return (b ** 2) * gamma(1 + 2 / a) - weibMean(a, b) ** 2


def weibSkew(a, b):
    ax_1 = gamma(1 + 3 / a) + 2 * (gamma(1 + 1 / a) ** 3)
    ax_2 = 3 * (gamma(1 + 2 / a) * gamma(1 + 1 / a))
    return (b ** 3) * (ax_1 - ax_2)

"""
Gamma function:

Lanczos approximation implementation.
Numerical Recipes in C (2nd ed. Cambridge University Press 1992)

"""
g = 5

p0 = float(1.000000000190015)
p1 = float(76.18009172947146)
p2 = float(-86.50532032941677)
p3 = float(24.01409824083091)
p4 = float(-1.231739572450155)
p5 = float(1.208650973866179e-3)
p6 = float(-5.395239384953e-6)

Pn_Coef = np.array([p0, p1, p2, p3, p4, p5, p6])


def gamma(z):
    z = complex(z)
    if z.real < 0.5:
        #Recursion method:
        return cm.pi / (cm.sin(cm.pi * z) * gamma(1 - z))
    else:
        ax_1 = cm.sqrt(cm.pi * 2) / z
        ax2_sum = Pn_Coef[0]
        for i in range(1, g + 2):
            ax2_sum += Pn_Coef[i] / (z + i)
        t = z + g + 0.5
        ax_3 = (t ** (z + 0.5)) * cm.exp(-t)
        return ((ax_1 * ax2_sum) * ax_3).real