R/NIST.R

In onls: Orthogonal Nonlinear Least-Squares Regression

NIST <- function(verbose = FALSE)
{
  oldpar <- par(no.readonly = TRUE) 
  on.exit(par(oldpar))
  par(mfrow = c(2, 2))
  
  ### Chwirut2
  x <- c(0.5, 1, 1.75, 3.75, 5.75, 0.875, 2.25, 3.25, 5.25, 0.75, 
       1.75, 2.75, 4.75, 0.625, 1.25, 2.25, 4.25, 0.5, 3, 0.75, 3, 1.5, 
       6, 3, 6, 1.5, 3, 0.5, 2, 4, 0.75, 2, 5, 0.75, 2.25, 3.75, 5.75, 
       3, 0.75, 2.5, 4, 0.75, 2.5, 4, 0.75, 2.5, 4, 0.5, 6, 3, 0.5, 
       2.75, 0.5, 1.75)
  y <- c(92.9, 57.1, 31.05, 11.5875, 8.025, 63.6, 21.4, 14.25, 8.475, 
       63.8, 26.8, 16.4625, 7.125, 67.3, 41, 21.15, 8.175, 81.5, 13.12, 
       59.9, 14.62, 32.9, 5.44, 12.56, 5.44, 32, 13.95, 75.8, 20, 10.42, 
       59.5, 21.67, 8.55, 62, 20.2, 7.76, 3.75, 11.81, 54.7, 23.7, 11.55, 
       61.3, 17.7, 8.74, 59.2, 16.3, 8.62, 81, 4.87, 14.62, 81.7, 17.17, 
       81.3, 28.9)
  ONLS1 <- onls(y ~ exp(-b1 * x)/(b2 + b3 * x), 
              start = list(b1 = 0.1, b2 = 0.005, b3 = 0.01), verbose = verbose)
  if (verbose) print(ONLS1)
  par(mar = c(2, 2, 2, 1))
  plot(ONLS1, main = "Chwirut2")
 
  ### Gauss1
  x <- 1:250
  y <- c(97.6223, 97.8072, 96.6225, 92.5902, 91.2387, 95.327, 90.3504, 
       89.4624, 91.7252, 89.8692, 86.8808, 85.9436, 87.6069, 86.2584, 
       80.7498, 83.0355, 88.2584, 82.0132, 82.741, 83.3003, 81.2785, 
       81.8551, 80.752, 80.0957, 81.0763, 78.8154, 78.386, 79.9339, 
       79.4847, 79.9594, 76.1069, 78.3983, 81.4306, 82.4887, 81.6546, 
       80.8432, 88.6866, 84.7444, 86.8393, 85.9774, 91.2851, 97.2241, 
       93.5173, 94.1016, 101.9176, 98.4313, 110.4214, 107.6628, 111.7288, 
       116.5115, 120.7609, 123.9553, 124.2437, 130.7996, 133.296, 130.7788, 
       132.0565, 138.6584, 142.9252, 142.7215, 144.1249, 147.4377, 148.2647, 
       152.0519, 147.3863, 149.2074, 148.9537, 144.5876, 148.1226, 148.0144, 
       143.8893, 140.9088, 143.4434, 139.3938, 135.9878, 136.3927, 126.7262,
       124.4487, 122.8647, 113.8557, 113.7037, 106.8407, 107.0034, 102.4629, 
       96.093, 94.5756, 86.9882, 84.9015, 81.1802, 76.4012, 67.092, 72.6716, 
       68.1085, 67.9909, 63.3409, 60.5525, 56.1869, 53.6448, 53.7031, 48.0789, 
       42.2126, 45.6518, 41.6973, 41.2495, 39.2135, 37.717, 36.684, 37.3039, 
       37.4328, 37.4501, 32.6465, 31.8435, 31.3995, 26.6891, 32.2532, 27.6101, 
       33.5865, 28.1071, 30.2643, 28.0165, 29.1102, 23.021, 25.6509, 28.503, 
       25.237, 26.1383, 33.5326, 29.252, 27.0985, 26.53, 25.524, 26.6922, 
       24.5527, 27.7176, 25.203, 25.6148, 25.0689, 27.6393, 24.9485, 25.8681, 
       22.4818, 26.9005, 25.3992, 17.9061, 23.7604, 25.8969, 27.6423, 22.861, 
       26.47, 23.7289, 27.5433, 30.5268, 28.0726, 34.9282, 28.2919, 34.1916, 
       35.4121, 37.0934, 40.9833, 39.5392, 47.8012, 47.4631, 51.0417, 54.5807, 
       57.53, 61.4209, 62.7903, 68.5146, 70.2305, 74.4278, 76.5991, 81.6205, 
       83.4221, 79.1745, 88.5699, 85.6653, 86.555, 90.6591, 84.2729, 85.7222, 
       83.107, 82.1688, 80.4257, 78.1569, 79.7969, 77.8438, 74.5033, 71.5729, 
       65.8803, 65.0139, 60.1958, 59.6673, 52.9548, 53.8779, 44.9127, 41.0991, 
       41.6802, 34.5338, 34.8642, 33.1479, 29.5886, 27.2946, 21.9144, 19.0816, 
       24.9029, 19.8234, 16.7555, 18.2456, 17.2355, 16.3493, 13.7129, 14.7568, 
       13.9717, 12.4287, 14.3552, 7.7033, 10.2344, 11.7832, 13.8777, 4.5357, 
       10.0593, 8.4248, 10.5331, 9.6023, 7.8775, 6.2581, 8.8999, 7.8778, 12.5119, 
       10.6621, 6.0354, 6.7907, 8.7835, 4.6003, 8.4009, 7.2166, 10.0174, 7.3313, 
       6.5279, 2.842, 10.3251, 4.791, 8.3771, 6.2644, 2.7062, 8.3623, 8.9837, 
       3.3626, 1.1827, 4.8754)
  ONLS2 <- onls(y ~ b1 * exp(-b2 * x) + b3 * exp(-(x - b4)^2/b5^2) + b6 * 
                exp(-(x - b7)^2/b8^2), start = list(b1 = 100, b2 = 0.01, 
                b3 = 100, b4 = 50, b5 = 20, b6 = 70, b7 = 200, b8 = 20), verbose = verbose)
  if (verbose) print(ONLS2)
  par(mar = c(2, 2, 2, 1))
  plot(ONLS2, main = "Gauss1")
  
  ### ENSO
  x <- 1:168
  y <- c(12.9, 11.3, 10.6, 11.2, 10.9, 7.5, 7.7, 11.7, 12.9, 14.3, 10.9, 
         13.7, 17.1, 14, 15.3, 8.5, 5.7, 5.5, 7.6, 8.6, 7.3, 7.6, 12.7, 11, 
         12.7, 12.9, 13, 10.9, 10.4, 10.2, 8, 10.9, 13.6, 10.5, 9.2, 12.4, 
         12.7, 13.3, 10.1, 7.8, 4.8, 3, 2.5, 6.3, 9.7, 11.6, 8.6, 12.4, 
         10.5, 13.3, 10.4, 8.1, 3.7, 10.7, 5.1, 10.4, 10.9, 11.7, 11.4, 13.7, 
         14.1, 14, 12.5, 6.3, 9.6, 11.7, 5, 10.8, 12.7, 10.8, 11.8, 12.6, 15.7, 
         12.6, 14.8, 7.8, 7.1, 11.2, 8.1, 6.4, 5.2, 12, 10.2, 12.7, 10.2, 14.7, 
         12.2, 7.1, 5.7, 6.7, 3.9, 8.5, 8.3, 10.8, 16.7, 12.6, 12.5, 12.5, 9.8, 
         7.2, 4.1, 10.6, 10.1, 10.1, 11.9, 13.6, 16.3, 17.6, 15.5, 16, 15.2, 
         11.2, 14.3, 14.5, 8.5, 12, 12.7, 11.3, 14.5, 15.1, 10.4, 11.5, 13.4, 
         7.5, 0.6, 0.3, 5.5, 5, 4.6, 8.2, 9.9, 9.2, 12.5, 10.9, 9.9, 8.9, 7.6, 
         9.5, 8.4, 10.7, 13.6, 13.7, 13.7, 16.5, 16.8, 17.1, 15.4, 9.5, 6.1, 
         10.1, 9.3, 5.3, 11.2, 16.6, 15.6, 12, 11.5, 8.6, 13.8, 8.7, 8.6, 8.6, 
         8.7, 12.8, 13.2, 14, 13.4, 14.8)
  ONLS3 <- onls(y ~ b1 + b2 * cos(2 * pi * x/12) + b3 * sin(2 * pi * x/12)
              + b5 * cos(2 * pi * x/b4) + b6 * sin(2 * pi * x/b4)
              + b8 * cos(2 * pi * x/b7) + b9 * sin(2 * pi * x/b7), 
              start = list(b1 = 10, b2 = 1, b3 = 1, b4 = 10, b5 = 0,
                           b6 = 1, b7 = 10, b8 = 1, b9 = 1), verbose = verbose)
  if (verbose) print(ONLS3)
  par(mar = c(2, 2, 2, 1))
  plot(ONLS3, main = "ENSO")
    
  ## Rat43
  x <- 1:15
  y <- c(16.08, 33.83, 65.80, 97.20, 191.55, 326.20, 386.87, 520.53, 
         590.03, 651.92, 724.93, 699.56, 689.96, 637.56, 717.41)
  ONLS4 <- onls(y ~ b1/(1 + exp(b2 - b3 * x))^(1/b4),   
                start = list(b1 = 1000, b2 = -1, b3 = 0.1, b4 = 0.1),
                extend = c(1, 1), verbose = verbose)
  if (verbose) print(ONLS4)
  plot(ONLS4, main = "Rat43") 
}