-
import numpy as np
-import matplotlib.pyplot as plt
+import matplotlib.pyplot as plt
+import numpy as np
+import random
from sklearn.linear_model import Ridge
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import make_pipeline
+from sklearn.linear_model import LinearRegression
-def f(x):
- """ function to approximate by polynomial interpolation"""
- return x*x*x
-
-# generate points used to plot
-x_plot = np.linspace(0, 10, 100)
-
-# generate points and keep a subset of them
-x = np.linspace(0, 10, 100)
-rng = np.random.RandomState(0)
-rng.shuffle(x)
-x = np.sort(x[:20])
-y = f(x)
-# create matrix versions of these arrays
-X = x[:, np.newaxis]
-X_plot = x_plot[:, np.newaxis]
-
-colors = ['teal', 'yellowgreen', 'gold']
-lw = 2
-plt.plot(x_plot, f(x_plot), color='cornflowerblue', linewidth=lw,
- label="ground truth")
-plt.scatter(x, y, color='navy', s=30, marker='o', label="training points")
-
-for count, degree in enumerate([3, 4, 5]):
- model = make_pipeline(PolynomialFeatures(degree), Ridge())
- model.fit(X, y)
- y_plot = model.predict(X_plot)
- plt.plot(x_plot, y_plot, color=colors[count], linewidth=lw,
- label="degree %d" % degree)
-
-plt.legend(loc='lower left')
+x=np.linspace(0.02,0.98,200)
+noise = np.asarray(random.sample((range(200)),200))
+y=x**3*noise
+yn=x**3*100
+poly3 = PolynomialFeatures(degree=3)
+X = poly3.fit_transform(x[:,np.newaxis])
+clf3 = LinearRegression()
+clf3.fit(X,y)
+Xplot=poly3.fit_transform(x[:,np.newaxis])
+poly3_plot=plt.plot(x, clf3.predict(Xplot), label='Cubic Fit')
+plt.plot(x,yn, color='red', label="True Cubic")
+plt.scatter(x, y, label='Data', color='orange', s=15)
+plt.legend()
plt.show()
+
+def error(a):
+ for i in y:
+ err=(y-yn)/yn
+ return abs(np.sum(err))/len(err)
+
+print (error(y))
+Similarly, using R, we can perform similar studies. The following R code illustrates this.
Non-Linear Least squares in R
@@ -703,6 +644,9 @@ text(0, 0.5
+
+In our lectures on regression analysis (and other ones as well), we will discuss in more details various R functionalities.
+
Another useful Python package is
pandas, which is an open source library
@@ -721,12 +665,149 @@ display(data_pandas)
-
Predator-Prey model from ecology
+Examples
+We present here several examples, with pertinent Python codes that we
+will us to illustrate various machine learning methods and ways to
+analyze, from simple to complex, various data sets. Many of these
+examples allow us to generate the data we want to analyze, following
+much of the same philosophy we discussed above when
+fitting various polynomials.
+
+
+We start with a simple exponential growth model that is meant to mimick an ecoli lab experiment.
+We can easily model this system and then produce the data used to train various machine learning algorithms.
+Another model from the life sciences is the so-called predator-prey model from ecology. Thereafter we present
+a simple model for financial transactions before moving to a random walk model and ending with
+the simulation of velocities of a non-interacting atom or molecule confined to move in a one-dimensional region.
+
+
Ecoli lab experiment
+
+
+A typical pattern seen in population models is that the population grows faster and faster. Why? Is there an underlying (general) mechanism?
+Here we will construct a model for cell growth based on a simple difference equation for the growth. We make the following assumptions
+
+
+- Cells divide after \( T \) seconds on average (one generation)
+- \( 2N \) celles divide into twice as many new cells \( \Delta N \) in a time
+ interval \( \Delta t \) as \( N \) cells would: \( \Delta N \propto N \)
+- \( N \) cells result in twice as many new individuals \( \Delta N \) in
+ time \( 2\Delta t \) as in time \( \Delta t \): \( \Delta N \propto\Delta t \)
+- Same proportionality wrt death
+- Proposed model: \( \Delta N = b\Delta t N - d\Delta tN \) for some unknown
+ constants \( b \) (births) and \( d \) (deaths)
+- Describe evolution in discrete time: \( t_n=n\Delta t \)
+- Program-friendly notation: \( N \) at \( t_n \) is \( N^n \)
+- Math model: \( N^{n+1} = N^n + r\Delta t\, N \) (with \( \ r=b-d \))
+- Program model:
N[n+1] = N[n] + r*dt*N[n]
+
+
+
+
+
+
+The difference equation can be programmed in a simple was, and in order to get started we
+set \( r=1.5 \), \( N^0=1 \), \( \Delta t=0.5 \). The program reads
+
+
+
+
+
import numpy as np
+
+t = np.linspace(0, 10, 21) # 20 intervals in [0, 10]
+dt = t[1] - t[0]
+N = np.zeros(t.size)
+
+N[0] = 1
+r = 0.5
+
+for n in range(0, N.size-1, 1):
+ N[n+1] = N[n] + r*dt*N[n]
+ print 'N[%d]=%.1f' % (n+1, N[n+1])
+
+
+and it generates the following output
+
+
+
+
N[1]=1.2
+N[2]=1.6
+N[3]=2.0
+N[4]=2.4
+N[5]=3.1
+N[6]=3.8
+N[7]=4.8
+N[8]=6.0
+N[9]=7.5
+N[10]=9.3
+N[11]=11.6
+N[12]=14.6
+N[13]=18.2
+N[14]=22.7
+N[15]=28.4
+N[16]=35.5
+N[17]=44.4
+N[18]=55.5
+N[19]=69.4
+N[20]=86.7
+
+
+This forms our data which later will define our training set.
+In this case we defined the value of the parameter \( r \). We could alternatively assume that we just received the
+above data file and where asked to use find \( r \). How can we estimate \( r \) from data?
+
+
+We can use the difference equation with the experimental data
+$$ N^{n+1} = N^n + r\Delta t N^n$$
+
+Suppose now that \( N^{n+1} \) and \( N^n \) are known from data. Then we could solve with respect to \( r \) as follows
+$$ r = \frac{N^{n+1}-N^n}{N^n\Delta t} $$
+
+Suppose we set \( t_1=600 \), \( t_2=1200 \),
+\( N^1=140 \) and \( N^2=250 \). We obtain then \( r=0.0013 \). The exact value is \( r = 0.000694 \)
+The following code plot
+
+
+
+
import numpy as np
+
+# Estimate r
+data = np.loadtxt('ecoli.csv', delimiter=',')
+t_e = data[:,0]
+N_e = data[:,1]
+i = 2 # Data point (i,i+1) used to estimate r
+r = (N_e[i+1] - N_e[i])/(N_e[i]*(t_e[i+1] - t_e[i]))
+print 'Estimated r=%.5f' % r
+# Can experiment with r values and see if the model can
+# match the data better
+
+T = 1200 # cell can divide after T sec
+t_max = 5*T # 5 generations in experiment
+t = np.linspace(0, t_max, 1000)
+dt = t[1] - t[0]
+N = np.zeros(t.size)
+
+N[0] = 100
+for n in range(0, len(t)-1, 1):
+ N[n+1] = N[n] + r*dt*N[n]
+
+import matplotlib.pyplot as plt
+plt.plot(t, N, 'r-', t_e, N_e, 'bo')
+plt.xlabel('time [s]'); plt.ylabel('N')
+plt.legend(['model', 'experiment'], loc='upper left')
+plt.show()
+
+
+We can then change the parameter \( r \) in the program and play around to make a better fit. By now we know that this
+'search bythe eye' approach is not the most optimal one.
+
+
Predator-Prey model from ecology
+
+
The population dynamics of a simple predator-prey system is a
classical example shown in many biology textbooks when ecological
systems are discussed. The system contains all elements of the
@@ -739,19 +820,7 @@ scientific method:
trying to extract general behaviors and extract eventual laws or patterns
develop mathematical relations for the uncovered regularities/laws and test these by per forming new experiments
-- - -
-
Lots of data about populations of hares and lynx collected from furs in Hudson Bay, Canada, are available. It is known that the populations oscillate. Why? Here we start by @@ -761,21 +830,7 @@ Here we start by
- - -
-
- -
Most mammalian predators rely on a variety of prey, which complicates mathematical modeling; however, a few predators have become highly specialized and seek almost exclusively a single prey species. An example of this simplified predator-prey interaction is seen in Canadian northern forests, where the populations of the lynx and the snowshoe hare are intertwined in a life and death struggle.
@@ -816,20 +871,6 @@ One reason that this particular system has been so extensively studied is that t
- - - - -
- - -
-
-
import numpy as np
@@ -852,26 +893,10 @@ plt.savefig(
plt.savefig('Hudson_Bay_data.png')
plt.show()
-
- - -
- - -
We see from the plot that there are indeed fluctuations. We would like to create a mathematical model that explains these population fluctuations. Ecologists have predicted that in a simple @@ -892,22 +917,9 @@ climate and other complicating factors. How significant are these?
- - -
-
The classical way (in all books) is to present the Lotka-Volterra equations: $$ @@ -926,43 +938,6 @@ Here,
- - -
- -
- - -
-
The population of hares evolves due to births and deaths exactly as a bacteria population: $$ @@ -979,19 +954,8 @@ So in fraction \( b\Delta t HL \), the lynx eat hares. This loss of hares must be accounted for. Subtracted in the equation for hares: $$ \Delta H = a\Delta t H^n - b \Delta t H^nL^n$$ -
- - -
-
We assume that the primary growth for the lynx population depends on sufficient food for raising lynx kittens, which implies an adequate source of nutrients from predation on hares. Thus, the growth of the lynx population does not only depend of how many lynx there are, but on how many hares they can eat. In a time interval \( \Delta t HL \) hares and lynx can meet, and in a fraction \( b\Delta t HL \) the lynx eats the hare. All of this does not @@ -1000,29 +964,10 @@ contribute to the growth of lynx, again just a fraction of \( d\Delta t HL \). In addition, lynx die just as in the population dynamics with one isolated animal population, leading to a loss \( -c\Delta t L \). -
-
The accounting of lynx then looks like $$ \Delta L = d\Delta t H^nL^n - c\Delta t L^n$$ -
- - -
-
By writing up the definition of \( \Delta H \) and \( \Delta L \), and putting all assumed known terms \( H^n \) and \( L^n \) on the right-hand side, we have @@ -1036,44 +981,16 @@ Note:
- - -
-
- -
- - -
@@ -1129,24 +1046,19 @@ plt.show()
- - -
-If we perform a least-square fitting, we can find optimal values for the parameters \( a \), \( b \), \( d \), \( c \). The optimal parameters are \( a=0.4807 \), \( b=0.02482 \), \( d=0.9272 \) and \( c=0.02756 \). These parameters result in a slightly modified initial conditions, namely \( H(0) = 34.91 \) and \( L(0)=3.857 \). With these parameters we are now ready to solve the equations and plot these data together with the experimental values. +We will later perform a least-square fitting. Then we can find optimal +values for the parameters \( a \), \( b \), \( d \), \( c \). In our calculations here +we set \( a=0.4807 \), \( b=0.02482 \), \( d=0.9272 \) and \( c=0.02756 \). These +parameters result in a slightly modified initial conditions, namely +\( H(0) = 34.91 \) and \( L(0)=3.857 \).
- - -
+The following Python demonstrates how we can use linear regression to fit for example the population of lynx. +Similarly, we have also used a decision tree algorithm to fit the lynx population data. As expected, the linear regression is not exactly impressive
@@ -1170,17 +1082,7 @@ plt.plot(x, y, label.show()
-
- - -
+The similar code for linear regression in R reads
@@ -1204,198 +1106,8 @@ plot(linearMod) confint(linearMod) predict(linearMod,data.frame(Year=c(1910,1914,1920)),interval="confidence")
-
- - -
-
-The population grows faster and faster. Why? Is there an underlying (general) mechanism? -
- -
N[n+1] = N[n] + r*dt*N[n]- - -
-
-Let us solve the difference equation in as simple way as possible, -just to train some programming: \( r=1.5 \), \( N^0=1 \), \( \Delta t=0.5 \) - -
- - -
import numpy as np
-
-t = np.linspace(0, 10, 21) # 20 intervals in [0, 10]
-dt = t[1] - t[0]
-N = np.zeros(t.size)
-
-N[0] = 1
-r = 0.5
-
-for n in range(0, N.size-1, 1):
- N[n+1] = N[n] + r*dt*N[n]
- print 'N[%d]=%.1f' % (n+1, N[n+1])
--
-% if FORMAT != 'ipynb': - - -
- - -
N[1]=1.2
-N[2]=1.6
-N[3]=2.0
-N[4]=2.4
-N[5]=3.1
-N[6]=3.8
-N[7]=4.8
-N[8]=6.0
-N[9]=7.5
-N[10]=9.3
-N[11]=11.6
-N[12]=14.6
-N[13]=18.2
-N[14]=22.7
-N[15]=28.4
-N[16]=35.5
-N[17]=44.4
-N[18]=55.5
-N[19]=69.4
-N[20]=86.7
--% endif - -
- - -
-
- -
-Use experimental data in the fraction, say \( t_1=600 \), \( t_2=1200 \), -\( N^1=140 \), \( N^2=250 \): \( r=0.0013 \). -
- - -
- - -
-
-
- - -
import numpy as np
-
-# Estimate r
-data = np.loadtxt('ecoli.csv', delimiter=',')
-t_e = data[:,0]
-N_e = data[:,1]
-i = 2 # Data point (i,i+1) used to estimate r
-r = (N_e[i+1] - N_e[i])/(N_e[i]*(t_e[i+1] - t_e[i]))
-print 'Estimated r=%.5f' % r
-# Can experiment with r values and see if the model can
-# match the data better
-
-T = 1200 # cell can divide after T sec
-t_max = 5*T # 5 generations in experiment
-t = np.linspace(0, t_max, 1000)
-dt = t[1] - t[0]
-N = np.zeros(t.size)
-
-N[0] = 100
-for n in range(0, len(t)-1, 1):
- N[n+1] = N[n] + r*dt*N[n]
-
-import matplotlib.pyplot as plt
-plt.plot(t, N, 'r-', t_e, N_e, 'bo')
-plt.xlabel('time [s]'); plt.ylabel('N')
-plt.legend(['model', 'experiment'], loc='upper left')
-plt.show()
-
-Change r in the program and play around to make a better fit!
-
- - -
The aim here is to simulate financial transactions among financial agents @@ -1469,8 +1181,7 @@ exponentially decreases with \( m' \). We assume that we have \( N=500 \) agents. In each simulation, we need a sufficiently large number of transactions, say \( 10^7 \). Our aim is find the final equilibrium distribution \( w_m \). In order to do that we would need several runs of the above simulations, at least \( 10^3-10^4 \) runs (experiments). -
Our task is to first set up an algorithm which simulates the above transactions with an initial amount \( m_0 \). The challenge here is to figure out a Monte Carlo simulation based on the @@ -1569,40 +1280,11 @@ $$ showing how money is conserved during a transaction. Select values of \( \lambda =0.25,0.5 \) and \( \lambda=0.9 \) and try to extract the corresponding - equilibrium distributions and compare these with the Gibbs distribution. Comment your results. -Extract a parametrization of the above curves, see for example Patriarca and collaborators and see if you can parametrize the high-end tails of the distributions in terms of power laws. Comment your results. + equilibrium distributions and compare these with the Gibbs distribution. We will use this model to +extract a parametrization of the above curves, see for example Patriarca and collaborators. -
-In the studies above the agents were selected randomly, irrespective of whether we allowed for -saving or not during a transaction. What is often observed is that various agents tend to make preferences for for whom to interact with. We will now study the evolution of the distribution of wealth \( w_m \) by assuming that there is a likelihood -$$ -p_{ij} \propto \vert m_i-m_j\vert^{-\alpha}, -$$ +
-Perform the analysis with and without a saving \( \lambda \) on each transaction and comment your results. -We add to the previous probability the possibility that two agents who interact have performed similar transactions earlier. That is, in addition to being financially close, we assume that the likelihood for interacting increases if two agents have interacted earlier. -We add this feature by modifying the previous likelihood to -$$ -p_{ij} \propto \vert m_i-m_j\vert^{-\alpha}\left(c_{ij}+1\right)^{\gamma}, -$$ - -where \( c_{ij} \) represents the number of previous interactions that have taken place between \( i \) and \( j \). The factor \( 1 \) is added in order to ensure that if they have not interacted earlier they can still interact. Perform similar studies as above with \( N=1000 \), \( \alpha=1.0 \) and \( \alpha=2.0 \) using \( \gamma = 0.0, 1.0, 2.0, 3.0 \) and \( 4.0 \). Plot the wealth distributions for these cases and try to extract eventual power law tails with and without a saving \( \lambda \) in each transaction. Comment your results and compare them with figures 5 and 6 of Goswami and Sen. - -
- - -
diff --git a/doc/pub/How2ReadData/html/How2ReadData-reveal.html b/doc/pub/How2ReadData/html/How2ReadData-reveal.html index 519fc4df6..d0879d153 100644 --- a/doc/pub/How2ReadData/html/How2ReadData-reveal.html +++ b/doc/pub/How2ReadData/html/How2ReadData-reveal.html @@ -611,46 +611,39 @@ a linear \( x \)-dependence we study now a cubic polynomial and use the polynomi
-
import numpy as np
-import matplotlib.pyplot as plt
+import matplotlib.pyplot as plt
+import numpy as np
+import random
from sklearn.linear_model import Ridge
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import make_pipeline
+from sklearn.linear_model import LinearRegression
-def f(x):
- """ function to approximate by polynomial interpolation"""
- return x*x*x
-
-# generate points used to plot
-x_plot = np.linspace(0, 10, 100)
-
-# generate points and keep a subset of them
-x = np.linspace(0, 10, 100)
-rng = np.random.RandomState(0)
-rng.shuffle(x)
-x = np.sort(x[:20])
-y = f(x)
-# create matrix versions of these arrays
-X = x[:, np.newaxis]
-X_plot = x_plot[:, np.newaxis]
-
-colors = ['teal', 'yellowgreen', 'gold']
-lw = 2
-plt.plot(x_plot, f(x_plot), color='cornflowerblue', linewidth=lw,
- label="ground truth")
-plt.scatter(x, y, color='navy', s=30, marker='o', label="training points")
-
-for count, degree in enumerate([3, 4, 5]):
- model = make_pipeline(PolynomialFeatures(degree), Ridge())
- model.fit(X, y)
- y_plot = model.predict(X_plot)
- plt.plot(x_plot, y_plot, color=colors[count], linewidth=lw,
- label="degree %d" % degree)
-
-plt.legend(loc='lower left')
+x=np.linspace(0.02,0.98,200)
+noise = np.asarray(random.sample((range(200)),200))
+y=x**3*noise
+yn=x**3*100
+poly3 = PolynomialFeatures(degree=3)
+X = poly3.fit_transform(x[:,np.newaxis])
+clf3 = LinearRegression()
+clf3.fit(X,y)
+Xplot=poly3.fit_transform(x[:,np.newaxis])
+poly3_plot=plt.plot(x, clf3.predict(Xplot), label='Cubic Fit')
+plt.plot(x,yn, color='red', label="True Cubic")
+plt.scatter(x, y, label='Data', color='orange', s=15)
+plt.legend()
plt.show()
+
+def error(a):
+ for i in y:
+ err=(y-yn)/yn
+ return abs(np.sum(err))/len(err)
+
+print (error(y))
+
+Similarly, using R, we can perform similar studies. The following R code illustrates this.
@@ -684,6 +677,9 @@ text(0, 0.5
+
+In our lectures on regression analysis (and other ones as well), we will discuss in more details various R functionalities.
+
Another useful Python package is
pandas, which is an open source library
@@ -703,11 +699,148 @@ display(data_pandas)
-Predator-Prey model from ecology
+Examples
+We present here several examples, with pertinent Python codes that we
+will us to illustrate various machine learning methods and ways to
+analyze, from simple to complex, various data sets. Many of these
+examples allow us to generate the data we want to analyze, following
+much of the same philosophy we discussed above when
+fitting various polynomials.
+
+
+We start with a simple exponential growth model that is meant to mimick an ecoli lab experiment.
+We can easily model this system and then produce the data used to train various machine learning algorithms.
+Another model from the life sciences is the so-called predator-prey model from ecology. Thereafter we present
+a simple model for financial transactions before moving to a random walk model and ending with
+the simulation of velocities of a non-interacting atom or molecule confined to move in a one-dimensional region.
+
+
Ecoli lab experiment
+
+
+A typical pattern seen in population models is that the population grows faster and faster. Why? Is there an underlying (general) mechanism?
+Here we will construct a model for cell growth based on a simple difference equation for the growth. We make the following assumptions
+
+- Cells divide after \( T \) seconds on average (one generation)
+- \( 2N \) celles divide into twice as many new cells \( \Delta N \) in a time
+ interval \( \Delta t \) as \( N \) cells would: \( \Delta N \propto N \)
+- \( N \) cells result in twice as many new individuals \( \Delta N \) in
+ time \( 2\Delta t \) as in time \( \Delta t \): \( \Delta N \propto\Delta t \)
+- Same proportionality wrt death
+- Proposed model: \( \Delta N = b\Delta t N - d\Delta tN \) for some unknown
+ constants \( b \) (births) and \( d \) (deaths)
+- Describe evolution in discrete time: \( t_n=n\Delta t \)
+- Program-friendly notation: \( N \) at \( t_n \) is \( N^n \)
+- Math model: \( N^{n+1} = N^n + r\Delta t\, N \) (with \( \ r=b-d \))
+- Program model:
N[n+1] = N[n] + r*dt*N[n]
+
+
+
+
+The difference equation can be programmed in a simple was, and in order to get started we
+set \( r=1.5 \), \( N^0=1 \), \( \Delta t=0.5 \). The program reads
+
+
+
+
+
import numpy as np
+
+t = np.linspace(0, 10, 21) # 20 intervals in [0, 10]
+dt = t[1] - t[0]
+N = np.zeros(t.size)
+
+N[0] = 1
+r = 0.5
+
+for n in range(0, N.size-1, 1):
+ N[n+1] = N[n] + r*dt*N[n]
+ print 'N[%d]=%.1f' % (n+1, N[n+1])
+
+
+and it generates the following output
+
+
+
+
N[1]=1.2
+N[2]=1.6
+N[3]=2.0
+N[4]=2.4
+N[5]=3.1
+N[6]=3.8
+N[7]=4.8
+N[8]=6.0
+N[9]=7.5
+N[10]=9.3
+N[11]=11.6
+N[12]=14.6
+N[13]=18.2
+N[14]=22.7
+N[15]=28.4
+N[16]=35.5
+N[17]=44.4
+N[18]=55.5
+N[19]=69.4
+N[20]=86.7
+
+
+This forms our data which later will define our training set.
+In this case we defined the value of the parameter \( r \). We could alternatively assume that we just received the
+above data file and where asked to use find \( r \). How can we estimate \( r \) from data?
+
+
+We can use the difference equation with the experimental data
+
+$$ N^{n+1} = N^n + r\Delta t N^n$$
+
+
+Suppose now that \( N^{n+1} \) and \( N^n \) are known from data. Then we could solve with respect to \( r \) as follows
+
+$$ r = \frac{N^{n+1}-N^n}{N^n\Delta t} $$
+
+
+Suppose we set \( t_1=600 \), \( t_2=1200 \),
+\( N^1=140 \) and \( N^2=250 \). We obtain then \( r=0.0013 \). The exact value is \( r = 0.000694 \)
+The following code plot
+
+
+
+
import numpy as np
+
+# Estimate r
+data = np.loadtxt('ecoli.csv', delimiter=',')
+t_e = data[:,0]
+N_e = data[:,1]
+i = 2 # Data point (i,i+1) used to estimate r
+r = (N_e[i+1] - N_e[i])/(N_e[i]*(t_e[i+1] - t_e[i]))
+print 'Estimated r=%.5f' % r
+# Can experiment with r values and see if the model can
+# match the data better
+
+T = 1200 # cell can divide after T sec
+t_max = 5*T # 5 generations in experiment
+t = np.linspace(0, t_max, 1000)
+dt = t[1] - t[0]
+N = np.zeros(t.size)
+
+N[0] = 100
+for n in range(0, len(t)-1, 1):
+ N[n+1] = N[n] + r*dt*N[n]
+
+import matplotlib.pyplot as plt
+plt.plot(t, N, 'r-', t_e, N_e, 'bo')
+plt.xlabel('time [s]'); plt.ylabel('N')
+plt.legend(['model', 'experiment'], loc='upper left')
+plt.show()
+
+
+We can then change the parameter \( r \) in the program and play around to make a better fit. By now we know that this
+'search bythe eye' approach is not the most optimal one.
+
+
Predator-Prey model from ecology
+
The population dynamics of a simple predator-prey system is a
classical example shown in many biology textbooks when ecological
@@ -721,17 +854,8 @@ scientific method:
trying to extract general behaviors and extract eventual laws or patterns
develop mathematical relations for the uncovered regularities/laws and test these by per forming new experiments
-
-
+ Lots of data about populations of hares and lynx collected from furs in Hudson Bay, Canada, are available. It is known that the populations oscillate. Why? Here we start by @@ -741,17 +865,8 @@ Here we start by
-
+ Most mammalian predators rely on a variety of prey, which complicates mathematical modeling; however, a few predators have become highly specialized and seek almost exclusively a single prey species. An example of this simplified predator-prey interaction is seen in Canadian northern forests, where the populations of the lynx and the snowshoe hare are intertwined in a life and death struggle.
@@ -786,17 +901,6 @@ One reason that this particular system has been so extensively studied is that t
-
@@ -820,23 +924,9 @@ plt.savefig('Hudson_Bay_data.pdf') plt.savefig('Hudson_Bay_data.png') plt.show()
We see from the plot that there are indeed fluctuations. We would like to create a mathematical model that explains these @@ -858,20 +948,10 @@ climate and other complicating factors. How significant are these?
-
+ The classical way (in all books) is to present the Lotka-Volterra equations:
@@ -893,37 +973,6 @@ Here,
-Most books quickly establish the model and then use considerable space on -discussing the qualitative properties of this nonlinear system of -ODEs (which cannot be solved) -
-
The population of hares evolves due to births and deaths exactly as a bacteria population:
@@ -944,16 +993,7 @@ loss of hares must be accounted for. Subtracted in the equation for hares:
$$ \Delta H = a\Delta t H^n - b \Delta t H^nL^n$$
-
-
We assume that the primary growth for the lynx population depends on sufficient food for raising lynx kittens, which implies an adequate source of nutrients from predation on hares. Thus, the growth of the lynx population does not only depend of how many lynx there are, but on how many hares they can eat. In a time interval \( \Delta t HL \) hares and lynx can meet, and in a @@ -963,26 +1003,11 @@ contribute to the growth of lynx, again just a fraction of \( d\Delta t HL \). In addition, lynx die just as in the population dynamics with one isolated animal population, leading to a loss \( -c\Delta t L \). -
-
The accounting of lynx then looks like
$$ \Delta L = d\Delta t H^nL^n - c\Delta t L^n$$
-
-
By writing up the definition of \( \Delta H \) and \( \Delta L \), and putting all assumed known terms \( H^n \) and \( L^n \) on the right-hand side, we have @@ -1000,38 +1025,17 @@ Note:
-
+
@@ -1083,24 +1087,20 @@ plt.show()
-If we perform a least-square fitting, we can find optimal values for the parameters \( a \), \( b \), \( d \), \( c \). The optimal parameters are \( a=0.4807 \), \( b=0.02482 \), \( d=0.9272 \) and \( c=0.02756 \). These parameters result in a slightly modified initial conditions, namely \( H(0) = 34.91 \) and \( L(0)=3.857 \). With these parameters we are now ready to solve the equations and plot these data together with the experimental values. -
+The following Python demonstrates how we can use linear regression to fit for example the population of lynx. +Similarly, we have also used a decision tree algorithm to fit the lynx population data. As expected, the linear regression is not exactly impressive
@@ -1123,15 +1123,8 @@ plt.plot(line, regline.predict(line), label= " plt.plot(x, y, label= "Linear Regression") plt.show()
+The similar code for linear regression in R reads
@@ -1156,185 +1149,7 @@ confint(linearMod) predict(linearMod,data.frame(Year=c(1910,1914,1920)),interval="confidence")
-
-The population grows faster and faster. Why? Is there an underlying (general) mechanism? -
N[n+1] = N[n] + r*dt*N[n]-
-Let us solve the difference equation in as simple way as possible, -just to train some programming: \( r=1.5 \), \( N^0=1 \), \( \Delta t=0.5 \) - -
- - -
import numpy as np
-
-t = np.linspace(0, 10, 21) # 20 intervals in [0, 10]
-dt = t[1] - t[0]
-N = np.zeros(t.size)
-
-N[0] = 1
-r = 0.5
-
-for n in range(0, N.size-1, 1):
- N[n+1] = N[n] + r*dt*N[n]
- print 'N[%d]=%.1f' % (n+1, N[n+1])
--% if FORMAT != 'ipynb': -
- - -
N[1]=1.2
-N[2]=1.6
-N[3]=2.0
-N[4]=2.4
-N[5]=3.1
-N[6]=3.8
-N[7]=4.8
-N[8]=6.0
-N[9]=7.5
-N[10]=9.3
-N[11]=11.6
-N[12]=14.6
-N[13]=18.2
-N[14]=22.7
-N[15]=28.4
-N[16]=35.5
-N[17]=44.4
-N[18]=55.5
-N[19]=69.4
-N[20]=86.7
--% endif -
-
- -We can use the difference equation with the experimental data - -
-$$ N^{n+1} = N^n + r\Delta t N^n$$
-
-
-Say \( N^{n+1} \) and \( N^n \) are known from data, solve wrt \( r \):
-
-
-$$ r = \frac{N^{n+1}-N^n}{N^n\Delta t} $$
-
-
-
-Use experimental data in the fraction, say \( t_1=600 \), \( t_2=1200 \), -\( N^1=140 \), \( N^2=250 \): \( r=0.0013 \). -
- - -
-
- - -
import numpy as np
-
-# Estimate r
-data = np.loadtxt('ecoli.csv', delimiter=',')
-t_e = data[:,0]
-N_e = data[:,1]
-i = 2 # Data point (i,i+1) used to estimate r
-r = (N_e[i+1] - N_e[i])/(N_e[i]*(t_e[i+1] - t_e[i]))
-print 'Estimated r=%.5f' % r
-# Can experiment with r values and see if the model can
-# match the data better
-
-T = 1200 # cell can divide after T sec
-t_max = 5*T # 5 generations in experiment
-t = np.linspace(0, t_max, 1000)
-dt = t[1] - t[0]
-N = np.zeros(t.size)
-
-N[0] = 100
-for n in range(0, len(t)-1, 1):
- N[n+1] = N[n] + r*dt*N[n]
-
-import matplotlib.pyplot as plt
-plt.plot(t, N, 'r-', t_e, N_e, 'bo')
-plt.xlabel('time [s]'); plt.ylabel('N')
-plt.legend(['model', 'experiment'], loc='upper left')
-plt.show()
-
-Change r in the program and play around to make a better fit!
-
The aim here is to simulate financial transactions among financial agents @@ -1420,8 +1235,7 @@ exponentially decreases with \( m' \). We assume that we have \( N=500 \) agents. In each simulation, we need a sufficiently large number of transactions, say \( 10^7 \). Our aim is find the final equilibrium distribution \( w_m \). In order to do that we would need several runs of the above simulations, at least \( 10^3-10^4 \) runs (experiments). -
Our task is to first set up an algorithm which simulates the above transactions with an initial amount \( m_0 \). The challenge here is to figure out a Monte Carlo simulation based on the @@ -1530,46 +1344,11 @@ $$ showing how money is conserved during a transaction. Select values of \( \lambda =0.25,0.5 \) and \( \lambda=0.9 \) and try to extract the corresponding - equilibrium distributions and compare these with the Gibbs distribution. Comment your results. -Extract a parametrization of the above curves, see for example Patriarca and collaborators and see if you can parametrize the high-end tails of the distributions in terms of power laws. Comment your results. + equilibrium distributions and compare these with the Gibbs distribution. We will use this model to +extract a parametrization of the above curves, see for example Patriarca and collaborators. -
-In the studies above the agents were selected randomly, irrespective of whether we allowed for -saving or not during a transaction. What is often observed is that various agents tend to make preferences for for whom to interact with. We will now study the evolution of the distribution of wealth \( w_m \) by assuming that there is a likelihood -
-$$
-p_{ij} \propto \vert m_i-m_j\vert^{-\alpha},
-$$
-
+
-$$
-w_m\propto m^{-1-\alpha}.
-$$
-
-
-What happens if \( \alpha \gg 1 \)?
-
-
-Perform the analysis with and without a saving \( \lambda \) on each transaction and comment your results. -We add to the previous probability the possibility that two agents who interact have performed similar transactions earlier. That is, in addition to being financially close, we assume that the likelihood for interacting increases if two agents have interacted earlier. -We add this feature by modifying the previous likelihood to -
-$$
-p_{ij} \propto \vert m_i-m_j\vert^{-\alpha}\left(c_{ij}+1\right)^{\gamma},
-$$
-
-
-where \( c_{ij} \) represents the number of previous interactions that have taken place between \( i \) and \( j \). The factor \( 1 \) is added in order to ensure that if they have not interacted earlier they can still interact. Perform similar studies as above with \( N=1000 \), \( \alpha=1.0 \) and \( \alpha=2.0 \) using \( \gamma = 0.0, 1.0, 2.0, 3.0 \) and \( 4.0 \). Plot the wealth distributions for these cases and try to extract eventual power law tails with and without a saving \( \lambda \) in each transaction. Comment your results and compare them with figures 5 and 6 of Goswami and Sen.
-
diff --git a/doc/pub/How2ReadData/html/How2ReadData-solarized.html b/doc/pub/How2ReadData/html/How2ReadData-solarized.html index fead58baf..dd677c506 100644 --- a/doc/pub/How2ReadData/html/How2ReadData-solarized.html +++ b/doc/pub/How2ReadData/html/How2ReadData-solarized.html @@ -70,46 +70,14 @@ div { text-align: justify; text-justify: inter-word; } None, '___sec5'), ('Non-Linear Least squares in R', 2, None, '___sec6'), - ('Predator-Prey model from ecology', 2, None, '___sec7'), - ('Case study from Hudson bay', 2, None, '___sec8'), - ('Hudson bay data', 2, None, '___sec9'), - ('Plotting the data', 2, None, '___sec10'), - ('Hares and lynx in Hudson bay from 1900 to 1920', - 2, + ('Examples', 2, None, '___sec7'), + ('Ecoli lab experiment', 3, None, '___sec8'), + ('Predator-Prey model from ecology', 3, None, '___sec9'), + ('Simulating financial transactions', 3, None, '___sec10'), + ('Particle in one dimension and velocity distribution', + 3, None, - '___sec11'), - ('Why now create a computer model for the hare and lynx ' - 'populations?', - 2, - None, - '___sec12'), - ('The traditional (top-down) approach', 2, None, '___sec13'), - ('Basic mathematics notation', 2, None, '___sec14'), - ('Basic dynamics of the population of hares', - 2, - None, - '___sec15'), - ('Basic dynamics of the population of lynx', 2, None, '___sec16'), - ('Evolution equations', 2, None, '___sec17'), - ('Adapt the model to the Hudson Bay case', 2, None, '___sec18'), - ('The program', 2, None, '___sec19'), - ('The plot', 2, None, '___sec20'), - ('Linear regression in Python', 2, None, '___sec21'), - ('Linear Least squares in R', 2, None, '___sec22'), - ('Example: ecoli lab experiment', 2, None, '___sec23'), - ('The program', 2, None, '___sec24'), - ('The output', 2, None, '___sec25'), - ('Parameter estimation', 2, None, '___sec26'), - ('A program relevant for the biological problem', - 2, - None, - '___sec27'), - ('Simulating financial transactions', 2, None, '___sec28'), - ('Simulation of Transactions', 3, None, '___sec29'), - ('Particle in one dimension an velocity distribution', - 2, - None, - '___sec30')]} + '___sec11')]} end of tocinfo -->
@@ -585,47 +553,39 @@ a linear \( x \)-dependence we study now a cubic polynomial and use the polynomi-
import numpy as np
-import matplotlib.pyplot as plt
+import matplotlib.pyplot as plt
+import numpy as np
+import random
from sklearn.linear_model import Ridge
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import make_pipeline
+from sklearn.linear_model import LinearRegression
-def f(x):
- """ function to approximate by polynomial interpolation"""
- return x*x*x
-
-# generate points used to plot
-x_plot = np.linspace(0, 10, 100)
-
-# generate points and keep a subset of them
-x = np.linspace(0, 10, 100)
-rng = np.random.RandomState(0)
-rng.shuffle(x)
-x = np.sort(x[:20])
-y = f(x)
-# create matrix versions of these arrays
-X = x[:, np.newaxis]
-X_plot = x_plot[:, np.newaxis]
-
-colors = ['teal', 'yellowgreen', 'gold']
-lw = 2
-plt.plot(x_plot, f(x_plot), color='cornflowerblue', linewidth=lw,
- label="ground truth")
-plt.scatter(x, y, color='navy', s=30, marker='o', label="training points")
-
-for count, degree in enumerate([3, 4, 5]):
- model = make_pipeline(PolynomialFeatures(degree), Ridge())
- model.fit(X, y)
- y_plot = model.predict(X_plot)
- plt.plot(x_plot, y_plot, color=colors[count], linewidth=lw,
- label="degree %d" % degree)
-
-plt.legend(loc='lower left')
+x=np.linspace(0.02,0.98,200)
+noise = np.asarray(random.sample((range(200)),200))
+y=x**3*noise
+yn=x**3*100
+poly3 = PolynomialFeatures(degree=3)
+X = poly3.fit_transform(x[:,np.newaxis])
+clf3 = LinearRegression()
+clf3.fit(X,y)
+Xplot=poly3.fit_transform(x[:,np.newaxis])
+poly3_plot=plt.plot(x, clf3.predict(Xplot), label='Cubic Fit')
+plt.plot(x,yn, color='red', label="True Cubic")
+plt.scatter(x, y, label='Data', color='orange', s=15)
+plt.legend()
plt.show()
+
+def error(a):
+ for i in y:
+ err=(y-yn)/yn
+ return abs(np.sum(err))/len(err)
+
+print (error(y))
+Similarly, using R, we can perform similar studies. The following R code illustrates this.
Non-Linear Least squares in R
@@ -659,6 +619,9 @@ text(0, 0.5
+
+In our lectures on regression analysis (and other ones as well), we will discuss in more details various R functionalities.
+
Another useful Python package is
pandas, which is an open source library
@@ -677,11 +640,147 @@ display(data_pandas)
-
Predator-Prey model from ecology
+Examples
+We present here several examples, with pertinent Python codes that we
+will us to illustrate various machine learning methods and ways to
+analyze, from simple to complex, various data sets. Many of these
+examples allow us to generate the data we want to analyze, following
+much of the same philosophy we discussed above when
+fitting various polynomials.
+
+
+We start with a simple exponential growth model that is meant to mimick an ecoli lab experiment.
+We can easily model this system and then produce the data used to train various machine learning algorithms.
+Another model from the life sciences is the so-called predator-prey model from ecology. Thereafter we present
+a simple model for financial transactions before moving to a random walk model and ending with
+the simulation of velocities of a non-interacting atom or molecule confined to move in a one-dimensional region.
+
+
Ecoli lab experiment
+
+
+A typical pattern seen in population models is that the population grows faster and faster. Why? Is there an underlying (general) mechanism?
+Here we will construct a model for cell growth based on a simple difference equation for the growth. We make the following assumptions
+
+
+
+- Cells divide after \( T \) seconds on average (one generation)
+- \( 2N \) celles divide into twice as many new cells \( \Delta N \) in a time
+ interval \( \Delta t \) as \( N \) cells would: \( \Delta N \propto N \)
+- \( N \) cells result in twice as many new individuals \( \Delta N \) in
+ time \( 2\Delta t \) as in time \( \Delta t \): \( \Delta N \propto\Delta t \)
+- Same proportionality wrt death
+- Proposed model: \( \Delta N = b\Delta t N - d\Delta tN \) for some unknown
+ constants \( b \) (births) and \( d \) (deaths)
+- Describe evolution in discrete time: \( t_n=n\Delta t \)
+- Program-friendly notation: \( N \) at \( t_n \) is \( N^n \)
+- Math model: \( N^{n+1} = N^n + r\Delta t\, N \) (with \( \ r=b-d \))
+- Program model:
N[n+1] = N[n] + r*dt*N[n]
+
+
+
+
+
+The difference equation can be programmed in a simple was, and in order to get started we
+set \( r=1.5 \), \( N^0=1 \), \( \Delta t=0.5 \). The program reads
+
+
+
+
+
import numpy as np
+
+t = np.linspace(0, 10, 21) # 20 intervals in [0, 10]
+dt = t[1] - t[0]
+N = np.zeros(t.size)
+
+N[0] = 1
+r = 0.5
+
+for n in range(0, N.size-1, 1):
+ N[n+1] = N[n] + r*dt*N[n]
+ print 'N[%d]=%.1f' % (n+1, N[n+1])
+
+
+and it generates the following output
+
+
+
+
N[1]=1.2
+N[2]=1.6
+N[3]=2.0
+N[4]=2.4
+N[5]=3.1
+N[6]=3.8
+N[7]=4.8
+N[8]=6.0
+N[9]=7.5
+N[10]=9.3
+N[11]=11.6
+N[12]=14.6
+N[13]=18.2
+N[14]=22.7
+N[15]=28.4
+N[16]=35.5
+N[17]=44.4
+N[18]=55.5
+N[19]=69.4
+N[20]=86.7
+
+
+This forms our data which later will define our training set.
+In this case we defined the value of the parameter \( r \). We could alternatively assume that we just received the
+above data file and where asked to use find \( r \). How can we estimate \( r \) from data?
+
+
+We can use the difference equation with the experimental data
+$$ N^{n+1} = N^n + r\Delta t N^n$$
+
+Suppose now that \( N^{n+1} \) and \( N^n \) are known from data. Then we could solve with respect to \( r \) as follows
+$$ r = \frac{N^{n+1}-N^n}{N^n\Delta t} $$
+
+Suppose we set \( t_1=600 \), \( t_2=1200 \),
+\( N^1=140 \) and \( N^2=250 \). We obtain then \( r=0.0013 \). The exact value is \( r = 0.000694 \)
+The following code plot
+
+
+
+
import numpy as np
+
+# Estimate r
+data = np.loadtxt('ecoli.csv', delimiter=',')
+t_e = data[:,0]
+N_e = data[:,1]
+i = 2 # Data point (i,i+1) used to estimate r
+r = (N_e[i+1] - N_e[i])/(N_e[i]*(t_e[i+1] - t_e[i]))
+print 'Estimated r=%.5f' % r
+# Can experiment with r values and see if the model can
+# match the data better
+
+T = 1200 # cell can divide after T sec
+t_max = 5*T # 5 generations in experiment
+t = np.linspace(0, t_max, 1000)
+dt = t[1] - t[0]
+N = np.zeros(t.size)
+
+N[0] = 100
+for n in range(0, len(t)-1, 1):
+ N[n+1] = N[n] + r*dt*N[n]
+
+import matplotlib.pyplot as plt
+plt.plot(t, N, 'r-', t_e, N_e, 'bo')
+plt.xlabel('time [s]'); plt.ylabel('N')
+plt.legend(['model', 'experiment'], loc='upper left')
+plt.show()
+
+
+We can then change the parameter \( r \) in the program and play around to make a better fit. By now we know that this
+'search bythe eye' approach is not the most optimal one.
+
+
Predator-Prey model from ecology
+
The population dynamics of a simple predator-prey system is a
classical example shown in many biology textbooks when ecological
@@ -695,18 +794,7 @@ scientific method:
trying to extract general behaviors and extract eventual laws or patterns
develop mathematical relations for the uncovered regularities/laws and test these by per forming new experiments
-
-
-
-
-
Lots of data about populations of hares and lynx collected from furs in Hudson Bay, Canada, are available. It is known that the populations oscillate. Why? Here we start by @@ -716,20 +804,7 @@ Here we start by
-
-
-
-
- -
Most mammalian predators rely on a variety of prey, which complicates mathematical modeling; however, a few predators have become highly specialized and seek almost exclusively a single prey species. An example of this simplified predator-prey interaction is seen in Canadian northern forests, where the populations of the lynx and the snowshoe hare are intertwined in a life and death struggle.
@@ -764,23 +839,10 @@ One reason that this particular system has been so extensively studied is that t
-
-
-
-
-
import numpy as np
+import numpy as np
from matplotlib import pyplot as plt
# Load in data file
@@ -800,24 +862,9 @@ plt.savefig('Hudson_Bay_data.pdf')
plt.savefig('Hudson_Bay_data.png')
plt.show()
-
-
-
-
-
-
-
-
We see from the plot that there are indeed fluctuations. We would like to create a mathematical model that explains these @@ -839,21 +886,9 @@ climate and other complicating factors. How significant are these?
-
-
-
-
The classical way (in all books) is to present the Lotka-Volterra equations: $$ @@ -872,41 +907,6 @@ Here,
-
-
-
- -
-
-
-
-
The population of hares evolves due to births and deaths exactly as a bacteria population: $$ @@ -923,17 +923,7 @@ So in fraction \( b\Delta t HL \), the lynx eat hares. This loss of hares must be accounted for. Subtracted in the equation for hares: $$ \Delta H = a\Delta t H^n - b \Delta t H^nL^n$$ -
-
-
-
-
We assume that the primary growth for the lynx population depends on sufficient food for raising lynx kittens, which implies an adequate source of nutrients from predation on hares. Thus, the growth of the lynx population does not only depend of how many lynx there are, but on how many hares they can eat. In a time interval \( \Delta t HL \) hares and lynx can meet, and in a @@ -943,26 +933,9 @@ contribute to the growth of lynx, again just a fraction of \( d\Delta t HL \). In addition, lynx die just as in the population dynamics with one isolated animal population, leading to a loss \( -c\Delta t L \). -
-
The accounting of lynx then looks like $$ \Delta L = d\Delta t H^nL^n - c\Delta t L^n$$ -
-
-
-
-
By writing up the definition of \( \Delta H \) and \( \Delta L \), and putting all assumed known terms \( H^n \) and \( L^n \) on the right-hand side, we have @@ -977,42 +950,16 @@ Note:
-
-
-
-
- -
-
-
-
@@ -1067,28 +1014,23 @@ plt.show()
-
-
-
-If we perform a least-square fitting, we can find optimal values for the parameters \( a \), \( b \), \( d \), \( c \). The optimal parameters are \( a=0.4807 \), \( b=0.02482 \), \( d=0.9272 \) and \( c=0.02756 \). These parameters result in a slightly modified initial conditions, namely \( H(0) = 34.91 \) and \( L(0)=3.857 \). With these parameters we are now ready to solve the equations and plot these data together with the experimental values. +We will later perform a least-square fitting. Then we can find optimal +values for the parameters \( a \), \( b \), \( d \), \( c \). In our calculations here +we set \( a=0.4807 \), \( b=0.02482 \), \( d=0.9272 \) and \( c=0.02756 \). These +parameters result in a slightly modified initial conditions, namely +\( H(0) = 34.91 \) and \( L(0)=3.857 \).
-
-
-
+The following Python demonstrates how we can use linear regression to fit for example the population of lynx. +Similarly, we have also used a decision tree algorithm to fit the lynx population data. As expected, the linear regression is not exactly impressive
-
import numpy as np
+import numpy as np
import matplotlib.pyplot as plt
from IPython.display import display
import sklearn
@@ -1107,21 +1049,12 @@ plt.plot(line, regline.predict(line), label= "
plt.plot(x, y, label= "Linear Regression")
plt.show()
-
-
-
-
-
+The similar code for linear regression in R reads
-
HudsonBay = read.csv("src/Hudson_Bay.csv",header=T)
+HudsonBay = read.csv("src/Hudson_Bay.csv",header=T)
fix(HudsonBay)
dim(HudsonBay)
names(HudsonBay)
@@ -1142,188 +1075,7 @@ confint(linearMod)
predict(linearMod,data.frame(Year=c(1910,1914,1920)),interval="confidence")
-
-
-
-
-
-The population grows faster and faster. Why? Is there an underlying (general) mechanism? -
- -
N[n+1] = N[n] + r*dt*N[n]
-
-
-
-
-Let us solve the difference equation in as simple way as possible, -just to train some programming: \( r=1.5 \), \( N^0=1 \), \( \Delta t=0.5 \) - -
- - -
import numpy as np
-
-t = np.linspace(0, 10, 21) # 20 intervals in [0, 10]
-dt = t[1] - t[0]
-N = np.zeros(t.size)
-
-N[0] = 1
-r = 0.5
-
-for n in range(0, N.size-1, 1):
- N[n+1] = N[n] + r*dt*N[n]
- print 'N[%d]=%.1f' % (n+1, N[n+1])
-
-% if FORMAT != 'ipynb':
-
-
-
- - -
N[1]=1.2
-N[2]=1.6
-N[3]=2.0
-N[4]=2.4
-N[5]=3.1
-N[6]=3.8
-N[7]=4.8
-N[8]=6.0
-N[9]=7.5
-N[10]=9.3
-N[11]=11.6
-N[12]=14.6
-N[13]=18.2
-N[14]=22.7
-N[15]=28.4
-N[16]=35.5
-N[17]=44.4
-N[18]=55.5
-N[19]=69.4
-N[20]=86.7
--% endif - -
-
-
-
-
- -
-Use experimental data in the fraction, say \( t_1=600 \), \( t_2=1200 \), -\( N^1=140 \), \( N^2=250 \): \( r=0.0013 \). -
-
-
-
- - -
-
-
- - -
import numpy as np
-
-# Estimate r
-data = np.loadtxt('ecoli.csv', delimiter=',')
-t_e = data[:,0]
-N_e = data[:,1]
-i = 2 # Data point (i,i+1) used to estimate r
-r = (N_e[i+1] - N_e[i])/(N_e[i]*(t_e[i+1] - t_e[i]))
-print 'Estimated r=%.5f' % r
-# Can experiment with r values and see if the model can
-# match the data better
-
-T = 1200 # cell can divide after T sec
-t_max = 5*T # 5 generations in experiment
-t = np.linspace(0, t_max, 1000)
-dt = t[1] - t[0]
-N = np.zeros(t.size)
-
-N[0] = 100
-for n in range(0, len(t)-1, 1):
- N[n+1] = N[n] + r*dt*N[n]
-
-import matplotlib.pyplot as plt
-plt.plot(t, N, 'r-', t_e, N_e, 'bo')
-plt.xlabel('time [s]'); plt.ylabel('N')
-plt.legend(['model', 'experiment'], loc='upper left')
-plt.show()
-
-Change r in the program and play around to make a better fit!
-
-
-
-
The aim here is to simulate financial transactions among financial agents @@ -1397,8 +1149,7 @@ exponentially decreases with \( m' \). We assume that we have \( N=500 \) agents. In each simulation, we need a sufficiently large number of transactions, say \( 10^7 \). Our aim is find the final equilibrium distribution \( w_m \). In order to do that we would need several runs of the above simulations, at least \( 10^3-10^4 \) runs (experiments). -
Our task is to first set up an algorithm which simulates the above transactions with an initial amount \( m_0 \). The challenge here is to figure out a Monte Carlo simulation based on the @@ -1497,40 +1248,11 @@ $$ showing how money is conserved during a transaction. Select values of \( \lambda =0.25,0.5 \) and \( \lambda=0.9 \) and try to extract the corresponding - equilibrium distributions and compare these with the Gibbs distribution. Comment your results. -Extract a parametrization of the above curves, see for example Patriarca and collaborators and see if you can parametrize the high-end tails of the distributions in terms of power laws. Comment your results. + equilibrium distributions and compare these with the Gibbs distribution. We will use this model to +extract a parametrization of the above curves, see for example Patriarca and collaborators. -
-In the studies above the agents were selected randomly, irrespective of whether we allowed for -saving or not during a transaction. What is often observed is that various agents tend to make preferences for for whom to interact with. We will now study the evolution of the distribution of wealth \( w_m \) by assuming that there is a likelihood -$$ -p_{ij} \propto \vert m_i-m_j\vert^{-\alpha}, -$$ +
-Perform the analysis with and without a saving \( \lambda \) on each transaction and comment your results. -We add to the previous probability the possibility that two agents who interact have performed similar transactions earlier. That is, in addition to being financially close, we assume that the likelihood for interacting increases if two agents have interacted earlier. -We add this feature by modifying the previous likelihood to -$$ -p_{ij} \propto \vert m_i-m_j\vert^{-\alpha}\left(c_{ij}+1\right)^{\gamma}, -$$ - -where \( c_{ij} \) represents the number of previous interactions that have taken place between \( i \) and \( j \). The factor \( 1 \) is added in order to ensure that if they have not interacted earlier they can still interact. Perform similar studies as above with \( N=1000 \), \( \alpha=1.0 \) and \( \alpha=2.0 \) using \( \gamma = 0.0, 1.0, 2.0, 3.0 \) and \( 4.0 \). Plot the wealth distributions for these cases and try to extract eventual power law tails with and without a saving \( \lambda \) in each transaction. Comment your results and compare them with figures 5 and 6 of Goswami and Sen. - -
-
-
-
diff --git a/doc/pub/How2ReadData/html/How2ReadData.html b/doc/pub/How2ReadData/html/How2ReadData.html index ded14b424..c401e5a8f 100644 --- a/doc/pub/How2ReadData/html/How2ReadData.html +++ b/doc/pub/How2ReadData/html/How2ReadData.html @@ -75,46 +75,14 @@ div { text-align: justify; text-justify: inter-word; } None, '___sec5'), ('Non-Linear Least squares in R', 2, None, '___sec6'), - ('Predator-Prey model from ecology', 2, None, '___sec7'), - ('Case study from Hudson bay', 2, None, '___sec8'), - ('Hudson bay data', 2, None, '___sec9'), - ('Plotting the data', 2, None, '___sec10'), - ('Hares and lynx in Hudson bay from 1900 to 1920', - 2, + ('Examples', 2, None, '___sec7'), + ('Ecoli lab experiment', 3, None, '___sec8'), + ('Predator-Prey model from ecology', 3, None, '___sec9'), + ('Simulating financial transactions', 3, None, '___sec10'), + ('Particle in one dimension and velocity distribution', + 3, None, - '___sec11'), - ('Why now create a computer model for the hare and lynx ' - 'populations?', - 2, - None, - '___sec12'), - ('The traditional (top-down) approach', 2, None, '___sec13'), - ('Basic mathematics notation', 2, None, '___sec14'), - ('Basic dynamics of the population of hares', - 2, - None, - '___sec15'), - ('Basic dynamics of the population of lynx', 2, None, '___sec16'), - ('Evolution equations', 2, None, '___sec17'), - ('Adapt the model to the Hudson Bay case', 2, None, '___sec18'), - ('The program', 2, None, '___sec19'), - ('The plot', 2, None, '___sec20'), - ('Linear regression in Python', 2, None, '___sec21'), - ('Linear Least squares in R', 2, None, '___sec22'), - ('Example: ecoli lab experiment', 2, None, '___sec23'), - ('The program', 2, None, '___sec24'), - ('The output', 2, None, '___sec25'), - ('Parameter estimation', 2, None, '___sec26'), - ('A program relevant for the biological problem', - 2, - None, - '___sec27'), - ('Simulating financial transactions', 2, None, '___sec28'), - ('Simulation of Transactions', 3, None, '___sec29'), - ('Particle in one dimension an velocity distribution', - 2, - None, - '___sec30')]} + '___sec11')]} end of tocinfo -->
@@ -590,47 +558,39 @@ a linear \( x \)-dependence we study now a cubic polynomial and use the polynomi-
import numpy as np
-import matplotlib.pyplot as plt
+import matplotlib.pyplot as plt
+import numpy as np
+import random
from sklearn.linear_model import Ridge
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import make_pipeline
+from sklearn.linear_model import LinearRegression
-def f(x):
- """ function to approximate by polynomial interpolation"""
- return x*x*x
-
-# generate points used to plot
-x_plot = np.linspace(0, 10, 100)
-
-# generate points and keep a subset of them
-x = np.linspace(0, 10, 100)
-rng = np.random.RandomState(0)
-rng.shuffle(x)
-x = np.sort(x[:20])
-y = f(x)
-# create matrix versions of these arrays
-X = x[:, np.newaxis]
-X_plot = x_plot[:, np.newaxis]
-
-colors = ['teal', 'yellowgreen', 'gold']
-lw = 2
-plt.plot(x_plot, f(x_plot), color='cornflowerblue', linewidth=lw,
- label="ground truth")
-plt.scatter(x, y, color='navy', s=30, marker='o', label="training points")
-
-for count, degree in enumerate([3, 4, 5]):
- model = make_pipeline(PolynomialFeatures(degree), Ridge())
- model.fit(X, y)
- y_plot = model.predict(X_plot)
- plt.plot(x_plot, y_plot, color=colors[count], linewidth=lw,
- label="degree %d" % degree)
-
-plt.legend(loc='lower left')
+x=np.linspace(0.02,0.98,200)
+noise = np.asarray(random.sample((range(200)),200))
+y=x**3*noise
+yn=x**3*100
+poly3 = PolynomialFeatures(degree=3)
+X = poly3.fit_transform(x[:,np.newaxis])
+clf3 = LinearRegression()
+clf3.fit(X,y)
+Xplot=poly3.fit_transform(x[:,np.newaxis])
+poly3_plot=plt.plot(x, clf3.predict(Xplot), label='Cubic Fit')
+plt.plot(x,yn, color='red', label="True Cubic")
+plt.scatter(x, y, label='Data', color='orange', s=15)
+plt.legend()
plt.show()
+
+def error(a):
+ for i in y:
+ err=(y-yn)/yn
+ return abs(np.sum(err))/len(err)
+
+print (error(y))
+Similarly, using R, we can perform similar studies. The following R code illustrates this.
Non-Linear Least squares in R
@@ -664,6 +624,9 @@ text(0, 0.5
+
+In our lectures on regression analysis (and other ones as well), we will discuss in more details various R functionalities.
+
Another useful Python package is
pandas, which is an open source library
@@ -682,11 +645,147 @@ display(data_pandas)
-
Predator-Prey model from ecology
+Examples
+We present here several examples, with pertinent Python codes that we
+will us to illustrate various machine learning methods and ways to
+analyze, from simple to complex, various data sets. Many of these
+examples allow us to generate the data we want to analyze, following
+much of the same philosophy we discussed above when
+fitting various polynomials.
+
+
+We start with a simple exponential growth model that is meant to mimick an ecoli lab experiment.
+We can easily model this system and then produce the data used to train various machine learning algorithms.
+Another model from the life sciences is the so-called predator-prey model from ecology. Thereafter we present
+a simple model for financial transactions before moving to a random walk model and ending with
+the simulation of velocities of a non-interacting atom or molecule confined to move in a one-dimensional region.
+
+
Ecoli lab experiment
+
+
+A typical pattern seen in population models is that the population grows faster and faster. Why? Is there an underlying (general) mechanism?
+Here we will construct a model for cell growth based on a simple difference equation for the growth. We make the following assumptions
+
+
+
+- Cells divide after \( T \) seconds on average (one generation)
+- \( 2N \) celles divide into twice as many new cells \( \Delta N \) in a time
+ interval \( \Delta t \) as \( N \) cells would: \( \Delta N \propto N \)
+- \( N \) cells result in twice as many new individuals \( \Delta N \) in
+ time \( 2\Delta t \) as in time \( \Delta t \): \( \Delta N \propto\Delta t \)
+- Same proportionality wrt death
+- Proposed model: \( \Delta N = b\Delta t N - d\Delta tN \) for some unknown
+ constants \( b \) (births) and \( d \) (deaths)
+- Describe evolution in discrete time: \( t_n=n\Delta t \)
+- Program-friendly notation: \( N \) at \( t_n \) is \( N^n \)
+- Math model: \( N^{n+1} = N^n + r\Delta t\, N \) (with \( \ r=b-d \))
+- Program model:
N[n+1] = N[n] + r*dt*N[n]
+
+
+
+
+
+The difference equation can be programmed in a simple was, and in order to get started we
+set \( r=1.5 \), \( N^0=1 \), \( \Delta t=0.5 \). The program reads
+
+
+
+
+
import numpy as np
+
+t = np.linspace(0, 10, 21) # 20 intervals in [0, 10]
+dt = t[1] - t[0]
+N = np.zeros(t.size)
+
+N[0] = 1
+r = 0.5
+
+for n in range(0, N.size-1, 1):
+ N[n+1] = N[n] + r*dt*N[n]
+ print 'N[%d]=%.1f' % (n+1, N[n+1])
+
+
+and it generates the following output
+
+
+
+
N[1]=1.2
+N[2]=1.6
+N[3]=2.0
+N[4]=2.4
+N[5]=3.1
+N[6]=3.8
+N[7]=4.8
+N[8]=6.0
+N[9]=7.5
+N[10]=9.3
+N[11]=11.6
+N[12]=14.6
+N[13]=18.2
+N[14]=22.7
+N[15]=28.4
+N[16]=35.5
+N[17]=44.4
+N[18]=55.5
+N[19]=69.4
+N[20]=86.7
+
+
+This forms our data which later will define our training set.
+In this case we defined the value of the parameter \( r \). We could alternatively assume that we just received the
+above data file and where asked to use find \( r \). How can we estimate \( r \) from data?
+
+
+We can use the difference equation with the experimental data
+$$ N^{n+1} = N^n + r\Delta t N^n$$
+
+Suppose now that \( N^{n+1} \) and \( N^n \) are known from data. Then we could solve with respect to \( r \) as follows
+$$ r = \frac{N^{n+1}-N^n}{N^n\Delta t} $$
+
+Suppose we set \( t_1=600 \), \( t_2=1200 \),
+\( N^1=140 \) and \( N^2=250 \). We obtain then \( r=0.0013 \). The exact value is \( r = 0.000694 \)
+The following code plot
+
+
+
+
import numpy as np
+
+# Estimate r
+data = np.loadtxt('ecoli.csv', delimiter=',')
+t_e = data[:,0]
+N_e = data[:,1]
+i = 2 # Data point (i,i+1) used to estimate r
+r = (N_e[i+1] - N_e[i])/(N_e[i]*(t_e[i+1] - t_e[i]))
+print 'Estimated r=%.5f' % r
+# Can experiment with r values and see if the model can
+# match the data better
+
+T = 1200 # cell can divide after T sec
+t_max = 5*T # 5 generations in experiment
+t = np.linspace(0, t_max, 1000)
+dt = t[1] - t[0]
+N = np.zeros(t.size)
+
+N[0] = 100
+for n in range(0, len(t)-1, 1):
+ N[n+1] = N[n] + r*dt*N[n]
+
+import matplotlib.pyplot as plt
+plt.plot(t, N, 'r-', t_e, N_e, 'bo')
+plt.xlabel('time [s]'); plt.ylabel('N')
+plt.legend(['model', 'experiment'], loc='upper left')
+plt.show()
+
+
+We can then change the parameter \( r \) in the program and play around to make a better fit. By now we know that this
+'search bythe eye' approach is not the most optimal one.
+
+
Predator-Prey model from ecology
+
The population dynamics of a simple predator-prey system is a
classical example shown in many biology textbooks when ecological
@@ -700,18 +799,7 @@ scientific method:
trying to extract general behaviors and extract eventual laws or patterns
develop mathematical relations for the uncovered regularities/laws and test these by per forming new experiments
-
-
-
-
-
Lots of data about populations of hares and lynx collected from furs in Hudson Bay, Canada, are available. It is known that the populations oscillate. Why? Here we start by @@ -721,20 +809,7 @@ Here we start by
-
-
-
-
- -
Most mammalian predators rely on a variety of prey, which complicates mathematical modeling; however, a few predators have become highly specialized and seek almost exclusively a single prey species. An example of this simplified predator-prey interaction is seen in Canadian northern forests, where the populations of the lynx and the snowshoe hare are intertwined in a life and death struggle.
@@ -769,19 +844,6 @@ One reason that this particular system has been so extensively studied is that t
-
-
-
-
@@ -805,24 +867,9 @@ plt.savefig( plt.savefig('Hudson_Bay_data.png') plt.show()
-
-
-
-
-
-
We see from the plot that there are indeed fluctuations. We would like to create a mathematical model that explains these @@ -844,21 +891,9 @@ climate and other complicating factors. How significant are these?
-
-
-
-
The classical way (in all books) is to present the Lotka-Volterra equations: $$ @@ -877,41 +912,6 @@ Here,
-
-
-
- -
-
-
-
-
The population of hares evolves due to births and deaths exactly as a bacteria population: $$ @@ -928,17 +928,7 @@ So in fraction \( b\Delta t HL \), the lynx eat hares. This loss of hares must be accounted for. Subtracted in the equation for hares: $$ \Delta H = a\Delta t H^n - b \Delta t H^nL^n$$ -
-
-
-
-
We assume that the primary growth for the lynx population depends on sufficient food for raising lynx kittens, which implies an adequate source of nutrients from predation on hares. Thus, the growth of the lynx population does not only depend of how many lynx there are, but on how many hares they can eat. In a time interval \( \Delta t HL \) hares and lynx can meet, and in a @@ -948,26 +938,9 @@ contribute to the growth of lynx, again just a fraction of \( d\Delta t HL \). In addition, lynx die just as in the population dynamics with one isolated animal population, leading to a loss \( -c\Delta t L \). -
-
The accounting of lynx then looks like $$ \Delta L = d\Delta t H^nL^n - c\Delta t L^n$$ -
-
-
-
-
By writing up the definition of \( \Delta H \) and \( \Delta L \), and putting all assumed known terms \( H^n \) and \( L^n \) on the right-hand side, we have @@ -982,42 +955,16 @@ Note:
-
-
-
-
- -
-
-
-
@@ -1072,24 +1019,19 @@ plt.show()
-
-
-
-If we perform a least-square fitting, we can find optimal values for the parameters \( a \), \( b \), \( d \), \( c \). The optimal parameters are \( a=0.4807 \), \( b=0.02482 \), \( d=0.9272 \) and \( c=0.02756 \). These parameters result in a slightly modified initial conditions, namely \( H(0) = 34.91 \) and \( L(0)=3.857 \). With these parameters we are now ready to solve the equations and plot these data together with the experimental values. +We will later perform a least-square fitting. Then we can find optimal +values for the parameters \( a \), \( b \), \( d \), \( c \). In our calculations here +we set \( a=0.4807 \), \( b=0.02482 \), \( d=0.9272 \) and \( c=0.02756 \). These +parameters result in a slightly modified initial conditions, namely +\( H(0) = 34.91 \) and \( L(0)=3.857 \).
-
-
-
+The following Python demonstrates how we can use linear regression to fit for example the population of lynx. +Similarly, we have also used a decision tree algorithm to fit the lynx population data. As expected, the linear regression is not exactly impressive
@@ -1112,17 +1054,8 @@ plt.plot(line, regline.plot(x, y, label= "Linear Regression") plt.show()
-
-
-
+The similar code for linear regression in R reads
@@ -1147,188 +1080,7 @@ confint(linearMod) predict(linearMod,data.frame(Year=c(1910,1914,1920)),interval="confidence")
-
-
-
-
-The population grows faster and faster. Why? Is there an underlying (general) mechanism? -
- -
N[n+1] = N[n] + r*dt*N[n]
-
-
-
-
-Let us solve the difference equation in as simple way as possible, -just to train some programming: \( r=1.5 \), \( N^0=1 \), \( \Delta t=0.5 \) - -
- - -
import numpy as np
-
-t = np.linspace(0, 10, 21) # 20 intervals in [0, 10]
-dt = t[1] - t[0]
-N = np.zeros(t.size)
-
-N[0] = 1
-r = 0.5
-
-for n in range(0, N.size-1, 1):
- N[n+1] = N[n] + r*dt*N[n]
- print 'N[%d]=%.1f' % (n+1, N[n+1])
-
-% if FORMAT != 'ipynb':
-
-
-
- - -
N[1]=1.2
-N[2]=1.6
-N[3]=2.0
-N[4]=2.4
-N[5]=3.1
-N[6]=3.8
-N[7]=4.8
-N[8]=6.0
-N[9]=7.5
-N[10]=9.3
-N[11]=11.6
-N[12]=14.6
-N[13]=18.2
-N[14]=22.7
-N[15]=28.4
-N[16]=35.5
-N[17]=44.4
-N[18]=55.5
-N[19]=69.4
-N[20]=86.7
--% endif - -
-
-
-
-
- -
-Use experimental data in the fraction, say \( t_1=600 \), \( t_2=1200 \), -\( N^1=140 \), \( N^2=250 \): \( r=0.0013 \). -
-
-
-
- - -
-
-
- - -
import numpy as np
-
-# Estimate r
-data = np.loadtxt('ecoli.csv', delimiter=',')
-t_e = data[:,0]
-N_e = data[:,1]
-i = 2 # Data point (i,i+1) used to estimate r
-r = (N_e[i+1] - N_e[i])/(N_e[i]*(t_e[i+1] - t_e[i]))
-print 'Estimated r=%.5f' % r
-# Can experiment with r values and see if the model can
-# match the data better
-
-T = 1200 # cell can divide after T sec
-t_max = 5*T # 5 generations in experiment
-t = np.linspace(0, t_max, 1000)
-dt = t[1] - t[0]
-N = np.zeros(t.size)
-
-N[0] = 100
-for n in range(0, len(t)-1, 1):
- N[n+1] = N[n] + r*dt*N[n]
-
-import matplotlib.pyplot as plt
-plt.plot(t, N, 'r-', t_e, N_e, 'bo')
-plt.xlabel('time [s]'); plt.ylabel('N')
-plt.legend(['model', 'experiment'], loc='upper left')
-plt.show()
-
-Change r in the program and play around to make a better fit!
-
-
-
-
The aim here is to simulate financial transactions among financial agents @@ -1402,8 +1154,7 @@ exponentially decreases with \( m' \). We assume that we have \( N=500 \) agents. In each simulation, we need a sufficiently large number of transactions, say \( 10^7 \). Our aim is find the final equilibrium distribution \( w_m \). In order to do that we would need several runs of the above simulations, at least \( 10^3-10^4 \) runs (experiments). -
Our task is to first set up an algorithm which simulates the above transactions with an initial amount \( m_0 \). The challenge here is to figure out a Monte Carlo simulation based on the @@ -1502,40 +1253,11 @@ $$ showing how money is conserved during a transaction. Select values of \( \lambda =0.25,0.5 \) and \( \lambda=0.9 \) and try to extract the corresponding - equilibrium distributions and compare these with the Gibbs distribution. Comment your results. -Extract a parametrization of the above curves, see for example Patriarca and collaborators and see if you can parametrize the high-end tails of the distributions in terms of power laws. Comment your results. + equilibrium distributions and compare these with the Gibbs distribution. We will use this model to +extract a parametrization of the above curves, see for example Patriarca and collaborators. -
-In the studies above the agents were selected randomly, irrespective of whether we allowed for -saving or not during a transaction. What is often observed is that various agents tend to make preferences for for whom to interact with. We will now study the evolution of the distribution of wealth \( w_m \) by assuming that there is a likelihood -$$ -p_{ij} \propto \vert m_i-m_j\vert^{-\alpha}, -$$ +
-Perform the analysis with and without a saving \( \lambda \) on each transaction and comment your results. -We add to the previous probability the possibility that two agents who interact have performed similar transactions earlier. That is, in addition to being financially close, we assume that the likelihood for interacting increases if two agents have interacted earlier. -We add this feature by modifying the previous likelihood to -$$ -p_{ij} \propto \vert m_i-m_j\vert^{-\alpha}\left(c_{ij}+1\right)^{\gamma}, -$$ - -where \( c_{ij} \) represents the number of previous interactions that have taken place between \( i \) and \( j \). The factor \( 1 \) is added in order to ensure that if they have not interacted earlier they can still interact. Perform similar studies as above with \( N=1000 \), \( \alpha=1.0 \) and \( \alpha=2.0 \) using \( \gamma = 0.0, 1.0, 2.0, 3.0 \) and \( 4.0 \). Plot the wealth distributions for these cases and try to extract eventual power law tails with and without a saving \( \lambda \) in each transaction. Comment your results and compare them with figures 5 and 6 of Goswami and Sen. - -
-
-
-
diff --git a/doc/pub/How2ReadData/ipynb/How2ReadData.ipynb b/doc/pub/How2ReadData/ipynb/How2ReadData.ipynb index 3e36c775f..106fc8862 100644 --- a/doc/pub/How2ReadData/ipynb/How2ReadData.ipynb +++ b/doc/pub/How2ReadData/ipynb/How2ReadData.ipynb @@ -545,51 +545,43 @@ }, "outputs": [], "source": [ - "import numpy as np\n", "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import random\n", "from sklearn.linear_model import Ridge\n", "from sklearn.preprocessing import PolynomialFeatures\n", "from sklearn.pipeline import make_pipeline\n", + "from sklearn.linear_model import LinearRegression\n", "\n", - "def f(x):\n", - " \"\"\" function to approximate by polynomial interpolation\"\"\"\n", - " return x*x*x\n", + "x=np.linspace(0.02,0.98,200)\n", + "noise = np.asarray(random.sample((range(200)),200))\n", + "y=x**3*noise\n", + "yn=x**3*100\n", + "poly3 = PolynomialFeatures(degree=3)\n", + "X = poly3.fit_transform(x[:,np.newaxis])\n", + "clf3 = LinearRegression()\n", + "clf3.fit(X,y)\n", "\n", - "# generate points used to plot \n", - "x_plot = np.linspace(0, 10, 100)\n", + "Xplot=poly3.fit_transform(x[:,np.newaxis])\n", + "poly3_plot=plt.plot(x, clf3.predict(Xplot), label='Cubic Fit')\n", + "plt.plot(x,yn, color='red', label=\"True Cubic\")\n", + "plt.scatter(x, y, label='Data', color='orange', s=15)\n", + "plt.legend()\n", + "plt.show()\n", "\n", - "# generate points and keep a subset of them \n", - "x = np.linspace(0, 10, 100)\n", - "rng = np.random.RandomState(0)\n", - "rng.shuffle(x)\n", - "x = np.sort(x[:20])\n", - "y = f(x)\n", - "# create matrix versions of these arrays \n", - "X = x[:, np.newaxis]\n", - "X_plot = x_plot[:, np.newaxis]\n", + "def error(a):\n", + " for i in y:\n", + " err=(y-yn)/yn\n", + " return abs(np.sum(err))/len(err)\n", "\n", - "colors = ['teal', 'yellowgreen', 'gold']\n", - "lw = 2\n", - "plt.plot(x_plot, f(x_plot), color='cornflowerblue', linewidth=lw,\n", - " label=\"ground truth\")\n", - "plt.scatter(x, y, color='navy', s=30, marker='o', label=\"training points\")\n", - "\n", - "for count, degree in enumerate([3, 4, 5]):\n", - " model = make_pipeline(PolynomialFeatures(degree), Ridge())\n", - " model.fit(X, y)\n", - " y_plot = model.predict(X_plot)\n", - " plt.plot(x_plot, y_plot, color=colors[count], linewidth=lw,\n", - " label=\"degree %d\" % degree)\n", - "\n", - "plt.legend(loc='lower left')\n", - "\n", - "plt.show()" + "print (error(y))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ + "Similarly, using **R**, we can perform similar studies. The following **R** code illustrates this.\n", "## Non-Linear Least squares in R" ] }, @@ -622,6 +614,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ + "In our lectures on regression analysis (and other ones as well), we will discuss in more details various **R** functionalities. \n", + "\n", + "\n", "Another useful Python package is\n", "[pandas](https://pandas.pydata.org/), which is an open source library\n", "providing high-performance, easy-to-use data structures and data\n", @@ -647,7 +642,199 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Predator-Prey model from ecology\n", + "## Examples\n", + "\n", + "We present here several examples, with pertinent Python codes that we\n", + "will us to illustrate various machine learning methods and ways to\n", + "analyze, from simple to complex, various data sets. Many of these\n", + "examples allow us to generate the data we want to analyze, following\n", + "much of the same philosophy we discussed above when\n", + "fitting various polynomials.\n", + "\n", + "We start with a simple exponential growth model that is meant to mimick an ecoli lab experiment.\n", + "We can easily model this system and then produce the data used to train various machine learning algorithms.\n", + "Another model from the life sciences is the so-called predator-prey model from ecology. Thereafter we present \n", + "a simple model for financial transactions before moving to a random walk model and ending with \n", + "the simulation of velocities of a non-interacting atom or molecule confined to move in a one-dimensional region.\n", + "\n", + "\n", + "### Ecoli lab experiment\n", + "\n", + "A typical pattern seen in population models is that the population grows faster and faster. [Why? Is there an underlying (general) mechanism](http://www.zo.utexas.edu/courses/Thoc/PopGrowth.html)?\n", + "Here we will construct a model for cell growth based on a simple difference equation for the growth. We make the following assumptions\n", + "1. Cells divide after $T$ seconds on average (one generation)\n", + "\n", + "2. $2N$ celles divide into twice as many new cells $\\Delta N$ in a time\n", + " interval $\\Delta t$ as $N$ cells would: $\\Delta N \\propto N$\n", + "\n", + "3. $N$ cells result in twice as many new individuals $\\Delta N$ in\n", + " time $2\\Delta t$ as in time $\\Delta t$: $\\Delta N \\propto\\Delta t$\n", + "\n", + "4. Same proportionality wrt death \n", + "\n", + "5. Proposed model: $\\Delta N = b\\Delta t N - d\\Delta tN$ for some unknown\n", + " constants $b$ (births) and $d$ (deaths)\n", + "\n", + "6. Describe evolution in discrete time: $t_n=n\\Delta t$\n", + "\n", + "7. Program-friendly notation: $N$ at $t_n$ is $N^n$\n", + "\n", + "8. Math model: $N^{n+1} = N^n + r\\Delta t\\, N$ (with $\\ r=b-d$)\n", + "\n", + "9. Program model: `N[n+1] = N[n] + r*dt*N[n]`\n", + "\n", + "\n", + "\n", + "The difference equation can be programmed in a simple was, and in order to get started we\n", + "set $r=1.5$, $N^0=1$, $\\Delta t=0.5$. The program reads" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "t = np.linspace(0, 10, 21) # 20 intervals in [0, 10]\n", + "dt = t[1] - t[0]\n", + "N = np.zeros(t.size)\n", + "\n", + "N[0] = 1\n", + "r = 0.5\n", + "\n", + "for n in range(0, N.size-1, 1):\n", + " N[n+1] = N[n] + r*dt*N[n]\n", + " print 'N[%d]=%.1f' % (n+1, N[n+1])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "and it generates the following output" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + " N[1]=1.2\n", + " N[2]=1.6\n", + " N[3]=2.0\n", + " N[4]=2.4\n", + " N[5]=3.1\n", + " N[6]=3.8\n", + " N[7]=4.8\n", + " N[8]=6.0\n", + " N[9]=7.5\n", + " N[10]=9.3\n", + " N[11]=11.6\n", + " N[12]=14.6\n", + " N[13]=18.2\n", + " N[14]=22.7\n", + " N[15]=28.4\n", + " N[16]=35.5\n", + " N[17]=44.4\n", + " N[18]=55.5\n", + " N[19]=69.4\n", + " N[20]=86.7\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This forms our data which later will define our training set. \n", + "In this case we defined the value of the parameter $r$. We could alternatively assume that we just received the \n", + "above data file and where asked to use find $r$. How can we estimate $r$ from data?\n", + "\n", + "We can use the difference equation with the experimental data" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "N^{n+1} = N^n + r\\Delta t N^n\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Suppose now that $N^{n+1}$ and $N^n$ are known from data. Then we could solve with respect to $r$ as follows" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "r = \\frac{N^{n+1}-N^n}{N^n\\Delta t}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Suppose we set $t_1=600$, $t_2=1200$,\n", + "$N^1=140$ and $N^2=250$. We obtain then $r=0.0013$. The exact value is $r = 0.000694$\n", + "The following code plot" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "# Estimate r\n", + "data = np.loadtxt('ecoli.csv', delimiter=',')\n", + "t_e = data[:,0]\n", + "N_e = data[:,1]\n", + "i = 2 # Data point (i,i+1) used to estimate r\n", + "r = (N_e[i+1] - N_e[i])/(N_e[i]*(t_e[i+1] - t_e[i]))\n", + "print 'Estimated r=%.5f' % r\n", + "# Can experiment with r values and see if the model can\n", + "# match the data better\n", + "\n", + "T = 1200 # cell can divide after T sec\n", + "t_max = 5*T # 5 generations in experiment\n", + "t = np.linspace(0, t_max, 1000)\n", + "dt = t[1] - t[0]\n", + "N = np.zeros(t.size)\n", + "\n", + "N[0] = 100\n", + "for n in range(0, len(t)-1, 1):\n", + " N[n+1] = N[n] + r*dt*N[n]\n", + "\n", + "import matplotlib.pyplot as plt\n", + "plt.plot(t, N, 'r-', t_e, N_e, 'bo')\n", + "plt.xlabel('time [s]'); plt.ylabel('N')\n", + "plt.legend(['model', 'experiment'], loc='upper left')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can then change the parameter $r$ in the program and play around to make a better fit. By now we know that this\n", + "'search bythe eye' approach is not the most optimal one. \n", + "\n", + "\n", + "### Predator-Prey model from ecology\n", "\n", "The population dynamics of a simple predator-prey system is a\n", "classical example shown in many biology textbooks when ecological\n", @@ -664,11 +851,6 @@ "\n", " * develop mathematical relations for the uncovered regularities/laws and test these by per forming new experiments\n", "\n", - "\n", - "\n", - "\n", - "## Case study from Hudson bay\n", - "\n", "Lots of data about populations of hares and lynx collected from furs in Hudson Bay, Canada, are available. It is known that the populations oscillate. Why?\n", "Here we start by\n", "\n", @@ -680,12 +862,6 @@ "\n", "4. using the model predict the evolution other predator-pray systems\n", "\n", - "\n", - "\n", - "## Hudson bay data\n", - "\n", - "\n", - "\n", "Most mammalian predators rely on a variety of prey, which complicates mathematical modeling; however, a few predators have become highly specialized and seek almost exclusively a single prey species. An example of this simplified predator-prey interaction is seen in Canadian northern forests, where the populations of the lynx and the snowshoe hare are intertwined in a life and death struggle.\n", "\n", "One reason that this particular system has been so extensively studied is that the Hudson Bay company kept careful records of all furs from the early 1800s into the 1900s. The records for the furs collected by the Hudson Bay company showed distinct oscillations (approximately 12 year periods), suggesting that these species caused almost periodic fluctuations of each other's populations. The table here shows data from 1900 to 1920.\n", @@ -718,18 +894,12 @@ "