Topics

- Additional applications of least squares to simple data sets
- Power laws and scaling
- Exponential and Arrhenius law
- The shape of a nautilus shell

Further reading

- Galileo’s discovery of scaling laws by M. Peterson, 2002.
- Scaling of cities
- Power law for pedestrian movement
- S. Jay Gould’s thoughts on Allometry
- See What price speed? Specific power required for propulsion of vehicles, by G. Gabrielli and T. von Karman, 1950

Linear least squares was developed for solving problems in surveying and astronomy, as we have seen, but it is now used far and wide. One of the reasons linear least squares is such a widely used technique is that the method is also useful for a variety of nonlinear problems. Two of the most common nonlinear modelling applications of least squares are power laws and exponential laws. It can even be used in polar coordinates to fit a model to the shape of a Nautilus’s shell. In the examples that follow, we’ll see just some examples of the ways you can use least squares.

A “power law”, or “scaling law”, is an equation of the form \[y=kx^n,\] where \(k\) is called a constant of proportionality and \(n\) is called the power, or “elasticity”.

The elasticity is frequently used in economics and ecology to describe the relationships between two variables, and can be calculated using the formula \[n = \frac{x}{y} \frac{dy}{dx}.\] The elasticity is dimensionless and has the convenient interpretation of the percent change in \(y\) per percent change in \(x\).

A power law is not a linear equation, but it can be transformed to a linear equation by taking logarithms of both sides, yielding \[\log y = \log k + n \log x.\] If the relationship between \(x\) and \(y\) is close to a power law, then the proportionality constant and the elasticity can be quickly estimated using linear least squares.

In the 4th century BC, the construction of seige engines and artillary was beginning to have a noticable impact on warfare. Ancient writtings tell of the importance played by artillary firing spears and stones in the seiges of Halicarnassus and Tyre by Alexander the Great and the later Roman seige of Syracuse. But seige-engine construction was a demanding engineering task, particularly because of the tools and materials available at the time.

According to Philon of Byzantium’s book on the design of siege engines in the cities of Alexandria and Rhodes in the 3rd century BC, engineers had determined that the dimensions of their catapults should not be in proportion to the weight of the stones being thrown, but in proportion to the cube-root of the weight of the stones. This formula was so important, Philon even presented a sophisticated geometric technique for accurately calculating this cube root. The reasoning is as follows. The weight of a stone is proportional to its volume. The volume \(V\) depends on the radius \(r\) of a sphere according to the formula \[V = \frac{4}{3} \pi r^3.\] So, the elasticity of the weight relative to the radius is \(3\), or alternatively, the radius’s elasticity with respect to weight is \(1/3\)’rd. To double the impact of a spherical projectile, we have to double the weight. Since the radius is proportional to the cube root of weight, the radius of the projectile and the seige engine itself must be increased to \(\sqrt[3]{2} \approx 1.26\) times the original length measurements.

The importance of power laws was probably repeatedly and independently discovered throughout history. But despite the occasional niche application, the general scientific importance of scaling relationships was not widely recognized until they were discussed by Galileo in *Dialogues Concerning Two New Sciences*. Galileo proposed the “square-cube” law to describe how bones sizes scale as animals get bigger. The naive theory is that if you double the mass of a person or animal, and you double mass of the skeleton correspondingly, everything will work out just fine (As suggested in Gulliver’s travels, for instance.) But Galileo argued that it isn’t so – the skeleton’s mass would have to increase more quickly to support this much weight.

One argument goes like this. Assume standard geometric scaling of lengths, areas, volumes, and mass holds: Mass is proportional to volume \(M \propto V\)), volume is proportional to the cube of the length scale (\(V \propto L^3\)), and area is proportional to the square of the length scale (\(A \propto L^2\)), for all systems. The strength of a bone in a body will be proportional to the cross-sectional area (\(A\)) of the bone. Since the bones have to support the full mass of the animal (\(M_{body}\)), \(A \propto M_{body}\). Now, the cross-sectional area of bone is proportional to the square of diameter, and the mass is proportional to the cube of the diameter, so the cross-sectional area must be proportional to the to the skeletal mass \(M_{skel}\) raised to the 2/3rds power (\(A \propto M_{skel}^{2/3}\)). It now follows that \(M_{skel} \propto M_{body}^{3/2}\).

While Galileo’s general argument against isometric skeleton mass scaling with body mass seems sound, it is only a hypothesis. We can test Galileo’s hypothesis by comparing it’s predictions to data. One such data set for a few common mammals was published by Kayser and Heusner in 1964. If we first log-transform this data set as described above, we can fit a line through it using the method of linear least squares.

[ * Data : hide , shown as table , shown as CSV shown as QR * ]

