Math 753/853 HW5: Least-squares fits to models

Solutions

Problem 1. Polynomial least-squares fit

(a) Fit a cubic polynomial to the following data, and make a plot that shows the data as dots and the cubic fit as a smooth curve.

In [1]:
xdata = [-1.0,  -0.8,  -0.6,  -0.4,  -0.2,  0.0,  0.2,  0.4,  0.6,  0.8,  1.0];
ydata = [3.41, 3.19,  2.57,  2.44,  1.90, 1.66, 1.17, 1.46, 1.07, 1.44, 2.2];
In [2]:
# Set up Vandermonde matrix 
m = length(xdata)
n = 4
X = zeros(m,n)
X[:,1] = 1
for j = 2:n
    X[:,j] = X[:,j-1] .* xdata
end

# Solve least-sqaures system X c = ydata. I'm going to use Julia's 
# backslash operater because I am confident of my understanding of 
# what's going on behind the scenes (least squares via QR decomp)
c = X\ydata
Out[2]:
4-element Array{Float64,1}:
  1.57247
 -1.61084
  1.18473
  0.97028
In [3]:
function polyeval(c,x)
    N = length(c);
    p = c[N]
    for n = N-1:-1:1
        p = p.*x + c[n]
    end
    p
end

using PyPlot
plot(xdata, ydata, "bo", label="datapoints")

x=linspace(-1.1, 1.1)
plot(x, polyeval(c,x), "r-", label="least-squares cubic fit")
xlabel("x")
ylabel("y")
grid("on")

(b) Write out the polynomial in the form $P(x) = c_0 + c_1 x + c_2 x^2 + c_3 x^3$ with the coefficients specified as numeric values with three digits.

In [4]:
c
Out[4]:
4-element Array{Float64,1}:
  1.57247
 -1.61084
  1.18473
  0.97028
In [5]:
P(x) = 1.57 - 1.61 x + 1.18 x^2 + 0.97 x^3
LoadError: syntax: extra token "x" after end of expression
while loading In[5], in expression starting on line 1

Problem 2. Exponential fit

(a) Given the following experimental measurements of alpha-particle emission of a radioactive substance, fit an exponential function $y = c \exp(a t)$ to the data using least squares. Make a plot that shows the data as dots and the exponential fit as a smooth curve. What are the values of $c$ and $a$?

In [6]:
tdata = [0,   4,  8, 12, 16, 20, 24, 28, 32, 36, 40, 44, 48];
ydata = [69, 64, 54, 41, 44, 34, 26, 29, 18, 22, 18, 19, 11];

First we take the log of $y = c e^{at}$ to get

\begin{equation*} \log y = \log c + a t \end{equation*}

The linear system for a bunch of $(t_i,y_i)$ data points is then

\begin{equation*} \left[\begin{array}{cc} 1 & t_1 \\ 1 & t_2 \\ 1 & t_3 \\ \vdots & \vdots \end{array} \right] \left[\begin{array}{c} \log c \\ a \end{array}\right] = \left[\begin{array}{c} \log y_1 \\ \log y_2 \\ \log y_3 \\ \vdots \end{array} \right] \end{equation*}

Call that $T x = \log y$, where $x_1 = \log c$ and $x_2 = a$. Do the calculation in Julia.

In [7]:
m = length(tdata)
T = zeros(m,2)
T[:,1] = 1
T[:,2] = tdata
x = T\log(ydata)
c = exp(x[1])
a = x[2]

plot(tdata, ydata, "bo", label="datapoints")
t = linspace(0, 60)
plot(t, c*exp(a*t), "r-", label="least-squares exponential")
xlabel("hours")
ylabel("alpha particle emission rate")
Out[7]:
PyObject <matplotlib.text.Text object at 0x7fef67d23f10>
In [8]:
c, a
WARNING: Base.writemime is deprecated.
  likely near /home/gibson/.julia/v0.5/IJulia/src/kernel.jl:31
WARNING: Base.writemime is deprecated.
  likely near /home/gibson/.julia/v0.5/IJulia/src/kernel.jl:31
WARNING: Base.writemime is deprecated.
  likely near /home/gibson/.julia/v0.5/IJulia/src/kernel.jl:31
WARNING: Base.writemime is deprecated.
  likely near /home/gibson/.julia/v0.5/IJulia/src/kernel.jl:31
