Fractal interpolation has been applied to reconstruct a continuous irregular function from its finite set of samples. Constructing a linear fractal interpolation function is the simplest way to obtain a fractal interpolation function and interpolate given sample data. This paper aims to provide a modified model of linear fractal interpolation functions to improve the quality of data fitting and keep the construction as simple and easy to apply as possible. We split the given raw data set into a training data set and a testing data set. A set of sample data is chosen to construct fractal interpolation functions. An objective function is defined and the training data are used to find the optimal values of parameters to minimize the objective function. The mean squared error on the testing data set is calculated to evaluate the fitness of the obtained fractal function to “unseen” data. The complexity of the raw data set and the obtained fractal functions are compared by their fractal dimensions. Optuna is used to solve the parameter identification problem. According to the results obtained in the examples, we see that the more data used for training, the smaller the test mean squared error will be. The fractal interpolation function established by the sum of the smooth regression function constructed from the sample data set and the linear fractal interpolation function constructed from the residual of the regression has the best performance. The effect of the difference between the fractal dimension of the training data set and that of the obtained fractal function on the test mean squared error is also discussed.