```
# Data extracted from Figure 7 of
# "Scaling in Biology: The consequences of size"
# by Knut Schmidt-Nielsen
# Journal of Experimental Zoology, 194: 287-308, 1975.
# https://dx.doi.org/10.1002/jez.1401940120
#
# Data originally published by Kayser and Heusner, 1964
#
# Data rows are sequentially, Mouse, Rat, Rabbit, Cat, Dog, Man, Elephant
# Columns ...
# Body mass in kilograms, Skeleton mass in kilograms
0.0202, 0.00104
0.198, 0.0118
2.13, 0.176
3.18, 0.302
20.0, 3.76
73.2, 9.23
7721., 1558.
```

[Show code]

```
#!/usr/bin/env python3
from scipy import *
from scipy import linalg
from matplotlib.pyplot import *
data = loadtxt('skeletonMassData.csv',delimiter=',')
# Best fit power law
ld = log(data[:,0])
A = vstack([ld**0, ld**1]).T
b = log(data[:,1])
pB = linalg.solve(A.T.dot(A), A.T.dot(b))
u1 = 10**linspace(-2,4,10)
v1 = exp(pB[0])*u1**pB[1]
# Linear theory fit
a = 1.
pL = sum(log(data[:,1]) - a*log(data[:,0]))/max(data.shape)
u2 = linspace(1e-2,1e4,10)
v2 = exp(pL)*u2**a
# Galileo's 3/2 theory fit
a = 1.5
pG = sum(log(data[:,1]) - a*log(data[:,0]))/max(data.shape)
u3 = linspace(1e-2,1e4,10)
v3 = exp(pG)*u3**a
# Now, do our plotting
loglog(data[:,0], data[:,1],'o',u1,v1,'r-',u2,v2, 'g--', u3, v3, 'b:')
legend(['Data', \
'Best fit ($y = %.3f x^{%.2f}$)'%(exp(pB[0]), pB[1]), \
'Linear theory ($y = %.3f x$)'%exp(pL), \
'Galilean theory ($y = %.3f x^{3/2}$)'%exp(pG)], loc=0, framealpha=0.)
xlabel('Mass of body (kg)', fontsize=18)
ylabel('Mass of skeleton (kg)', fontsize=18)
subplots_adjust(bottom=0.13,top=0.95) # margin control
savefig('skeletonMassFit.pdf')
savefig('skeletonMassFit.png', transparent=True)
```

We find that the best fit power law does not have an elastic of 1.5 as Galileo hypothesized, but instead has an elasticity of 1.12, a value very close to the linear scaling Galileo argued against. Other studies like Prange, Anderson, and Rahn, 1979 and Dumont, 2010 (see Exercises below) support this result. They consistently find elasticities that are slightly larger than 1, but not as large as the 1.5 Galileo predicted. Got any ideas about where Galileo went wrong?

In biology, the study of how size effects things like shape, anatomy and physiology is called “allometry”. Allometry studies have revealed many interesting examples of how biology depends on size. Two more examples (egg size and trotting speed) are given below. These power laws have inspired many attempts to discern how such simple rules emerge from inherently complex biological systems (e.g. West and Brown, 2005). While the data clearly show the existence of scaling relationships, they continue to inspire scientific conversation as to the reasons.

[ * Data : hide , shown as table , shown as CSV shown as QR * ]

```
# mass of adult, mass of egg
# Brody, 1945, partially aggregated
# http://archive.org/download/problemsofrelati033234mbp/problemsofrelati033234mbp.pdf
113380,1700
4536,165.4
3629,94.5
1020,34
283,14
3.6,0.6
```

[Show code]

```
#!/usr/bin/env python2.7
"""
This python script illustrates how to find a scaling
law that matches data using log-transformed least squares
regression.
We are fitting the log-linear relationship
log(v) = k + m log(u)
which corresponds to the scaling law
m k
v(u) = u e
"""
from numpy import *
from pylab import loglog, show, xlabel, ylabel, savefig, text, subplot, plot
# data file at http://www.math.psu.edu/treluga/450/egg.data
egg_data = loadtxt('egg.data', delimiter=',')
print egg_data
# First, we log-transform our data set
log_data = log(egg_data)/log(10) # use log base 10
# Now, we construct A and b for A x = b.
b = log_data[:, 1:2]
A = hstack([ones((6, 1)), log_data[:, 0:1]])
# solve A x = b for x using gaussian elimination with partial pivoting
x = linalg.solve( A.T.dot(A), A.T.dot(b) )
k, m = float(x[0]), float(x[1]) # float calls convert 0-d arrays to scalars
print (k, m)
## Alternatively, we could use numpy's built-in least-squares
## algorithm. If you look close, you'll see that the answer
## is a little different from the one above. This is because
## the algorithm is using a different calculation method based
## on SVD decomposition of the matrix. SVD is more powerful and
## more robust, but also slower for large problems.
#
(x_alt, residues, a_rank, svd_vals) = linalg.lstsq(A, b)
print (float(x_alt[0]), float(x_alt[1]))
####################
#
# Now, we display our results
#
####################
z = 10**linspace(0, 6, 10)
loglog(egg_data[:, 0], egg_data[:, 1], 'o', z, 10**k*z**m,'r-')
xlabel('Adult bird mass ($M$)', fontsize=18)
ylabel('Egg mass ($M_e$)', fontsize=18)
Equation = "$M_e = %.3f M^{%.3f}$"%(10**k, m)
text(300, 0.5, Equation, fontsize=18)
savefig('bird_to_egg_scaling.png')
savefig('bird_to_egg_scaling.pdf')
#show()
```