WARNING: Base.writemime is deprecated.
  likely near /home/gibson/.julia/v0.5/IJulia/src/kernel.jl:31
WARNING: Base.writemime is deprecated.
  likely near /home/gibson/.julia/v0.5/IJulia/src/kernel.jl:31
WARNING: Base.writemime is deprecated.
  likely near /home/gibson/.julia/v0.5/IJulia/src/kernel.jl:31
WARNING: Base.writemime is deprecated.
  likely near /home/gibson/.julia/v0.5/IJulia/src/kernel.jl:31
WARNING: Base.writemime is deprecated.
  likely near /home/gibson/.julia/v0.5/IJulia/src/kernel.jl:31
WARNING: Base.writemime is deprecated.
  likely near /home/gibson/.julia/v0.5/IJulia/src/kernel.jl:31
WARNING: Base.writemime is deprecated.
  likely near /home/gibson/.julia/v0.5/IJulia/src/kernel.jl:31
WARNING: Base.writemime is deprecated.
  likely near /home/gibson/.julia/v0.5/IJulia/src/kernel.jl:31
WARNING: Base.writemime is deprecated.
  likely near /home/gibson/.julia/v0.5/IJulia/src/kernel.jl:31
in show at /home/gibson/.julia/v0.5/PyCall/src/PyCall.jl
WARNING: Base.writemime is deprecated.
  likely near /home/gibson/.julia/v0.5/IJulia/src/kernel.jl:31
in show at /home/gibson/.julia/v0.5/PyCall/src/PyCall.jl
WARNING: Base.writemime is deprecated.
  likely near /home/gibson/.julia/v0.5/IJulia/src/kernel.jl:31
in show at /home/gibson/.julia/v0.5/PyCall/src/PyCall.jl
WARNING: Base.writemime is deprecated.
  likely near /home/gibson/.julia/v0.5/IJulia/src/kernel.jl:31
in show at /home/gibson/.julia/v0.5/PyCall/src/PyCall.jl
WARNING: Base.writemime is deprecated.
  likely near /home/gibson/.julia/v0.5/IJulia/src/kernel.jl:31
in show at /home/gibson/.julia/v0.5/PyCall/src/PyCall.jl
WARNING: Base.writemime is deprecated.
  likely near /home/gibson/.julia/v0.5/IJulia/src/kernel.jl:31
in show at /home/gibson/.julia/v0.5/PyCall/src/PyCall.jl
WARNING: Base.writemime is deprecated.
  likely near /home/gibson/.julia/v0.5/IJulia/src/kernel.jl:31
in show at /home/gibson/.julia/v0.5/PyCall/src/PyCall.jl
WARNING: Base.writemime is deprecated.
  likely near /home/gibson/.julia/v0.5/IJulia/src/kernel.jl:31
in show at /home/gibson/.julia/v0.5/PyCall/src/PyCall.jl
WARNING: Base.writemime is deprecated.
  likely near /home/gibson/.julia/v0.5/IJulia/src/kernel.jl:31
in show at /home/gibson/.julia/v0.5/PyCall/src/PyCall.jl
WARNING: Base.writemime is deprecated.
  likely near /home/gibson/.julia/v0.5/IJulia/src/kernel.jl:31
in show at /home/gibson/.julia/v0.5/PyCall/src/PyCall.jl
WARNING: Base.writemime is deprecated.
  likely near /home/gibson/.julia/v0.5/IJulia/src/kernel.jl:31
in show at /home/gibson/.julia/v0.5/PyCall/src/PyCall.jl
WARNING: Base.writemime is deprecated.
  likely near /home/gibson/.julia/v0.5/IJulia/src/kernel.jl:31
in show at /home/gibson/.julia/v0.5/PyCall/src/PyCall.jl
Out[8]:
(69.23564314813459,-0.03475031092464803)
In [9]:
# So the exponential fit is 
f(t) = 69.2 * exp(-0.348*t)
Out[9]:
f (generic function with 1 method)

(b) What is the substance's half-life? (i.e. the time $t$ for which $y(t)/y(0) = 1/2$)

In [10]:
# for  y(t) = c exp(at), y(0) = c, so we have
#      y(t) = y(0) exp(at)
# y(t)/y(0) = exp(at) 

# Thus if y(t)/y(0) = 1/2 = exp(at),
# log 1/2 = at, and 
# the halflife is given by t = 1/a log 1/2
t = 1/a * log(1/2)
Out[10]:
19.946502984187784
In [11]:
# Looks like the half-life of the substance is about 20 hrs.

