Nonlinear Regression

Rumford data contains 13 observations from an experiment by Count Rumford (1798) on the rate of cooling.

Time – a numeric vector giving the time since the beginning of the experiment (hr)

Temperature – the temperature (degrees Fahrenheit) of the cannon

\[\frac{dT}{dt}=-\theta(T-T_\infty)\;\;\;\;\;\;T(0)=130F,\;T_\infty=T(\infty)=60F\] \[T(t)=T_\infty+b e^{-\theta t}\;\;\;\;\;\;T(0)=T_\infty+b=130,\;\;T_\infty=60\]

\[T(t)=60+70e^{-\theta t}\]

Linear model

\[Z(t)=-\ln\frac{T(t)-60}{70}=\theta t\] \[Z_i=\theta t_i + \epsilon_i,\;\;\;\;\;\epsilon_i \sim N(0,\sigma^2)\]

## 
## Call:
## lm(formula = -log((Temperature - 60)/70) ~ Time - 1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.04143 -0.00753  0.02197  0.04058  0.04354 
## 
## Coefficients:
##       Estimate Std. Error t value Pr(>|t|)    
## Time 0.0092172  0.0003765   24.48  1.3e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03293 on 12 degrees of freedom
## Multiple R-squared:  0.9804, Adjusted R-squared:  0.9787 
## F-statistic: 599.3 on 1 and 12 DF,  p-value: 1.302e-11

## Warning: the model is without an intercept, hence Multiple R-squared is differnt!

One-parameter nonlinear model

\[T_i=60+70e^{-\theta t_i}+\epsilon_i,\;\;\;\;\;\epsilon_i \sim N(0,\sigma^2)\]

Gauss-Newton: \[\theta_{[j+1]}=\theta_{[j]}- \frac{\sum_{i=1}^n (T_i-60-70e^{-\theta_{[j]}t_i})~ t_i e^{-\theta_{[j]}t_i}}{70\sum_{i=1}^n e^{-2\theta_{[j]}t_i}~t_i^2}\]

\(\theta_{[0]}=.1\)

   par(mfrow=c(1,2),pty="s")
        rumnolin <- nls(Temperature ~ 60. + 70. * exp( - theta * Time), start
         = list(theta = 0.1), trace = T)

## 23628.7 :  0.1
## 2341.244 :  0.0002786525
## 79.42445 :  0.008132221
## 44.18935 :  0.009375149
## 44.1558 :  0.009414963
## 44.15579 :  0.00941549

    print(summary(rumnolin))

## 
## Formula: Temperature ~ 60 + 70 * exp(-theta * Time)
## 
## Parameters:
##        Estimate Std. Error t value Pr(>|t|)    
## theta 0.0094155  0.0004201   22.41 3.67e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.918 on 12 degrees of freedom
## 
## Number of iterations to convergence: 5 
## Achieved convergence tolerance: 4.591e-06

\(\theta_{[0]}=\hat{\theta}_{lin}\)

## 44.97028 :  0.009217199
## 44.15599 :  0.009412413
## 44.15579 :  0.009415457
## 44.15579 :  0.009415496

## 
## Formula: Temperature ~ 60 + 70 * exp(-theta * Time)
## 
## Parameters:
##        Estimate Std. Error t value Pr(>|t|)    
## theta 0.0094155  0.0004201   22.41 3.67e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.918 on 12 degrees of freedom
## 
## Number of iterations to convergence: 3 
## Achieved convergence tolerance: 3.53e-07

Two-parameter nonlinear model

\[T(t)=a+(130-a)e^{-\theta t}\] \[T_i=a+(130-a)e^{-\theta t_i}+\epsilon_i,\;\;\;\;\;\epsilon_i \sim N(0,\sigma^2)\] Gauss-Newton: \(a_{[0]}=60,\;\theta_{[0]}=\hat{\theta}_{nl}\)

   par(mfrow=c(1,2),pty="s")
    rumnolin <- nls(Temperature ~ a + (130. - a) * exp( - theta * Time),
        start = list(a = 60., theta = rumlin.fit$coef), trace = T)

## 44.97028 :  60.000000000  0.009217199
## 43.3613 :  66.6796776  0.0102399
## 41.97549 :  71.93454498  0.01124889
## 40.63664 :  76.16068273  0.01224128
## 40.43069 :  83.08517885  0.01418893
## 39.60845 :  87.89591243  0.01607545
## 38.01841 :  91.38848229  0.01788565
## 35.87664 :  94.01469601  0.01961225
## 34.44551 :  98.07767667  0.02289303
## 30.64914 :  100.60031215   0.02585525
## 26.24409 :  103.93925509   0.03113382
## 16.97053 :  106.95942060   0.03934768
## 2.218421 :  107.85122264   0.04793675
## 1.617296 :  107.51124229   0.04797537
## 1.617296 :  107.51162470   0.04797651

    print(summary(rumnolin))

## 
## Formula: Temperature ~ a + (130 - a) * exp(-theta * Time)
## 
## Parameters:
##        Estimate Std. Error t value Pr(>|t|)    
## a     1.075e+02  6.068e-01  177.19  < 2e-16 ***
## theta 4.798e-02  2.589e-03   18.53 1.21e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3834 on 11 degrees of freedom
## 
## Number of iterations to convergence: 14 
## Achieved convergence tolerance: 5.863e-06

    matplot(Time, cbind(Temperature, fitted(rumnolin)), type = "pl", pch = 
        "*", lty = 1., xlab = "Time", ylab = "Temrerature",col=c("black","blue"))
    plot(predict(rumnolin), residuals(rumnolin), xlab = "Predicted Volume",
        ylab = "Residuals")

    qqnorm(residuals(rumnolin), ylab = "Residuals")

\(H_0: T_\infty=60\; (a=60)\)

\(F=\frac{\Delta RSS/1}{s^2}=\frac{44.2-1.6}{1.6/11}=292.9,\;\;\;p-value \sim 0\)

\(t=\frac{|\hat{a}-60|}{s}=\frac{107.5-60}{0.607}=78.3 \neq \sqrt{292.9},\;\;\; p-value \sim 0\)