[ * Data : hide , shown as table , shown as CSV shown as QR * ]

```
# Body mass (kg), Stride frequency (per minute), type of animal
# Trot-gallop transition
# Data extracted from figure 12 page 304 of Schmidt-Nielsen, 1975
# Scaling in Biology: The Consequences of Size
# https://dx.doi.org/10.1002/jez.1401940120
#
0.0114147294105,527.712615249,mice
0.0119388305755,543.270061821,mice
0.0324501676971,456.358342072,mice
0.0405911949976,436.896483228,rat
0.398121423131,278.480604293,dog
1.14914465892,262.759459364,dog
9.57038373087,208.262936379,dog
28.0199093282,158.029501577,dog
30.6447987596,155.750346704,dog
694.018672869,119.912471239,horse
```

Another family of models that can be fit using linear least squares are exponential models of the form \[y = e^{f(x)}.\] Taking the logarithm of both sides, \[\log y = f(x).\] If \(f(x)\) can be written as a weighted sum of terms depending on \(x\), then linear least squares can estimate these weights.

A classic example of this is population growth. Naturalists have long observed that when a small population becomes established in a resource-rich environment, the number of individuals \(N\) often grows exponentially according to \[N = N _ 0 e^{r t}\] where \(N _ 0\) is the initial population size, \(r\) is the growth rate, and \(t\) is the time elapsed since initiation. Log-transformed, \[\ln N = \ln N _ 0 + {r t}.\] For example, from 1700 to 1850, the population of the United States

grew at a steady exponential rate according to census estimates.

[ * Data : hide , shown as table , shown as CSV shown as QR * ]

```
# year, number of people
#
# estimated sizes of the population of the United States
# of America, 1610 - 2010.
2010,308745538
2000,281421906
1990,248709873
1980,226542199
1970,203302031
1960,179323175
1950,151325798
1940,132164569
1930,123202624
1920,106021537
1910,92228496
1900,76212168
1890,62979766
1880,50189209
1870,38558371
1860,31443321
1850,23191876
1840,17063353
1830,12860702
1820,9638453
1810,7239881
1800,5308483
1790,3929214
1780,2780369
1770,2148076
1760,1593625
1750,1170760
1740,905563
1730,629445
1720,466185
1710,331711
1700,250888
1690,210372
1680,151507
1670,111935
1660,75058
1650,50368
1640,26634
1630,4646
1620,2302
1610,350
```

Based on linear least squares, the population grew according to \[\ln N \approx -38.859 + 0.030184 t, \; \text{and} \; N \approx 2^{(t-1287)/22.96}.\] So, the data indicate the population has been doubling every 23 years, and, based on extrapolation, European’s started colonizing North America around 1287. That is about 200 years too early, we know, probably because we haven’t accounted for rapid immigration.

Another exponential law that least squares can fit is Arrhenius’s law. Arrhenius’s law is a model that proposes the rate or frequency \(k\) of an event depends on the absolute temperature \(T\) according to the equation \[k(T) = e^{a-c/T}\] where \(a\) and \(m\) are empirically determined constants. Svante Arrhenius originally proposed this formula for applications in chemistry, but it appears to apply to some biological phenomena as well.

The constants of Arrhenius’s law can be estimated using linear least squares, by taking the logarithm of both sides of the equation. \[\log k = a - \frac{c}{T}\] Thus, any data set obeying Arrhenius’s law should appear as a straight line when the logarithm of the rate is plotted against the reciprocal of the temperature – we call this an “Arrhenius plots”.

Be warned – one particular catch with Arrhenius’s law is that it relates rate to **absolute temperature**. Thus, relative temperature scales like Celsius and Fahrenheit which allow for negative temperatures *will not work*. Instead, we have to measure temperature using a scale like degrees Kelvin, where \(0\) is the absolute minimum temperature.

In 1915, Charles and Aleida Snyder observed and recorded the flashing rates of male fireflies in Maryland. Their data relating temperature to flashing rate is shown below.

