“In theory there is no difference between theory and practice. In practice, there is.”
Yogi Berra
By “motivated curve fitting”, we are going to focus on two different types of motivated curve fitting. The first is a from the ground up modeling approach, where we try and predict what we will see based on what we know about the world. The second is to use more statistical tools to learn from data, with our choice of model being less constricted than the first approach.
For the “ground up approach”, this stems from the idea that we have an underlying mechanism we think exists that motivates the functions we use to match our data. In physics, these mechanisms are derived from “first principles”, which means we logically build a model from the ground up based on what we know about the real world. An example of these are Newton’s laws, derived from a “falling apple” (supposedly).
Alternatively, we could fit a statistical model to some sort of experimental data, using a more broad class of functions. To maintain the “motivated” portion, we would like to maintain some constraints or prior information (such as when modeling a disease spread we cannot have a negative population).
In an ideal world, the former approach can be considered a more interpretable model, as the parameters in the model are causal and derived from first principles. Intervention in a system via changing a parameter will change the output in an expected way. From a more detrimental perspective, the motivated curve fitting paradigm could be viewed as constrained statistical modeling. This view abandons the generosity and grace of the first worldview and looks at parameterized models as constrictions of the data fitting process to only consider certain types of functions. If we know a process should only reside in these spaces, great! If not, well then, to be frank, we have really limited the predictive capabilities of our model for no real good reason. And the causality we infer may not even be accurate. TLDR, interpretability is great, but our model may not be realistic, in which case, sure we understand causality but it may not be right. And constraints are great…again if they are actual constraints!
The following fun image helps illustrate these approaches.
The image describes airbnb bookings during the 2024 American solar eclipse. Building a model of the moon’s orbit, even a simplification like Newton’s laws, can predict where the moon will be with great precision. Knowing this led astronomers to relay this information to the greater public, which led to massive public interest in this unique and seemingly biblical event.
On the other hand, if an analyst were studying this data after the fact, they’d probably struggle to explain what was going on. Without knowing which features to include in their model, they may not be able to learn the signal in the data, and would probably not be able to explain the pattern they are seeing. Simply include the location of the moon, an example of clever feature engineering could provide insight into what the pattern is. While the analyst need not be familiar with the detailed physics under the hood, having the intuition to include physical data could provide both increased predictive capabilities as well as explanatory insights to the analyst. Otherwise the job becomes searching for a needle in a haystack.
If a model that dictates a certain functional form relating inputs to outputs is justified, then the approaches suggested here are very useful. For example, we understand governing dynamics in weather patterns fairly well1. However, we do not always know the functional form of interest or anything particularly close. This is part of the reason we tend to advocate for non-parametric models in later chapter. We can still glean statistical properties of certain estimands, whether or not a specific “curve” is pre-specified to fit the data. In the meantime, read this article on non-parametric identification. We replicate some insightful excerpts below:
1 Unfortunately, the systems are non-linear dynamical systems which are subject to chaotic behavior, impeding the ability to forecast out very far with any confidence.
Why does this work? Because we’ve specified a relationship between the behavior of Y conditional on Y being positive and the way Y behaves when it’s negative. The data doesn’t reveal this relationship. Our parametric assumption imposes it.
Maybe this model is motivated by some reasonable model of the world, but it’s probably not. It’s just a functional form. There is a mapping between the data and the parameter of interest because the functional form assumptions fill the gap between the data and what we want to know.
So, let’s look at why nonparametric identification is so important:
When we observe (\(Y,X\)), we are nonparametrically identified. Even if we just run a linear model in estimation, we know that the regression is trying to fit variation that does precisely reveal what we want to know, independent of the functional form assumption.
When we only observe (\(Y*,X\) ), \(t(x)\) (which is the conditional mean \(E(Y\mid X=x)\)) is not nonparametrically identified. Our functional form tells us where the parameter of interest is in the data. If we choose another parametric form, another part of the data will determine the parameter of interest. The functional form determines the variation we use to get at the parameter of interest.
Nonparametric identification matters, and we should doubt analysis based on parametrically identified parameters.
…
Parametric estimation is fine, but our models should not require the specific parametric form for the data to be informative about the parameter of interest.
Zach Flynn
While we discuss motivated curve fitting in this section, most of this book will not espouse it as the end all be all. Generally, we will show a strong preference for non-parametric modeling, which we will discuss in more detail later. Why should we prefer a non-parametric approach? 1) We do not always know the rules of how the variables behave together. 2) If we can collect data on more variables that could plausibly be impactful, we should. Just because we cannot derive the rules of the game for these variables, does not mean we should abandon all hope. In fact, this is where letting the “machine” sort out the relationships is a really handy device
Draw a spring going back and forth
First, simulate a spring that eventually dies out (think a slinky). We then fit the data using an equation derived from a physics model (solved using Green’s functions). Software will pick the parameters that align the fit of the curve to the data as well as possible. In this case, knowing the physics of the model allows us to extrapolate our fit out beyond the “training” region.
import numpy as np
import matplotlib.pyplot as plt
N = 250
mu = 1
m = 0.2
k = 200
gamma = mu/(2*m)
t=np.linspace(0,1,N)
t_test=np.linspace(0,2,2*N)
def spring_eq(gamma, k, t):
w = np.sqrt(k/m)
omega = np.sqrt(w**2-gamma**2)
phi = np.arctan(-gamma/omega)
A = 1/(2*np.cos(phi))
x = 2*np.exp(-gamma*t)*A*np.cos(phi+omega*t)
return(x)
plt.scatter(t, spring_eq(gamma, k,t), color='#073d6d')
plt.show()
Now, we use “curve fit” from scipy. Of course, we are fitting the actual equation and our measurements are not at all corrupted by noise, so this is very much a proof of concept and a validation of the scipy curve fit software.
from scipy.optimize import curve_fit
def fit_func(t,gamma, A, phi, omega):
return(2*np.exp(-gamma*t)*A*np.cos(phi+omega*t))
popt, pcov = curve_fit(fit_func,t, spring_eq(gamma, k,t))
plt.scatter(t,spring_eq(gamma, k,t), color='#073d6d')
plt.plot(t_test, fit_func(t_test, *popt), color='#55AD89', linewidth=2.75)
plt.show()
A famous example of a curve fitting is the windchill equation. It is the equation that determines what the temperature feels like for someone given how windy it is. This equation is almost certainly not as rigorously derived as the damped harmonic oscillator, which is a well established and tested model which was spun up from first principles. The windchill equation, on the other hand, more resembles a linear regression, but with certain transformations of the input variables included that may or may not be physically motivated. The equation that is used by weather services is (with temperature in Fahrenheit and wind in miles per hour):
\[\begin{align}\text{windchill} &= 35.74 + 0.6215*\text{temperature}- \\ &35.75*\text{wind}^{0.16} +0.4275*\text{temperature}^{0.16} \end{align} \]
Let’s fit the equation directly in R:
options(warn=-1)
suppressMessages(library(NatParksPalettes))
suppressMessages(library(tidyverse))
temp <- rep(seq(from=-50, to=50, length.out = 15), 15)
winds <- seq(from=-5, to=50, length.out=15)
wind <- unlist(lapply(1:15, function(i)rep(winds[i], 15)))
vals <- unique(scales::rescale(c(volcano)))
o <- order(vals, decreasing = FALSE)
cols <- scales::col_numeric("PuOr", domain = NULL)(vals)
colz <- setNames(data.frame(vals[o], cols[o]), NULL)
wc_df = data.frame(temp = temp, wind = wind)
wc_df$windchill <- 35.74+ 0.6215*wc_df$temp - 35.75*(wc_df$wind^(0.16))+0.4275*wc_df$temp*wc_df$wind^(0.16)
wc_df = data.frame(temp = temp, wind = wind)
wc_df$windchill <- 35.74+ 0.6215*wc_df$temp - 35.75*(wc_df$wind^(0.16))+0.4275*wc_df$temp*wc_df$wind^(0.16)
pal <- natparks.pals("Acadia",n=9,type="discrete")
pal <- pal[c(1, 2, 4, 6, 8,9)]
wc_true = ggplot(wc_df, aes(temp, wind, fill= windchill))+
geom_tile()+scale_fill_gradientn(colors=pal)+#natparks.pals("Acadia"))+
geom_text(data=wc_df,aes(x=temp, y=wind,label=round(windchill,0)),
size=4, color='#f8f9fa', fontface='bold')+
scale_x_continuous(name="Temperature", limits=c(-55, 55),
breaks=round(seq(from=-50, to=50, length.out=21)))+
scale_y_continuous(name="Wind (mph)", limits=c(0, 55),
breaks=round(seq(from=0, to=50, length.out=11)))+
theme_minimal()A more interesting plot might look at the difference in temperature between the actual and the feels like at each wind level:
wc_diff = ggplot(wc_df, aes(temp, wind, fill= windchill-temp))+
geom_tile()+scale_fill_gradientn(colors=pal)+#natparks.pals("Acadia"))+
geom_text(data=wc_df,aes(x=temp, y=wind,label=round(windchill-temp,0)),
size=4, color='#f8f9fa', fontface='bold')+
scale_x_continuous(name="Temperature", limits=c(-55, 55),
breaks=round(seq(from=-50, to=50, length.out=21)))+
scale_y_continuous(name="Wind (mph)", limits=c(0, 55),
breaks=round(seq(from=0, to=50, length.out=11)))+
theme_minimal()
wc_diff
For an interactive version of the plot:
suppressMessages(library(reactablefmtr))
wind = c(5,10,15,20,25,30,35,40,45,50)
temp = c(-20,-15,-10,-5,0,5,10,15,20,25,30,35,40,45,50)
windchill = matrix(NA, ncol=length(temp),nrow=length(wind))
for (j in 1:length(wind)){
for (k in 1:length(temp)){
windchill[j,k] = 35.74+ 0.6215*temp[k] -
35.75*(wind[j]^(0.16))+0.4275*temp[k]*wind[j]^(0.16)
}
}
windchill = as.data.frame(windchill)
colnames(windchill) = paste0(temp, '(\u00B0F)')
#temp
windchill = cbind(wind, windchill)
colnames(windchill)[1] = 'wind'
reactable(
windchill,
theme = espn(),
compact = TRUE,
pagination = FALSE,
showSortIcon = T,
defaultSorted = "wind",
defaultSortOrder = 'asc',
defaultColDef = colDef(
maxWidth = 55,
cell = color_tiles(windchill,
number_fmt = scales::comma,
colors = MetBrewer::met.brewer('Morgenstern',8),
text_size = 22, span = TRUE)
),
columns = list(
wind = colDef(maxWidth = 70)
)
) %>%
add_title('Windchill calculator') %>%
add_source('Demetrios Papakostas: https://www.weather.gov/safety/cold-faqs')Registered S3 method overwritten by 'MetBrewer':
method from
print.palette NatParksPalettesDemetrios Papakostas: https://www.weather.gov/safety/cold-faqs
The SIR model, which stands for “susceptible, infected, and recovered” dates back almost 100 years. It’s early uses included modeling the spread of disease in a population, an all too familiar concept in the post 2020 world.
The SIR model is a nice mechanistic model that generally matches our expectation for how a disease would spread. The dynamics of the spread of a various are governed by 3 “compartments” of people in a population, susceptible, infected, and recovered individuals. People get infected, spread the disease to other people, become immune to the disease, and move into a recovered compartment. Eventually, the disease runs out of people to spread and dissipates.
This model is nice because it is intuitively sound and relatively simple. Real life infection curves mostly look like the shapes that solutions to the SIR model grant. The model can be studied mathematically and the parameters have some meaning. If we change a parameter, then the outcome will change in ways we can easily test2.
2 Just wait about a paragraph though.
Issues:
The base SIR model is fairly simple and may not capture the necessary interactions. As we will see in the simulations below, the shapes permitted from the SIR model are not extremely flexible; they are mono-modal and generally symmetric or with slight skews for realistic \(R_0\) values (a measure of disease infectiousness). More complex models can capture more complex behaviours, but they are adding more assumptions about the correct formulation of the disease spread. For example, including compartments to attempt to account for human behaviour seems like a good idea. But nailing the form of that behaviour is hard and doing it a-priori reduces the flexibility of a model. Say this new model assumes people stay home if the disease spreads quickly. We know from experience with COVID that this varied a lot!
Additionally, the disease itself evolves in (more predictable) ways. Therefore, we can extend the SIR model to include multiple strains (Pell et al. 2023), but this requires careful thought about interactions between strains and consideration of immune dynamics, and maybe iffy assumptions. It’s unwieldy with more than two strains. (Dukic, Lopes, and Polson 2012) embed the SIR into a “state-space” modeling regime which allows for parameters to vary over time and for noise in the observations to be modeled explicitly. In fact, many of the best performing COVID forecasts incorporated this approach, including (Gibson, Reich, and Sheldon 2023), who had the cleanest webpage I think I have ever seen and the model of Youyang Gu (Gu 2020).
The question arises if the implied constraints of the SEIR (which expands the SIR by including an “exposed” class) lead to the quality forecasts. In one sense, seir constraints seem logically generally valid, if somewhat simple. The SEIR implies a skewed unimodal curve, which is a reasonable constraining condition. However, the SEIR is not the only data generating system to yield curves that look like! Many functions look like that, such as the functions generating beta or log-normal distributions. And those may be preferable to the SEIR formulation, which stipulates conditions/rules/constraints (whatever you wanna call them) between unobserved things3 which makes it almost impossible to grok if those conditions/rules/constraints are actually valid in practice.
3 For example, we do not actually observe a “susceptible” or “exposed” class, and arguably do not know the infected or recovered class either due to under-representation of cases relative to true disease spread.
Some dynamical system models, such as the fox-rabbit model are defined in terms of observable variables. In these models, the “causal levers” that the parameters define in the real world are more meaningful. Additionally, the simplicity of this model makes practical identification more likely. However, to avoid the over-simplification of the fox-rabbit model makes it unlikely to match real world data. Incorporating more complexity through a model that defines the structures of the data explicitly means that we’d likely have to delve into unobserved variables again or make strong assumptions about behaviours of the animals and the environment around them.
Additionally, the fitting of dynamical systems can be challenging, particularly if a system consists of partial differential equations, which we will ignore here. Some parameters are inherently correlated (deaths affect the contact rate) and the formulation/fitting assume independence. As you add more parameters, they are less likely to be identifiable. The procedure generally becomes more prone to overfitting as well.
Finally, data quality also affects our analysis. This is a problem that we will become well aquainted with. For example, for a a study on different strains, validation is thrown out the window if detailed data on strains (and thus waves) are not available.
Using a dot to denote differentiation with respect to time, the SIR model with standard incidence is given by \[ \begin{align} \dot{S}&= -\tilde{\beta} S I,\\ \dot{I}&= \tilde{\beta} S I-\tilde{\gamma} I,\\ \dot{R}&= \tilde{\gamma} I. \end{align} \]
The model stipulates people go from susceptible class (\(S\)) to infectious (\(I\)) to recovered (\(R\)), with there being a fixed number of people \(N\) in the population equaling the sum of these three compartments. The simplicity of the standard SIR model implies several assumptions. The population is closed, and the model as written in the above system of equations does not account for demographics or disease-induced mortality. Furthermore, individuals have equal probability of coming into contact with one another, which is often an inaccurate reflection of contact networks. It should be noted that the waiting time in each compartment is exponentially distributed; for example, the average duration of infectiousness in the SIR formulation s \(1/\tilde{\gamma}\) days.
In code, simulating from the SIR can be done as follows:
options(warn=-1)
## Load deSolve package
suppressMessages(library(deSolve))
## Create an SIR function
sir <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
dS <- -beta * S * I
dI <- beta * S * I - gamma * I
dR <- gamma * I
return(list(c(dS, dI, dR)))
})
}
### Set parameters
## Proportion in each compartment: Susceptible 0.999999, Infected 0.000001, Recovered 0
init <- c(S = 1-1e-6, I = 1e-6, R = 0.0)
## beta: infection parameter; gamma: recovery parameter
parameters <- c(beta = 1.4247, gamma = 0.14286)
parameters <- c()
## Time frame
times <- seq(0, 100, by = .1)
out <- c()
R0_list <- c()
n_vals <- 8
for (i in 1:n_vals){
parameters[[i]] <- c(beta = seq(from=.5, to=1, length.out=n_vals)[i], gamma =
seq(from=0.1, to=.75, length.out=n_vals)[i])
R0_list[[i]] <- parameters[[i]][1]/parameters[[i]][2]
parameter2 <- parameters[[i]]
## Solve using ode (General Solver for Ordinary Differential Equations)
out[[i]] <- ode(y = init, times = times, func = sir, parms = parameter2)
}
out <- sapply(1:n_vals, function(i)out[[i]][,3])
## change to data frame
out <- as.data.frame(out)
## Delete time variable
out$time <- NULL
suppressMessages(library(NatParksPalettes))
#use Arches, Acadia, or Yosemite
pal <- natparks.pals("Acadia", n=10)
matplot(x = times, y = out[,], type = "l",
xlab = "Time", ylab = "Infected", main = "SIR with different R0",
lwd = 4, lty = 1, bty = "l", col = pal)
## Add legend
legend('topright', legend=round(unlist(R0_list),2),lty=1,lwd=4,
col = pal, bty = "n", title='R0')
Simulate data using this code from stackover flow ellipse generator answer, and using the convenient cartesian to polar translation equations given.
options(warn=-1)
N = 250
x_0 = 4
y_0 = -4
a = 1 # major axis length
b = 2 # minor axis length
phi = pi # angle of major axis with x axis
eps = rnorm(N, 0 ,0.05) #noise
t = seq(from=0, to=2*pi, length.out = N)
x = x_0 + a*cos(t)*cos(phi) - b*sin(t)*sin(phi) + eps
y = y_0 + a*cos(t)*cos(phi) + b*sin(t)*cos(phi) + eps
plot(x,y, col='#012296', pch=16)
We need to convert to Cartesian coordinates to minimize the equation of the ellipse:
\[ ax^2+by^2+cx*y+dx+ey+1=0 \]
where we set \(f=1\) as a constraint so that not all parameters are just 0. To check the fit, you can sum the fitted coefficients and check they add to 0:
options(warn=-1)
mod = lm(zero ~ x2+xy+y2+x+y-1, data=data.frame(x2=x^2,xy=x*y,y2=y^2, x,y,zero=rep(-1, length(y))),singular.ok = F)
a = mod$coefficients[1]
b = mod$coefficients[2]
c = mod$coefficients[3]
d = mod$coefficients[4]
e = mod$coefficients[5]
f = 1
# check the solution sums to 0
print(paste0('The sum of the fitted coefficients is: ', round(sum(a*x^2+c*y^2+b*x*y+d*x+e*y+1),3)))[1] "The sum of the fitted coefficients is: 0.001"The following code converts the cartesian coordinates back to polar. Admittedly, this is copied from this stackoverflow link (the same one as earlier), so the code is hidden since the details are beyond the scope of these notes because they require more advanced linear algebra.
#| echo: true
#| code-fold: true
#| code-summary: "expand for full code"
#https://stackoverflow.com/questions/41820683/how-to-plot-ellipse-given-a-general-equation-in-r
plot.ellipse <- function (a, b, c, d, e, f, n.points = 1000) {
## solve for centre
A <- matrix(c(a, c / 2, c / 2, b), 2L)
B <- c(-d / 2, -e / 2)
mu <- solve(A, B)
## generate points on circle
r <- sqrt(a * mu[1] ^ 2 + b * mu[2] ^ 2 + c * mu[1] * mu[2] - f)
theta <- seq(0, 2 * pi, length = n.points)
v <- rbind(r * cos(theta), r * sin(theta))
## transform for points on ellipse
z <- backsolve(chol(A), v) + mu
x_fit = z[1,]
y_fit = z[2,]
df = data.frame(x_fit, y_fit)
return(df)
}
vals = plot.ellipse(a,c,b,d,e,1)
x_fit = vals$x_fit
y_fit = vals$y_fit
plot(x,y, pch=16, col='#012296')
lines(x_fit, y_fit, lwd =3,col='#FD8700')
In this subsection, we will try and show how to fit a flexible model to data that is still well motivated. To illustrate this idea, we use the “optim” function in R, which allows for a variety of optimization techniques in an easy-ish way. We simulate from a sine curve with noise to mimic seasonal precipitation totals in a seasonal state in the U.S., but fit with a 4-parameter Beta distribution, which is nice because the parameters are meaningful. In particular, we can control the total amount of the outcome contributed in each curve, the shape of the curves, the number of curves, and the start and end points. We manually set the number of points, \(m\), from a visual inspection of the data. The data also guide our initial estimates of the \(\alpha\), \(\beta\), and start and end parameters for each beta curve.
Unsurprisingly, initial estimates are very important. If we start in a local minimum, we are doomed from the get-go. We first use the “Nelder-Mead” (Nelder and Mead 1965) approach to find the minimum/maximum (depending on the form of the function) of the objective function4. We will not go into detail here, but it is a simplex based method (i.e. a point in 0 dimensions, a line in 1 dimension, a triangle in 2 dimensions, a triangular pyramid in 3 dimensions, and so on). Broadly speaking, the algorithm orders the outputs \(f(x_i)\) for \(i\in \{1,2, \ldots, d+1\}\), where \(d\) is the number of dimensions. The points are ordered by what the evaluation of \(f(x)\), such that \(f(x_1)\) is the “best” evaluation (the closest to the minimum/maximum amongst the points evaluated) and \(f(x_{d+1})\) is the worst (farthest from the minimum/maximum amongst the points evaluated).
4 The mean squared error between the observed data and the “mixture” of beta-functions model we propose.
The simplex is initialized by the \(d+1\) \(x_i\) points (recall, in two dimensions this is a triangle) . We begin by calculating the centroid of all points but the maximum \(x_i\). The centroid is crucial for the calculations of the next steps. These next steps involve going through reflections, expansions, or contractions around the point \(x_{d+1}\), evaluate the new point “\(r\)” $, which is \(f(x_{r})\), and evaluate some criterion for whether or not that move is accepted. All test points are “shrunk” except the best (\(x_1\)). See the wikipedia page for more detail and cool visuals.
options(warn=-1)
suppressMessages(library(gt))
set.seed(012024)
m = 4
N = 200
x = seq(from=1, to=200, length.out = N)
y = ifelse(x<100, 10*sin(x/16)^2, 20*sin(x/16)^2) + rnorm(N, 0, 2)
y = y/max(y)
dbeta4param <- function(x,shp1,shp2,bnd1,bnd2){
mask = (x>=bnd1 & x<=bnd2)
xScale = (x-bnd1)/(bnd2-bnd1)
mask * dbeta(xScale,shape1=shp1,shape2=shp2)/(bnd2-bnd1)
}
f_obj = function(theta){
fit <- rowSums(sapply(1:m, function(j)
theta[j+4*m]*dbeta4param(x, theta[j],theta[j+m]*theta[j], theta[j+2*m], theta[j+3*m])))
return(sum((fit-y)^2))
}
alpha_shape = rep(1,m);alpha_lower = rep(1,m);alpha_upper = rep(50,m)
beta_shape = rep(1,m);beta_lower = rep(1,m);beta_upper = rep(50,m);
start = c(1,51, 101, 151); start_lower = rep(1,m); start_upper = rep(200, m);
end = c(50, 100, 150, 225); end_lower = rep(20, m); end_upper = rep(350, m);
scale = c(0.5,0.5,1,1); scale_lower = rep(.1,m); scale_upper = rep(1000,m);
outs = optim(c(alpha_shape,beta_shape,start,end, scale), f_obj,
method='Nelder-Mead')
theta_beta = outs$par
preds <- rowSums(sapply(1:(m), function(j)
theta_beta[j+4*m]*dbeta4param(x,theta_beta[j],theta_beta[j+m]*theta_beta[j], theta_beta[j+2*m], theta_beta[j+3*m])))
plot(x,y, pch=16, col='#012296', xlab='time', ylab = 'rainfall totals', cex=1.25, main='Nelder Mead')
lines(x, preds, col="#FD8700", lwd=4)
Given the function we are fitting, constraints our important because of the allowed inputs to the modified beta function we are fitting. So now we try the ‘L-BFGS-B’ method, which allows for constraints on parameters encoded through intervals from which potential values are explored from. This makes use of a modification of the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm, a method similar in spirit to Newton’s method and put forward in these 4 seemingly independent (!) papers in 1970, see (Broyden 1970), (Fletcher 1970), (Goldfarb 1970), and (Shanno 1970).
outs = optim(c(alpha_shape,beta_shape,start,end, scale), f_obj,
lower=c(alpha_lower, beta_lower, start_lower, end_lower, scale_lower),
upper=c(alpha_upper, beta_upper, start_upper, end_upper, scale_upper),
method='L-BFGS-B')
theta_beta = outs$par
preds <- rowSums(sapply(1:(m), function(j)
theta_beta[j+4*m]*dbeta4param(x,theta_beta[j],theta_beta[j+m]*theta_beta[j], theta_beta[j+2*m], theta_beta[j+3*m])))
plot(x,y, pch=16, col='#012296', xlab='time', ylab = 'rainfall totals', cex=1.25, main='BFGS-B')
lines(x, preds, col="#FD8700", lwd=4)
data.frame(starts=round(theta_beta[9:12],3), ends=round(theta_beta[13:16],3), totals = round(theta_beta[17:20],3))%>%
gt()| starts | ends | totals |
|---|---|---|
| 2.034 | 52.972 | 10.874 |
| 50.207 | 100.115 | 10.940 |
| 99.171 | 152.279 | 20.328 |
| 144.300 | 223.874 | 21.440 |