Nonlinear Regression

In the last lecture, we started discussing nonlinear regression models. As a motivating example, we considered the following dataset (from FRED) on Annual Estimates of the Resident Population of California (units are thousands of persons) from 1900 to 2024.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm

capop = pd.read_csv('CAPOP_11Sept2025.csv')
print(capop.head(10))
print(capop.tail(10))
tme = np.arange(1900, 2025)
plt.plot(tme, capop['CAPOP'], label='California Population')
plt.xlabel('Year')
plt.ylabel('Population (in thousands)')
plt.title('California Annual Resident Population (in thousands)')
plt.show()

We work with the logarithms of the population data as this will lead to models with better interpretability.

y = np.log(capop['CAPOP'])
n = len(y)
plt.plot(tme, y)
plt.xlabel('Year')
plt.ylabel('Log(Population in thousands)')
plt.title('Logarithm of California Annual Resident Population (in thousands)')
plt.show()

If we fit a linear regression model to this data, we will get an estimate of the overall growth rate.

x = np.arange(1, n+1)
X = np.column_stack((np.ones(n), x))
linreg = sm.OLS(y, X).fit()
print(linreg.summary())
plt.plot(tme, y, label='Log(Population in thousands)')
plt.plot(tme, linreg.fittedvalues, color='red', label='Fitted least squares line')
plt.xlabel('Year')
plt.ylabel('Log(Population in thousands)')
plt.title('Log Data with fitted least squares line')
plt.show()

The fitted slope coefficient here is 0.0269. The interpretation is that the population increases by 2.69% each year.

The simple linear regression model does not provide a good fit to the data. The quality of the fit can be assessed using the Residual Sum of Squares (RSS).

rss_lm = np.sum(linreg.resid ** 2)
print(rss_lm)

6.188666219005153

The simple linear regression model is not a good fit to the data. It is clear that the population growth rate is not 2.69% uniformly. In the initial years, the growth rate seems to higher than 2.69%, and in recent years, it seems to be lower. The simple linear regression model cannot pick up these variable growth rates. We can instead consider the following model:

y_t = \beta_0 + \beta_1 t + \beta_2 (t - c)_+ + \epsilon_t

This model uses $\beta_1$ for the slope before $c$ , and $\beta_1 + \beta_2$ for the slope after $c$ .

If $c$ is known, then this is again linear regression (but now it is multiple linear regression as opposed to simple linear regression), and we can fit this model as follows.

c = 25 #some arbitrary value
x = np.arange(1, n+1)
x_c = np.maximum(0, x - c)
X = np.column_stack((np.ones(n), x, x_c))
linreg_c_fixed = sm.OLS(y, X).fit()
print(linreg_c_fixed.summary())
plt.plot(tme, y, label='Log(Population in thousands)')
plt.plot(tme, linreg_c_fixed.fittedvalues, color='red', label='Fitted least squares line')
plt.xlabel('Year')
plt.ylabel('Log(Population in thousands)')
plt.title('Log Data with fitted least squares line (c fixed)')
plt.show()

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                  CAPOP   R-squared:                       0.976
Model:                            OLS   Adj. R-squared:                  0.976
Method:                 Least Squares   F-statistic:                     2487.
Date:                Tue, 15 Sep 2026   Prob (F-statistic):           1.33e-99
Time:                        19:05:55   Log-Likelihood:                 56.364
No. Observations:                 125   AIC:                            -106.7
Df Residuals:                     122   BIC:                            -98.24
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          7.1511      0.058    124.023      0.000       7.037       7.265
x1             0.0589      0.003     20.971      0.000       0.053       0.064
x2            -0.0355      0.003    -11.498      0.000      -0.042      -0.029
==============================================================================
Omnibus:                       12.276   Durbin-Watson:                   0.011
Prob(Omnibus):                  0.002   Jarque-Bera (JB):                5.611
Skew:                          -0.286   Prob(JB):                       0.0605
Kurtosis:                       2.134   Cond. No.                         368.
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

To assess the quality of fit of this model, we again look at the residual sum of squares.

rss_lm_c_fixed = np.sum(linreg_c_fixed.resid ** 2)
print(rss_lm_c_fixed, rss_lm)

2.97016881050577 6.188666219005153

The RSS now is much smaller than the RSS for the simple linear regression model. If we change the value of $c$ , there is the possibility of getting even smaller RSS.