[ * Data : hide , shown as table , shown as CSV shown as QR * ]

```
#Temperature (C), flashing rate (per minute)
#
# Snyder and Snyder, 1920 data on
# @article{charles1920flashing,
# title={THE FLASHING INTERVAL OF FIREFLIES—ITS TEMPERATURE COEFFICIENT—AN EXPLANATION OF SYNCHRONOUS FLASHING},
# author={Charles, D and Snyder, Aleida v't H},
# journal={American Journal of Physiology--Legacy Content},
# volume={51},
# number={3},
# pages={536--542},
# year={1920},
# publisher={Am Physiological Soc}
# }
#
19.4,8.1
22.3,9.9
22.6,10.0
23.2,11.1
24.1,11.5
26.0,12.6
26.5,12.1
28.3,15.0
28.8,15.4
```

Just from looking at the data, we can see that the rate of flashing increases as the temperature increases, but perhaps not linearly (though see the problem below). It has been proposed in Laidler (Journal of Chemical Education, 1972) that the relationship between temperature and flashing rate actually obeys the Arrhenius Law.

Using linear least-squares and an Arrhenius plot, we see that the Snyder data are well-explained by the Arrhenius law \[k = e^{9.5 - 2510/T}.\]

[Show code]

```
#!/usr/bin/env python2
# Plot best fit Arrhenius Law to Snyder firefly data
from scipy import *
from scipy import linalg
from matplotlib.pyplot import *
data = loadtxt('data/fireflyflashing.csv',delimiter=',')
def Celsius2Kelvin(x):
return 273.15 + x
A = vstack([data[:,0]*0+1, 1/Celsius2Kelvin(data[:,0])]).T
b = log10(data[:,1])
est_parms = linalg.solve(A.T.dot(A), A.T.dot(b))
print "Arrhenius: ", est_parms
u = A[:,1]
v = est_parms[0] + est_parms[1]*u
plot(u,v,'r-',u,b,'o')
xlabel('Inverse temperature (1/Kelvin)',fontsize=18)
ylabel('$log_{10}$ flashing rate (per minute)',fontsize=18)
savefig('figures/fireflyflashing.png',transparent=True)
```

The chambered nautilus is one of the most charismatic “living fossils” on the planet. At first glance, they appear to be some unusual combination of snail, fish, and jellyfish. But the chambered nautili are actually cephelopods, like squid, cuttle fish, and octopi. Nautili have been around since the Cambrian explosion, 500 million years ago and look much the same today as they did then. The famous documentary series “The Undersea World of Jacque Cousteau” memorably filmed this undersea creature, and in 2015 the “fuzzy nautilus” *Allonautilus scrobiculatus* was finally also filmed.

One of the many mysterious things about a chambered nautilus is the shape of its shell. A nautilus shell coils around itself, and when cut in half, reveals a spiral with the regular partitions that gave rise to the moniker “chambered”. (Unfortunately, the nautilus is again threatened with extinction, in part because of our human attraction to this shell). The spiral itself is so attractive that one almost can not help but ask if it might be revealing some deeper law of nature. One way of getting at this question is asking what equation, if any, describes this spiral shape. And fortunately, the math of polar coordinates is just the thing for drawing spirals.

A spiral is a curve that winds a path around a center in a regular manner. One of the oldest definitions of a spiral comes from Sicily.

“If a straight line one extremity of which remains fixed be made to revolve at a uniform rate in a plane until it returns to the position from which it started and if, at the same time as the straight line is revolving, a point move at a uniform rate along the straight line, starting from the fixed extremity, the point will describe a spiral in the plane.” – Archimedes, On spirals, circa 250 BC

This is what we now call an “Archimedean” or linear spiral. If we let the symbols \(t\) represent time, \(\theta\) be the angle of the line, and \(\rho\) be the distance of a point out from the center, then Archimedes’ definition is equivalent to the parametric polar system \[\theta(t) = a t, \quad \rho(t) = b t \] where \(a\) and \(b\) are the velocities of change of the angle and the radius respectively. By solving the first equation for \(t\) and then substituting into the second, we find that this means the Archimedean spiral has its radius proportional to its angle, specifically \[\rho = \frac{b}{a} \theta .\] Archimedes wrote a treatsie on this spiral and its uses.

A little more thought about polar equations reveals that spirals come in many forms. A spiral might grow quadratically (\(\rho = b \theta^2\)), exponentially (\(\rho = e^{b \theta}\), sometimes called “logarithmic spiral” when written \(\ln \rho = b \theta\).), or according to some other crazy function (e.g. \(\rho = 2 \theta + \sin(6 \theta)\)).