Problem 3. Power-law fit

(a) Fit a power-law curve $y = c t^a$ to the following data, and make a plot showing the datapoints as dots and the fit as a smooth curve.

In [12]:
tdata =  [   2,    3,    4,    5,    6,    7,    8,    9,   10]; 
ydata =  [12.7, 11.2, 8.99, 8.62, 8.12, 8.47, 7.39, 7.24, 6.99];

Take the log of $y = c t^a$ to get

\begin{equation*} \log y = \log c + a \log t \end{equation*}

The linear system for a bunch of $(t_i,y_i)$ data points is then

\begin{equation*} \left[\begin{array}{cc} 1 & \log t_1 \\ 1 & \log t_2 \\ 1 & \log t_3 \\ \vdots & \vdots \end{array} \right] \left[\begin{array}{c} \log c \\ a \end{array}\right] = \left[\begin{array}{c} \log y_1 \\ \log y_2 \\ \log y_3 \\ \vdots \end{array} \right] \end{equation*}

Call that $T x = \log y$, where $x_1 = \log c$ and $x_2 = a$.

Now do the calculation in Julia.

In [13]:
m = length(tdata)
T = zeros(m,2)
T[:,1] = 1
T[:,2] = log(tdata)
x = T\log(ydata)
c = exp(x[1])
a = x[2]

plot(tdata, ydata, "bo", label="datapoints")
t = linspace(1.5, 12)
plot(t, c*t.^a, "r-", label="least-squares power law")
xlabel("t")
ylabel("y")
grid("on")

(b) Write out the least-squares power-law fit $y = c t^a$ with $c$ and $a$ specified as numeric values with three digits.

In [14]:
c, a
Out[14]:
(16.04696670297259,-0.3652492064261485)
In [15]:
y(t) = 16.05 * t^-0.365
Out[15]:
y (generic function with 1 method)

Problem 4. c t exp(at) fit

(a) The following data represent measurements of blood concentration of a drug after intravenous injection as a function of time. Fit a function of the form $y = c \, t \, e^{at}$ to the data using least squares. Make a plot that shows the data as dots and the exponential fit as a smooth curve.

In [16]:
# t is time in hours, y is concentration in ng/ml
tdata = [4,   8,  12,  16,  20,  24];
ydata = [21, 31,  25,  21,  15,  16];

Rewrite $y = c t e^{at}$ as $y/t = c e^{at}$ and then take the log of both sides to get

\begin{equation*} \log y/t = \log c + a t \end{equation*}

The linear system for a bunch of $(t_i,y_i)$ data points is then

\begin{equation*} \left[\begin{array}{cc} 1 & t_1 \\ 1 & t_2 \\ 1 & t_3 \\ \vdots & \vdots \end{array} \right] \left[\begin{array}{c} \log c \\ a \end{array}\right] = \left[\begin{array}{c} \log y_1/t_1 \\ \log y_2/t_2 \\ \log y_3/t_3 \\ \vdots \end{array} \right] \end{equation*}

Call that $A x = \log y/t$, where $x_1 = \log c$ and $x_2 = a$.

Now do the calculation in Julia.

In [17]:
m = length(tdata)
A = zeros(m,2)
A[:,1] = 1
A[:,2] = tdata
x = A\log(ydata./tdata)
c = exp(x[1])
a = x[2]

plot(tdata, ydata, "bo", label="datapoints")
t = linspace(0, 48)
plot(t, c*t.*exp(a*t), "r-", label="least-squares ct exp(at) fit")
xlabel("t")
ylabel("y")
grid("on")

(b) Write out the model $y = c \, t \, e^{at}$ with numeric values specified to three digits.

In [18]:
c, a
Out[18]:
(8.372674496790392,-0.11219417563518795)
In [19]:
y(t) = 8.37*t * exp(-0.112*t)
WARNING: Method definition y(Any) in module Main at In[15]:1 overwritten at In[19]:1.
Out[19]:
y (generic function with 1 method)

(c) Based on the model, at what time do you expect the concentration to reach 5 ng/ml?

In [20]:
# Looks like about 37 hrs from the graph. Could try to solve the nonlinear equation for t,
# i.e. 5 = 8.37 t exp(-0.112 t), but that would be overkill, given the approximate nature
# of the least-squares fit.