c = 75 #some arbitrary value
x = np.arange(1, n+1)
x_c = np.maximum(0, x - c)
X = np.column_stack((np.ones(n), x, x_c))
linreg_c_fixed = sm.OLS(y, X).fit()
rss_lm_c_fixed = np.sum(linreg_c_fixed.resid ** 2)
print(rss_lm_c_fixed)

0.5429018431026454

A natural way of estimating $c$ is the following:

For each fixed value of $c$ , calculate RSS
Use the value of $c$ with the smallest RSS as the estimate

The following function calculates the value of RSS for each fixed value of $c$ .

def rss(c):
    x = np.arange(1, n+1)
    x_c = np.maximum(0, x - c)
    X = np.column_stack((np.ones(n), x, x_c))
    md = sm.OLS(y, X).fit()
    rss = np.sum(md.resid ** 2)
    return rss

We compute $RSS(c)$ for each $c \in \{1, \dots, n\}$ as follows.

num_c_vals = 1000 #this is the number of different values of c we will try
allcvals = np.linspace(1, n, num_c_vals)
rssvals = np.array([rss(c) for c in allcvals])
plt.plot(allcvals, rssvals)
plt.show()

The estimate $\hat{c}$ is obtained by minimizing $RSS(c)$ as follows.

c_hat = allcvals[np.argmin(rssvals)]
print(c_hat)
print(c_hat - 1 + tme[0]) #this is the estimated year when the slope changes
rss_smallest = np.min(rssvals)
print(rss_smallest, rss_lm, rss_lm_c_fixed)

66.28928928928929
1965.2892892892892
0.34904607760661777 6.188666219005153 0.5429018431026454

The fitted values will now look much better than before.

c = c_hat
x = np.arange(1, n+1)
x_c = np.maximum(0, x - c)
X = np.column_stack((np.ones(n), x, x_c))
md = sm.OLS(y, X).fit()
plt.plot(tme, y)
plt.plot(tme, md.fittedvalues, color = 'red')
plt.axvline(c_hat - 1 + tme[0], color='green', linestyle='--')
plt.show()

Point estimates of the other parameters are obtained as follows.

#Estimates of other parameters: 
c = c_hat
x = np.arange(1, n+1)
x_c = np.maximum(0, x - c)
X = np.column_stack((np.ones(n), x, x_c))
md = sm.OLS(y, X).fit()
print(md.params) #this gives estimates of beta_0, beta_1, beta_2

const    7.369373
x1       0.037996
x2      -0.024062
dtype: float64

The estimate of $\hat{\beta}_1$ is 0.038 and the estimate of $\hat{\beta}_2$ is -0.024. This means that the growth rate before 1965 was 3.8% while the growth rate after 1965 is $3.8 - 2.4 = 1.4$ %.

Bayesian Posterior¶

Next step is uncertainty quantification. For uncertainty quantification, we do Bayesian analysis. The posterior density for $c$ is given by:

\begin{align*} \text{posterior}(c) \propto I\{1 < c < n\} |X_c^T X_c|^{-1/2} \left(\frac{1}{RSS(c)} \right)^{(n-p)/2} \end{align*}

(1)

Note that the formula on the right hand side gives the unnormalized posterior (because the proportionality sign hides the normalizing constant). The function below computes this posterior on the log-scale for fixed $c$ . Computing the logarithm is numerically much more stable.

#Bayesian log posterior
def logpost(c):
    x = np.arange(1, n+1)
    x_c = np.maximum(0, x - c)
    X = np.column_stack([np.ones(n), x, x_c])
    p = X.shape[1]
    md = sm.OLS(y, X).fit()
    rss = np.sum(md.resid ** 2)
    sgn, log_det = np.linalg.slogdet(np.dot(X.T, X)) #sgn gives the sign of the determinant (in our case, this should 1)
    #log_det gives the logarithm of the absolute value of the determinant
    logval = ((p-n)/2) * np.log(rss) - (0.5)*log_det
    return logval

We now evaluate the log posterior on a grid of values of $c$ . We will not consider values of $c$ that are too close to 1 or $n$ to avoid numerical issues arising from near singularity of the determinant term.

print(allcvals)