Are any of these spirals is a good description of the spiral we see in a nautilus shell? One crude approach to this is to measure the width of each sequential revolution along a common ray from the center of the nautilus, and look for a pattern in how the width changes. For an Archimedean spiral, each revolution should have the same width. When we measure, what we actually find is that the radius between the first to the second spiral is twice the radius of the first spiral, and the radius between the second and third spiral is about 4 times the radius between the first and the second. So the shape of a chambered nautilus shell a closer match to a logarithmic spiral than Archimede’s spiral. More systematic curve fitting supports this.

[ * Data : hide , shown as table , shown as CSV shown as QR * ]

```
#x-coord,y-coord
-0.528897,-0.0648674
-0.525538,0.0469049
-0.508178,-0.206924
-0.492023,0.140459
-0.47521,-0.294637
-0.422145,0.270158
-0.421131,-0.394498
-0.334304,0.345421
-0.327787,-0.497499
-0.258484,-0.558131
-0.237464,0.399509
-0.187465,-0.0507932
-0.187291,0.00660813
-0.18009,-0.609728
-0.160457,-0.111298
-0.153812,0.088078
-0.143745,0.420374
-0.107665,-0.637137
-0.100245,-0.180966
-0.0660265,0.145215
-0.0634161,0.00925503
-0.045472,-0.0512221
-0.0322192,-0.664555
-0.0198792,0.42
-0.0179807,0.0483927
-0.0157719,-0.220496
-0.00559501,0.148053
0.00582319,-0.0725251
0.0372118,-0.682892
0.0541885,-0.0636078
0.0607875,0.120662
0.062759,-0.226776
0.0646392,0.395575
0.0905606,-0.0244428
0.0938829,0.0752451
0.105794,0.0178073
0.112722,-0.689162
0.150444,-0.202871
0.179186,-0.689363
0.182253,0.325733
0.210976,-0.166801
0.233439,0.268177
0.260766,-0.686588
0.265539,-0.106542
0.275561,0.210648
0.296481,0.135056
0.302085,-0.00997609
0.308337,0.0594914
0.336339,-0.671711
0.40588,-0.653794
0.460324,-0.632811
0.520828,-0.605803
0.57227,-0.578768
0.62978,-0.542688
0.672195,-0.503541
0.720651,-0.464413
0.793432,-0.373998
0.851107,-0.283538
0.905807,-0.177963
0.930213,-0.099487
0.951689,0.0092094
0.9611,0.123985
```

[Show code]

