[simulations] figure for static nonlinearity

simulations/code/normaldata.m

@@ -5,7 +5,7 @@ randn(2, 4)
 
 % simulate tiger weights:
 mu = 220.0;       % mean and ...
-sigma = 30.0;     % ... standard deviation of the tigers in kg
+sigma = 40.0;     % ... standard deviation of the tigers in kg
 for n = [100, 10000]
     fprintf('\nn=%d:\n', n)
     for i = 1:5

simulations/code/normaldata.out

@@ -17,15 +17,15 @@ ans =
 
 
 n=100:
-  m=218kg, std= 31kg
-  m=227kg, std= 29kg
-  m=224kg, std= 28kg
-  m=214kg, std= 26kg
-  m=219kg, std= 30kg
+  m=218kg, std= 39kg
+  m=223kg, std= 39kg
+  m=217kg, std= 40kg
+  m=216kg, std= 36kg
+  m=220kg, std= 40kg
 
 n=10000:
-  m=220kg, std= 30kg
-  m=220kg, std= 30kg
-  m=220kg, std= 30kg
-  m=220kg, std= 30kg
-  m=220kg, std= 30kg
+  m=220kg, std= 40kg
+  m=220kg, std= 40kg
+  m=221kg, std= 40kg
+  m=220kg, std= 40kg
+  m=220kg, std= 41kg

simulations/lecture/normaldata.py

@@ -9,22 +9,22 @@ if __name__ == "__main__":
     # Generally, males vary in total length from 250 to 390 cm and
     # weigh between 90 and 306 kg
     n = 300
-    mu = 200.0
-    sigma = 50.0
-    rng = np.random.RandomState(22281)
+    mu = 220.0
+    sigma = 40.0
+    rng = np.random.RandomState(32281)
     indices = np.arange(n)
-    data = 50.0*rng.randn(len(indices))+200.0
+    data = sigma*rng.randn(len(indices))+mu
 
-    fig = plt.figure(figsize=cm_size(16.0, 8.0))
+    fig = plt.figure(figsize=cm_size(16.0, 6.0))
     spec = gridspec.GridSpec(nrows=1, ncols=2, width_ratios=[3, 1],
-                             left=0.12, bottom=0.17, right=0.97, top=0.98, wspace=0.08)
+                             left=0.12, bottom=0.23, right=0.97, top=0.96, wspace=0.08)
     ax1 = fig.add_subplot(spec[0, 0])
     show_spines(ax1, 'lb')
     ax1.scatter(indices, data, c=colors['blue'], edgecolor='white', s=50)
     ax1.set_xlabel('index')
     ax1.set_ylabel('Weight / kg')
     ax1.set_xlim(-10, 310)
-    ax1.set_ylim(0, 350)
+    ax1.set_ylim(0, 370)
     ax1.set_yticks(np.arange(0, 351, 100))
     
     ax2 = fig.add_subplot(spec[0, 1])
@@ -32,14 +32,14 @@ if __name__ == "__main__":
     xx = np.arange(0.0, 350.0, 0.5)
     yy = st.norm.pdf(xx, mu, sigma)
     ax2.plot(yy, xx, color=colors['red'])
-    bw = 25.0
-    h, b = np.histogram(data, np.arange(0, 351, bw))
+    bw = 20.0
+    h, b = np.histogram(data, np.arange(0, 401, bw))
     ax2.barh(b[:-1], h/np.sum(h)/(b[1]-b[0]), fc=colors['yellow'], height=0.9*bw, align='edge')
     ax2.set_xlabel('pdf / 1/kg')
-    ax2.set_xlim(0, 0.01)
+    ax2.set_xlim(0, 0.012)
     ax2.set_xticks([0, 0.005, 0.01])
     ax2.set_xticklabels(['0', '0.005', '0.01'])
-    ax2.set_ylim(0, 350)
+    ax2.set_ylim(0, 370)
     ax2.set_yticks(np.arange(0, 351, 100))
     ax2.set_yticklabels([])
 

simulations/lecture/simulations.tex

@@ -49,10 +49,11 @@ mean we just add the desired mean $\mu$ to the random numbers:
   \includegraphics[width=1\textwidth]{normaldata}
   \titlecaption{\label{normaldatafig} Generating normally distributed
     data.}{With the help of a computer the weight of 300 tigers can be
-    measured in no time using the \code{randn()} function (left).  We
-    then even now the population distribution, its mean and standard
-    deviation from which the simulated data values were drawn (red
-    line, right).}
+    measured in no time using the \code{randn()} function (left).  By
+    construction we then even know the population distribution (red
+    line, right), its mean (here 220\,kg) and standard deviation
+    (40\,kg) from which the simulated data values were drawn (yellow
+    histogram).}
 \end{figure}
 
 \begin{exercise}{normaldata.m}{normaldata.out}