[  1.           1.12412412   1.24824825   1.37237237   1.4964965
   1.62062062   1.74474474   1.86886887   1.99299299   2.11711712
   2.24124124   2.36536537   2.48948949   2.61361361   2.73773774
   2.86186186   2.98598599   3.11011011   3.23423423   3.35835836
   3.48248248   3.60660661   3.73073073   3.85485485   3.97897898
   4.1031031    4.22722723   4.35135135   4.47547548   4.5995996
   4.72372372   4.84784785   4.97197197   5.0960961    5.22022022
   5.34434434   5.46846847   5.59259259   5.71671672   5.84084084
   5.96496496   6.08908909   6.21321321   6.33733734   6.46146146
   6.58558559   6.70970971   6.83383383   6.95795796   7.08208208
   7.20620621   7.33033033   7.45445445   7.57857858   7.7027027
   7.82682683   7.95095095   8.07507508   8.1991992    8.32332332
   8.44744745   8.57157157   8.6956957    8.81981982   8.94394394
   9.06806807   9.19219219   9.31631632   9.44044044   9.56456456
   9.68868869   9.81281281   9.93693694  10.06106106  10.18518519
  10.30930931  10.43343343  10.55755756  10.68168168  10.80580581
  10.92992993  11.05405405  11.17817818  11.3023023   11.42642643
  11.55055055  11.67467467  11.7987988   11.92292292  12.04704705
  12.17117117  12.2952953   12.41941942  12.54354354  12.66766767
  12.79179179  12.91591592  13.04004004  13.16416416  13.28828829
  13.41241241  13.53653654  13.66066066  13.78478478  13.90890891
  14.03303303  14.15715716  14.28128128  14.40540541  14.52952953
  14.65365365  14.77777778  14.9019019   15.02602603  15.15015015
  15.27427427  15.3983984   15.52252252  15.64664665  15.77077077
  15.89489489  16.01901902  16.14314314  16.26726727  16.39139139
  16.51551552  16.63963964  16.76376376  16.88788789  17.01201201
  17.13613614  17.26026026  17.38438438  17.50850851  17.63263263
  17.75675676  17.88088088  18.00500501  18.12912913  18.25325325
  18.37737738  18.5015015   18.62562563  18.74974975  18.87387387
  18.997998    19.12212212  19.24624625  19.37037037  19.49449449
  19.61861862  19.74274274  19.86686687  19.99099099  20.11511512
  20.23923924  20.36336336  20.48748749  20.61161161  20.73573574
  20.85985986  20.98398398  21.10810811  21.23223223  21.35635636
  21.48048048  21.6046046   21.72872873  21.85285285  21.97697698
  22.1011011   22.22522523  22.34934935  22.47347347  22.5975976
  22.72172172  22.84584585  22.96996997  23.09409409  23.21821822
  23.34234234  23.46646647  23.59059059  23.71471471  23.83883884
  23.96296296  24.08708709  24.21121121  24.33533534  24.45945946
  24.58358358  24.70770771  24.83183183  24.95595596  25.08008008
  25.2042042   25.32832833  25.45245245  25.57657658  25.7007007
  25.82482482  25.94894895  26.07307307  26.1971972   26.32132132
  26.44544545  26.56956957  26.69369369  26.81781782  26.94194194
  27.06606607  27.19019019  27.31431431  27.43843844  27.56256256
  27.68668669  27.81081081  27.93493493  28.05905906  28.18318318
  28.30730731  28.43143143  28.55555556  28.67967968  28.8038038
  28.92792793  29.05205205  29.17617618  29.3003003   29.42442442
  29.54854855  29.67267267  29.7967968   29.92092092  30.04504505
  30.16916917  30.29329329  30.41741742  30.54154154  30.66566567
  30.78978979  30.91391391  31.03803804  31.16216216  31.28628629
  31.41041041  31.53453453  31.65865866  31.78278278  31.90690691
  32.03103103  32.15515516  32.27927928  32.4034034   32.52752753
  32.65165165  32.77577578  32.8998999   33.02402402  33.14814815
  33.27227227  33.3963964   33.52052052  33.64464464  33.76876877
  33.89289289  34.01701702  34.14114114  34.26526527  34.38938939
  34.51351351  34.63763764  34.76176176  34.88588589  35.01001001
  35.13413413  35.25825826  35.38238238  35.50650651  35.63063063
  35.75475475  35.87887888  36.003003    36.12712713  36.25125125
  36.37537538  36.4994995   36.62362362  36.74774775  36.87187187
  36.995996    37.12012012  37.24424424  37.36836837  37.49249249
  37.61661662  37.74074074  37.86486486  37.98898899  38.11311311
  38.23723724  38.36136136  38.48548549  38.60960961  38.73373373
  38.85785786  38.98198198  39.10610611  39.23023023  39.35435435
  39.47847848  39.6026026   39.72672673  39.85085085  39.97497497
  40.0990991   40.22322322  40.34734735  40.47147147  40.5955956
  40.71971972  40.84384384  40.96796797  41.09209209  41.21621622
  41.34034034  41.46446446  41.58858859  41.71271271  41.83683684
  41.96096096  42.08508509  42.20920921  42.33333333  42.45745746
  42.58158158  42.70570571  42.82982983  42.95395395  43.07807808
  43.2022022   43.32632633  43.45045045  43.57457457  43.6986987
  43.82282282  43.94694695  44.07107107  44.1951952   44.31931932
  44.44344344  44.56756757  44.69169169  44.81581582  44.93993994
  45.06406406  45.18818819  45.31231231  45.43643644  45.56056056
  45.68468468  45.80880881  45.93293293  46.05705706  46.18118118
  46.30530531  46.42942943  46.55355355  46.67767768  46.8018018
  46.92592593  47.05005005  47.17417417  47.2982983   47.42242242
  47.54654655  47.67067067  47.79479479  47.91891892  48.04304304
  48.16716717  48.29129129  48.41541542  48.53953954  48.66366366
  48.78778779  48.91191191  49.03603604  49.16016016  49.28428428
  49.40840841  49.53253253  49.65665666  49.78078078  49.9049049
  50.02902903  50.15315315  50.27727728  50.4014014   50.52552553
  50.64964965  50.77377377  50.8978979   51.02202202  51.14614615
  51.27027027  51.39439439  51.51851852  51.64264264  51.76676677
  51.89089089  52.01501502  52.13913914  52.26326326  52.38738739
  52.51151151  52.63563564  52.75975976  52.88388388  53.00800801
  53.13213213  53.25625626  53.38038038  53.5045045   53.62862863
  53.75275275  53.87687688  54.001001    54.12512513  54.24924925
  54.37337337  54.4974975   54.62162162  54.74574575  54.86986987
  54.99399399  55.11811812  55.24224224  55.36636637  55.49049049
  55.61461461  55.73873874  55.86286286  55.98698699  56.11111111
  56.23523524  56.35935936  56.48348348  56.60760761  56.73173173
  56.85585586  56.97997998  57.1041041   57.22822823  57.35235235
  57.47647648  57.6006006   57.72472472  57.84884885  57.97297297
  58.0970971   58.22122122  58.34534535  58.46946947  58.59359359
  58.71771772  58.84184184  58.96596597  59.09009009  59.21421421
  59.33833834  59.46246246  59.58658659  59.71071071  59.83483483
  59.95895896  60.08308308  60.20720721  60.33133133  60.45545546
  60.57957958  60.7037037   60.82782783  60.95195195  61.07607608
  61.2002002   61.32432432  61.44844845  61.57257257  61.6966967
  61.82082082  61.94494494  62.06906907  62.19319319  62.31731732
  62.44144144  62.56556557  62.68968969  62.81381381  62.93793794
  63.06206206  63.18618619  63.31031031  63.43443443  63.55855856
  63.68268268  63.80680681  63.93093093  64.05505506  64.17917918
  64.3033033   64.42742743  64.55155155  64.67567568  64.7997998
  64.92392392  65.04804805  65.17217217  65.2962963   65.42042042
  65.54454454  65.66866867  65.79279279  65.91691692  66.04104104
  66.16516517  66.28928929  66.41341341  66.53753754  66.66166166
  66.78578579  66.90990991  67.03403403  67.15815816  67.28228228
  67.40640641  67.53053053  67.65465465  67.77877878  67.9029029
  68.02702703  68.15115115  68.27527528  68.3993994   68.52352352
  68.64764765  68.77177177  68.8958959   69.02002002  69.14414414
  69.26826827  69.39239239  69.51651652  69.64064064  69.76476476
  69.88888889  70.01301301  70.13713714  70.26126126  70.38538539
  70.50950951  70.63363363  70.75775776  70.88188188  71.00600601
  71.13013013  71.25425425  71.37837838  71.5025025   71.62662663
  71.75075075  71.87487487  71.998999    72.12312312  72.24724725
  72.37137137  72.4954955   72.61961962  72.74374374  72.86786787
  72.99199199  73.11611612  73.24024024  73.36436436  73.48848849
  73.61261261  73.73673674  73.86086086  73.98498498  74.10910911
  74.23323323  74.35735736  74.48148148  74.60560561  74.72972973
  74.85385385  74.97797798  75.1021021   75.22622623  75.35035035
  75.47447447  75.5985986   75.72272272  75.84684685  75.97097097
  76.0950951   76.21921922  76.34334334  76.46746747  76.59159159
  76.71571572  76.83983984  76.96396396  77.08808809  77.21221221
  77.33633634  77.46046046  77.58458458  77.70870871  77.83283283
  77.95695696  78.08108108  78.20520521  78.32932933  78.45345345
  78.57757758  78.7017017   78.82582583  78.94994995  79.07407407
  79.1981982   79.32232232  79.44644645  79.57057057  79.69469469
  79.81881882  79.94294294  80.06706707  80.19119119  80.31531532
  80.43943944  80.56356356  80.68768769  80.81181181  80.93593594
  81.06006006  81.18418418  81.30830831  81.43243243  81.55655656
  81.68068068  81.8048048   81.92892893  82.05305305  82.17717718
  82.3013013   82.42542543  82.54954955  82.67367367  82.7977978
  82.92192192  83.04604605  83.17017017  83.29429429  83.41841842
  83.54254254  83.66666667  83.79079079  83.91491491  84.03903904
  84.16316316  84.28728729  84.41141141  84.53553554  84.65965966
  84.78378378  84.90790791  85.03203203  85.15615616  85.28028028
  85.4044044   85.52852853  85.65265265  85.77677678  85.9009009
  86.02502503  86.14914915  86.27327327  86.3973974   86.52152152
  86.64564565  86.76976977  86.89389389  87.01801802  87.14214214
  87.26626627  87.39039039  87.51451451  87.63863864  87.76276276
  87.88688689  88.01101101  88.13513514  88.25925926  88.38338338
  88.50750751  88.63163163  88.75575576  88.87987988  89.004004
  89.12812813  89.25225225  89.37637638  89.5005005   89.62462462
  89.74874875  89.87287287  89.996997    90.12112112  90.24524525
  90.36936937  90.49349349  90.61761762  90.74174174  90.86586587
  90.98998999  91.11411411  91.23823824  91.36236236  91.48648649
  91.61061061  91.73473473  91.85885886  91.98298298  92.10710711
  92.23123123  92.35535536  92.47947948  92.6036036   92.72772773
  92.85185185  92.97597598  93.1001001   93.22422422  93.34834835
  93.47247247  93.5965966   93.72072072  93.84484484  93.96896897
  94.09309309  94.21721722  94.34134134  94.46546547  94.58958959
  94.71371371  94.83783784  94.96196196  95.08608609  95.21021021
  95.33433433  95.45845846  95.58258258  95.70670671  95.83083083
  95.95495495  96.07907908  96.2032032   96.32732733  96.45145145
  96.57557558  96.6996997   96.82382382  96.94794795  97.07207207
  97.1961962   97.32032032  97.44444444  97.56856857  97.69269269
  97.81681682  97.94094094  98.06506507  98.18918919  98.31331331
  98.43743744  98.56156156  98.68568569  98.80980981  98.93393393
  99.05805806  99.18218218  99.30630631  99.43043043  99.55455455
  99.67867868  99.8028028   99.92692693 100.05105105 100.17517518
 100.2992993  100.42342342 100.54754755 100.67167167 100.7957958
 100.91991992 101.04404404 101.16816817 101.29229229 101.41641642
 101.54054054 101.66466466 101.78878879 101.91291291 102.03703704
 102.16116116 102.28528529 102.40940941 102.53353353 102.65765766
 102.78178178 102.90590591 103.03003003 103.15415415 103.27827828
 103.4024024  103.52652653 103.65065065 103.77477477 103.8988989
 104.02302302 104.14714715 104.27127127 104.3953954  104.51951952
 104.64364364 104.76776777 104.89189189 105.01601602 105.14014014
 105.26426426 105.38838839 105.51251251 105.63663664 105.76076076
 105.88488488 106.00900901 106.13313313 106.25725726 106.38138138
 106.50550551 106.62962963 106.75375375 106.87787788 107.002002
 107.12612613 107.25025025 107.37437437 107.4984985  107.62262262
 107.74674675 107.87087087 107.99499499 108.11911912 108.24324324
 108.36736737 108.49149149 108.61561562 108.73973974 108.86386386
 108.98798799 109.11211211 109.23623624 109.36036036 109.48448448
 109.60860861 109.73273273 109.85685686 109.98098098 110.10510511
 110.22922923 110.35335335 110.47747748 110.6016016  110.72572573
 110.84984985 110.97397397 111.0980981  111.22222222 111.34634635
 111.47047047 111.59459459 111.71871872 111.84284284 111.96696697
 112.09109109 112.21521522 112.33933934 112.46346346 112.58758759
 112.71171171 112.83583584 112.95995996 113.08408408 113.20820821
 113.33233233 113.45645646 113.58058058 113.7047047  113.82882883
 113.95295295 114.07707708 114.2012012  114.32532533 114.44944945
 114.57357357 114.6976977  114.82182182 114.94594595 115.07007007
 115.19419419 115.31831832 115.44244244 115.56656657 115.69069069
 115.81481481 115.93893894 116.06306306 116.18718719 116.31131131
 116.43543544 116.55955956 116.68368368 116.80780781 116.93193193
 117.05605606 117.18018018 117.3043043  117.42842843 117.55255255
 117.67667668 117.8008008  117.92492492 118.04904905 118.17317317
 118.2972973  118.42142142 118.54554555 118.66966967 118.79379379
 118.91791792 119.04204204 119.16616617 119.29029029 119.41441441
 119.53853854 119.66266266 119.78678679 119.91091091 120.03503504
 120.15915916 120.28328328 120.40740741 120.53153153 120.65565566
 120.77977978 120.9039039  121.02802803 121.15215215 121.27627628
 121.4004004  121.52452452 121.64864865 121.77277277 121.8968969
 122.02102102 122.14514515 122.26926927 122.39339339 122.51751752
 122.64164164 122.76576577 122.88988989 123.01401401 123.13813814
 123.26226226 123.38638639 123.51051051 123.63463463 123.75875876
 123.88288288 124.00700701 124.13113113 124.25525526 124.37937938
 124.5035035  124.62762763 124.75175175 124.87587588 125.        ]