```
import numpy
from scipy import pi, array, cos, sin, arctan2, exp, sqrt, log
from numpy.linalg import norm, solve
from matplotlib.pyplot import *
rcParams.update({'font.size':16})
atan2 = arctan2
# First we load our data points, which were calculated assuming
# the origin has been choosen correctly (it hasn't, there is error).
fname="Nautilus_Spiral-points2.csv"
C = numpy.loadtxt(fname, delimiter=',', skiprows=1)
n = max(C.shape) # n = the number of data points -- useful below
# We need to transform our data points from cartesian to polar
# coordinates if we want to fit a line with least squares to estimate the
# spiral parameters.
rho = sqrt(C[:,1]*C[:,1]+C[:,0]*C[:,0])
theta = atan2(C[:,1], C[:,0])
# It so happens that the data points were not well-sorted
# ahead of time. We would like to sort them, in order from closest
# to the center to farthest from the center, so we have to do a
# little transforming, sorting, than transforming back.
B = [ (rho[i],theta[i]) for i in range(n) ]
B.sort()
B = array(B)
rho, theta = B[:,0], B[:,1]
#for i in B:
# print i
# Because arctan2 automatically unwinds to cannonical radian angles,
# we need to add a few 2*pi shifts to the angles to make sure
# everything will line up. Assuming the data points are sequential allong
# the spiral, we can iterate through the data, and every time
# there is a large-enough jump in the angle (from 2*pi to 0),
# we add back in the appropriate multiple of 2*pi.
maxgap = 2.
shift = 0*theta
for i in range(1,n):
shift[i] = shift[i-1] + (1 if theta[i-1] - theta[i] > maxgap else 0)
theta = theta + 2*pi*shift
# Now, we solve the least-squares fit for the
# logarithmic spiral coefficients
# log(r) = p[0] + p[1]*theta
A = numpy.vstack([numpy.ones(n), theta]).T
b = log(rho)
p = solve( A.T.dot(A), A.T.dot(b))
# We can find quite close agreement
# between p[1] and the golden ratio
# phi according to the following formulas.
# Is this a mystically meaningful relationship?
phi = 0.5+0.5*(5**0.5) # golden ratio
vp1 = exp(pi*p[1]/2) # formula involving p[1]
vphi = phi**(2/pi) # formula involving phi
print p, vp1, vphi, "Error fraction: %.3f"%(abs(vp1/vphi - 1))
# For comparison, here we fit a quadratic spiral
# r = a[0] + a[1]*theta + a[2]*theta**2
# This includes the Archimedian spiral as a special case
A = numpy.vstack([numpy.ones(n), theta, theta*theta]).T
b = rho
p_quad = solve( A.T.dot(A), A.T.dot(b))
print p_quad
# Now, let us plot things
figure(1, figsize=(16,9))
subplot(2, 3, 2)
xlabel(r'x-coordinate')
ylabel(r'y-coordinate')
x, y = exp(p[0]+p[1]*theta)*cos(theta), exp(p[0]+p[1]*theta)*sin(theta)
plot(C[:,0], C[:,1], 'ro', x, y, 'b-')
subplot(2, 3, 1)
title('Logarithmic spiral fit')
xlabel(r'Angle ($\theta$)')
ylabel(r'Radius ($\rho$)')
#plot(theta, log(rho), 'ro', theta, p[0]+p[1]*theta, 'k-')
plot(theta, rho, 'ro', theta, exp(p[0]+p[1]*theta), 'b-', linewidth=2)
try:
legend(['Data', 'Fitted spiral'], loc=2, framealpha=0.)
except:
legend(['Data', 'Fitted spiral'], loc=2)
subplot(2, 3, 3)
res = rho - exp(p[0]+p[1]*theta)
plot(theta, res, 'ko')
ylabel('Residuals')
xlabel(r'Angle ($\theta$)')
title('Error = %f'%sum(abs(res)))
subplot(2, 3, 5)
t = numpy.linspace(min(theta),max(theta), 1000)
r = p_quad[0]+p_quad[1]*t + p_quad[2]*t*t
x, y = r*cos(t), r*sin(t)
xlabel(r'x-coordinate')
ylabel(r'y-coordinate')
plot(C[:,0], C[:,1], 'ro', x, y, 'b-')
subplot(2, 3, 4)
title('Quadratic spiral fit')
xlabel(r'Angle ($\theta$)')
ylabel(r'Radius ($\rho$)')
plot(theta, rho, 'ro', t, r, 'b-', linewidth=2)
try:
legend(['Data', 'Fitted spiral'], loc=2, framealpha=0.)
except:
legend(['Data', 'Fitted spiral'], loc=2)
subplot(2, 3, 6)
ylabel('Residuals')
xlabel(r'Angle ($\theta$)')
res = rho - (p_quad[0]+p_quad[1]*theta + p_quad[2]*theta*theta)
plot(theta, res, 'ko')
title('Error = %f'%sum(abs(res)))
tight_layout()
savefig('Nautilus_fit.png',bbox_inches = 'tight', transparent=True)
#show()
```

Both the quadratic and logarithmic spirals fit the Nautilus shell’s shape reasonably well, but the logarithmic spiral’s fit seems to be about twice as good, and visually seems to be a more satisfying fit, particularly at small angles where the quadratic fit seems to fold over itself. The best fit logarithmic spiral has formula \[\rho = e^{0.172 \theta -3.34}.\]

In the examples above, we’ve seen a number of cases of how the method of linear least squares can be used to help us study problems where the phenomena of interest is traditionally thought of as nonlinear. A word of caution should raised though, to point out that we are abusing the linear least squares method. Proper probabilistic model fitting for each of the problems above should account for these nonlinear effects in the distributions of results and errors. For example, standard linear least squares implicitly assumes all errors are drawn from the same distribution. For power laws, however, it is very common to see situations where the errors scale with size in ways that do not correspond with constant variation after logarithmic transformation. We are fortunate that, for many problems, such nonlinear effects do not fundamentally change the results, and that more complicated statistical models often give essentially the same results as those obtained by a very fast and simple linear least squares analysis.

The argument against isometric skeletal scaling we present above is a clumsy one. For one thing, it assumes isometric scaling of a body implies isometric scaling of the skeleton also. Suppose instead, the cross-sectional area of bones scales proportionally to body mass, while the length of the bones scales proportional to the cube root of body mass. What body-to-skeleton mass relationship would this imply?

Kepler’s harmonic law of planetary motion (his 3rd law) states that the square of the orbital period is proportional to the cube of the characteristic distance of the planet from the sun. Test this by fitting a power law to the data below. Give your answer to 4 decimal places.

[ * Data : hide , shown as table , shown as CSV shown as QR * ]

```
#period (days), distance (AU)
#
# Data supporting Kepler's third law
#
# distance is measured relative to the earth's mean distance,
# which is about 9.26e7 miles and 1.4960e8 meters
# Planets, in order: Mercury,Venus,Earth,Mars,Jupiter,Saturn
# from Timoshenko and Young, page 92
#
87.96,0.3871
224.6,0.7233
356.2,1.000
686.9,1.523
4332,5.200
10759,9.540
```