@@ -60,7 +61,7 @@ mean we just add the desired mean $\mu$ to the random numbers:
   check its output for some (small) input arguments.  Write a little
   script that generates $n=100$ normally distributed data simulating
   the weight of Bengal tiger males with mean 220\,kg and standard
-  deviation 30\,kg. Check the actual mean and standard deviation of
+  deviation 40\,kg. Check the actual mean and standard deviation of
   the generated data. Do this, let's say, five times using a
   for-loop. Then increase $n$ to 10\,000 and run the code again. It is
   so simple to measure the weight of 10\,000 tigers for getting a
@@ -69,10 +70,9 @@ mean we just add the desired mean $\mu$ to the random numbers:
   values as a function of their index.
 \end{exercise}
 
-\subsection{Uniformly distributed data}
+\subsection{Other probability densities}
 \code{rand()}
-
-\subsection{Other distributions}
+gamma
 
 \subsection{Random integers}
 \code{randi()}
@@ -80,6 +80,17 @@ mean we just add the desired mean $\mu$ to the random numbers:
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{Static nonlinearities}
 
+\begin{figure}[t]
+  \includegraphics[width=1\textwidth]{staticnonlinearity}
+  \titlecaption{\label{staticnonlinearityfig} Generating data
+    fluctuating around a function.}{The open probability of the
+    mechontransducer channel in hair cells of the inner ear is a
+    saturating function of the deflection of hairs (left, red line).
+    Measured data will fluctuate around this function (blue dots).
+    Ideally the residuals (yellow histogram) are normally distributed
+    (right, red line).}
+\end{figure}
+
 Example: mechanotransduciton!
 
 draw (and plot) random functions (in statistics chapter?)

simulations/lecture/staticnonlinearity.py

@@ -0,0 +1,55 @@
+import numpy as np
+import scipy.stats as st
+import matplotlib.pyplot as plt
+import matplotlib.gridspec as gridspec
+from plotstyle import colors, cm_size, show_spines
+
+def boltzmann(x, x0, k):
+    return 1.0/(1.0+np.exp(-k*(x-x0)))
+
+if __name__ == "__main__":
+    n = 70
+    xmin = -18.0
+    xmax = 18.0
+    x0 = 2.0
+    k = 0.3
+    sigma = 0.08
+    rng = np.random.RandomState(38281)
+    x = (xmax-xmin)*rng.rand(n) + xmin
+    y = boltzmann(x, x0, k) + sigma*rng.randn(len(x))
+    xx = np.linspace(xmin, xmax, 200)
+    yy = boltzmann(xx, x0, k)
+
+    fig = plt.figure(figsize=cm_size(16.0, 6.0))
+    spec = gridspec.GridSpec(nrows=1, ncols=2,
+                             left=0.10, bottom=0.23, right=0.97, top=0.96, wspace=0.4)
+    ax1 = fig.add_subplot(spec[0, 0])
+    show_spines(ax1, 'lb')
+    ax1.plot(xx, yy, colors['red'], lw=2)
+    ax1.scatter(x, y, c=colors['blue'], edgecolor='white', s=50)
+    ax1.set_xlabel('Hair deflection / nm')
+    ax1.set_ylabel('Open probability')
+    ax1.set_xlim(-20, 20)
+    ax1.set_ylim(-0.17, 1.17)
+    ax1.set_xticks(np.arange(-20.0, 21.0, 10.0))
+    ax1.set_yticks(np.arange(0, 1.1, 0.2))
+    
+    ax2 = fig.add_subplot(spec[0, 1])
+    show_spines(ax2, 'lb')
+    xg = np.linspace(-1.0, 1.01, 200)
+    yg = st.norm.pdf(xg, 0.0, sigma)
+    ax2.plot(xg, yg, colors['red'], lw=2)
+    bw = 0.05
+    h, b = np.histogram(y-boltzmann(x, x0, k), np.arange(-1.0, 1.01, bw))
+    ax2.bar(b[:-1], h/np.sum(h)/(b[1]-b[0]), fc=colors['yellow'], width=0.9*bw, align='edge')
+    ax2.set_xlabel('residuals')
+    ax2.set_ylabel('pdf')
+    ax2.set_xlim(-0.3, 0.3)
+    #ax2.set_ylim(0, 370)
+    #ax2.set_xticks([0, 0.005, 0.01])
+    #ax2.set_xticklabels(['0', '0.005', '0.01'])
+    #ax2.set_yticks(np.arange(0, 351, 100))
+    #ax2.set_yticklabels([])
+
+    fig.savefig("staticnonlinearity.pdf")
+    plt.close()