allcvals_modified = allcvals[5:-5] #we are dropping candidate c values that are very near the edges
print(allcvals_modified)
logpostvals = np.array([logpost(c) for c in allcvals_modified])
print(logpostvals)
plt.plot(allcvals_modified, logpostvals)
plt.xlabel('Change of Slope Point')
plt.ylabel('Value')
plt.title('Logarithm of (unnormalized) posterior density')
plt.show()

Fetching long content....

Next we normalize the posterior density.

postvals_unnormalized = np.exp(logpostvals - np.max(logpostvals))
postvals = postvals_unnormalized/(np.sum(postvals_unnormalized))
plt.plot(allcvals_modified, postvals)
plt.xlabel('Changepoint')
plt.ylabel('Probability')
plt.title('Posterior distribution of change points')
plt.show()

Because of the presence of the term $|X_c^T X_c|^{-1/2}$ , the maximizer of the posterior might be different from the least squares estimator (which minimizes RSS). However the two estimators will be close to each other:

print(allcvals_modified[np.argmax(postvals)]) #this is the posterior mode or maximizer
print(c_hat) #this is the least squares estimator

66.28928928928929
66.28928928928929

Using this posterior distribution, we can compute a 95% uncertainty interval for $c$ in the following way. We calculate the probability of regions around the posterior maximizer, and find the region with probability at least 0.95. The code for this is given below.