- The simple pendulum was one of the first methods of measuring the progression of time. Below is a simple data set relating the length of a pendulum to the period of its swing with all other variables held constant. Use linear least squares to estimate a power-law scaling relationship between the length and the period.

[ * Data : hide , shown as table , shown as CSV shown as QR * ]

```
# length (cm), period (seconds)
10.0,0.560
20.0,0.860
30.0,1.050
40.0,1.280
50.0,1.390
60.0,1.520
70.0,1.640
```

- In 1890, Charles Steinmetz, one of the first Americans to work on the mathematical theory of electromagnetism, used least squares to interpret the power-law relationship between the maximum magnetization of an iron wire subject to alternating current and the energy lost due to magnetic hysteresis. Use python to estimate this scaling relationship.

[ * Data : hide , shown as table , shown as CSV shown as QR * ]

```
# magnetomotive force, maximum magnetization, hysteretic loss (ergs per cm cubic)
#
# Ewing's observations of hysteretic energy loss in soft iron wire
# according to Steinmetz, 1890, "Note on the law of hysteresis"
#
# maximum magnetization measured in lines per square centimeter
# hysteretic loss measured in ergs per cm cubic
#
01.50,01974,00410
01.95,03830,01160
02.56,05950,02190
03.01,07180,02940
03.76,08790,03990
04.96,10590,05560
06.62,11480,06160
07.04,11960,06590
26.50,13700,08690
75.20,15560,10190
```

Table 2, Page 45 of J. Huxley’s 1932 book Problems of Relative Growth provides a table of red deer antler statistics relative to body mass. Use python to recalculated a scaling relationship between the deer mass and the antler mass from this data. In your answer, present the data in a new table, your parameterization for the scaling law, your fitted parameter values, and a plot comparing the data to the fitted scaling law.

Theropods were an order of dinosaurs that included Tyranosaurous Rex and where also the ancestors of birds. Below is a data set containing records of length and mass of a variety of species of theropods. Fit a power-law model to these data. Estimate the parameters to 3 decimal places.

[ * Data : hide , shown as table , shown as CSV shown as QR * ]

```
# "Dinosaur Species","Body Length (m)","Body Mass (kg)"
"Staurikosaurus pricei",2.08,19
"Coelophysis bauri",2.68,15.3
"Coelophysis rhodesiensis",2.15,13
"Elaphrosaurus bambergi",6.2,210
"Liliensternus liliensterni",5.15,127
"Dilophosaurus wetherilli",6.03,283
"Ceratosaurus nasicoris",5.69,524
"Eustreptospondylus oxoniensis",4.63,218
"Metriacanthosaurus? sp.",3.8,130
"Allosaurus fragilis",7.4,1010
"Allosaurus atrox",7.9,1320
"Gorgosaurus (Paul's Albertosaurus ) libratus",5.8,700
"Gorgosaurus (Paul's Albertosaurus ) libratus",8.6,2500
"Daspletosaurus (Paul's Tyrannosaurus ) torosus",9,2300
"Tarbosaurus (Paul's Tyrannosaurus ) bataar",5.8,760
"Tarbosaurus (Paul's Tyrannosaurus ) bataar",7.7,2100
"Tyrannosaurus rex",10.6,5700
"Deinonychus (Paul's Velociraptor ) antirrhopus",3.06,45
"Deinonychus (Paul's Velociraptor ) antirrhopus",3.43,73
"Velociraptor mongoliensis",2.07,15
"Ornithomimus edmontonicus",3.3,110
"Ornithomimus brevitertius (=edmontonicus)",3.66,144
"Gallimimus (Paul's Ornithomimus ) bullatus",2.15,27
"Gallimimus (Paul's Ornithomimus ) bullatus",6,440
```

In a 2009 article, Head and other authors use a fossil vertebrae to deduce that about 60 million years ago, there existed a neotropical snake that grew to 13 meters in length and weighted more than 1,000 kilograms (42 feet and more than 2,000 pounds). They called this snake

*Titanoboa cerrejonensis*. This discovery and others like it inspired the PBS Secrets of the Dead episode Graveyard of the giant beasts in 2016, also found on (youtube). Let’s see if we agree with the hype by looking through this scientific paper and check their numbers.What were the two regression lines the team obtained for the relationship between between vertabra width and total body length using the extreme values of 60% and 65% for vertebral position?

Create a plot of width vs total body length using known 60% and 65% position data. Add to this plot the regression lines from part (a). Make sure your axes scales are large enough to show the predicted body length of

*Titanoboa*, given that the discovered vertabra was 12 cm wide.Do you think Head et al’s conclusion about body length are supported by the evidence, are contradicted by the evidence, or that the evidence is inconclusive? Defend your opinion.