#95% credible interval for c:
def PostProbAroundMax(m):
    est_ind = np.argmax(postvals)
    ans = np.sum(postvals[(est_ind-m):(est_ind+m)])
    return(ans)
m = 0
while PostProbAroundMax(m) <= 0.95:
    m = m+1
est_ind = np.argmax(postvals)
c_est = allcvals_modified[est_ind]
#95% credible interval for f:
ci_c_low = allcvals_modified[est_ind - m]
ci_c_high = allcvals_modified[est_ind + m]
print(np.array([c_est, ci_c_low, ci_c_high]))
print(np.array([c_est - 1+ tme[0], ci_c_low - 1 + tme[0], ci_c_high - 1 + tme[0] ]))

[66.28928929 64.3033033  68.27527528]
[1965.28928929 1963.3033033  1967.27527528]

Below we draw posterior samples for $c$ , and plot the posterior samples along with the data.

#Drawing posterior samples for c: 
N = 4000
rng = np.random.default_rng(seed = 42)
cpostsamples = rng.choice(allcvals_modified, N, replace = True, p = postvals)

#Let us plot the posterior samples for c on the original data:
plt.figure(figsize = (8, 7))
plt.plot(tme, y)
for i in range(N):
    plt.axvline(x = cpostsamples[i] - 1 + tme[0], color = 'gray')
plt.plot(tme, y, color = 'blue')
plt.axvline(x = c_hat - 1 + tme[0], color = 'black')
plt.show()