[ * Data : hide , shown as table , shown as CSV shown as QR * ]

```
# VW 65% (mm), VW 60% (mm),TBL(mm),SVL(mm)
# These are the data column headings ..
#
# Extant snake data used to estimate Titanoboa's TBL and SVL
# VW = (V)ertabral (W)idth at a fixed precentage done spine
# TBL = (t)otal (b)ody (l)ength
# SVL = (s)nout (v)ent (l)ength
#
# @article{bib:Head2009,
# title={Giant boid snake from the Palaeocene neotropics reveals hotter past equatorial temperatures},
# volume={457},
# ISSN={1476-4687},
# url={http://dx.doi.org/10.1038/nature07671},
# DOI={10.1038/nature07671},
# number={7230},
# journal={Nature},
# publisher={Springer Nature},
# author={Head, Jason J. and Bloch, Jonathan I. and Hastings, Alexander K. and Bourque, Jason R. and Cadena, Edwin A. and Herrera, Fabiany A. and Polly, P. David and Jaramillo, Carlos A.},
# year={2009},
# month={Feb},
# pages={715-717},
# }
#
# Data extracted from Supplement, Table 2
#
6.530,6.380,1535,1423
19.04,19.66,2040,1900
21.42,23.50,2480,2390
12.64,12.92,1606,1355
27.08,27.69,3220,2970
28.01,28.20,3434,3129
10.96,11.22,867,803
11.80,12.20,1450,1216
10.54,10.79,1734,1360
17.97,17.61,2330,1960
11.27,11.54,1380,1195
10.22,10.65,1210,1070
11.58,11.66,1700,1450
10.85,11.20,1770,1385
15.08,15.09,2250,1950
10.66,10.80,1720,1490
17.26,17.10,2470,2110
23.99,24.32,3320,2910
17.24,17.67,2510,2190
23.12,24.33,2690,2310
15.69,16.15,1760,1610
```

- A string was wrapped around a metal bar. A 2 ounce weight was hung at one end of the string, and weights W at the other, so as to counterbalance the weight of 2 oz. The more revolutions of string used, the greater an imbalance the friction between string and bar can sustain. Find a law relating the number of revolutions to the maximum sustained weight.

[ * Data : hide , shown as table , shown as CSV shown as QR * ]

```
# # of Revolutions, weight W (oz.)
#
0.25, 02.875
0.50, 04.000
0.75, 05.706
1.00, 08.901
1.25, 12.437
1.50, 14.700
1.75, 19.062
2.00, 26.50
2.25, 33.75
2.50, 40.00
2.75, 52.00
3.00, 76.00
```

Looking at Snyder’s original firefly data relating temperature to flashing rate, you might suspect that the data can be explained just as well or better with a straight line as with the Arrhenius law. Try this and make a case for which provides a better explanation of the data.

An old reference book on ballooning includes the data below to describe the amount of hydrogen gas a balloon would need to reach various altitudes. Find an equation summarizing this data. (Bonus: What was book where these data come from?) (This result will re-appear in our study of ballistic motion)

[ * Data : hide , shown as table , shown as CSV shown as QR * ]

```
# Altitude (meters), Volume of hydrogen (cubic meters)
12900, 1.25
18400, 10
29500, 640
42300, 80000
49700, 1250000
```

As mentioned above, several groups of researchers have collected data on how skeleton mass actually scales with body mass. Let’s try to replicate their results.

The data of Dumont, 2010 is here. Use the data to produce a corrected version of figure 1, giving the power laws for each fit.

Extract the data from Prange, Anderson, and Rahn, 1979 as a csv file, and reproduce figure 1 of the paper. Do you get the same fit?

- A particularly famous and controversial case has been West, Brown, and Enquist’s explanation of metabolic scaling based on fractal geometry.
- Replicate Kleiber’s 3/4-power scaling law discovery for the relationship between body mass and metabolic rate.
- Read West et al.’s paper and discuss it’s explanation of Kleiber’s law.

Find an equation to describe the scaling of a lightbulb’s brightness as a function of the wattage applied to the bulb, based on the historical data below.

[ * Data : hide , shown as table , shown as CSV shown as QR * ]

```
# Luminoscity (candle power), Power (Watts)
#
# Luminosity characteristic of 40-watt tungsten incandescent lamp.
# Steinmetz, 1911, page 254
# https://archive.org/details/matheengineering00steirich
2,12.25
4,16.33
8,21.35
12,25.60
16,28.91
20,31.64
24,34.55
28,37.29
32,39.26
36,41.47
40,44.14
44,45.42
48,47.05
64,54,.31
96,65.73
128,76.77
192,95.24
256,109.0
291,118.2
320,123.1
382,135.6
460,145.2
```