Let's assume that I have a linear model with $k$ variables:

$y = \beta_0 + \beta_1\cdot x_1 + \dots + \beta_k \cdot x_k$.

Now, I want to add variable $x_{k+1}$, but, according to domain knowledge, the dependency of $y$ onto $x_{k+1}$ is not linear, but rather "S-shaped". To capture this dependency, a well-established method is to use parametrised arcus tangens function: $D \cdot \arctan{\frac{x_{k+1} - A}{B}} + C$. To include it in the linear model we can safely ignore $D$ and $C$ parameters (as they will contribute to $\beta_{k+1}$ and $\beta_0$ respectively), so eventually I'd end up with the model:

$y = \beta_0 + \beta_1\cdot x_1 + \dots + \beta_k \cdot x_k + \beta_{k+1}\cdot \arctan{\frac{x_{k+1} - A}{B}}$.

What I would like to do is to find a way to test different transformations (with different $A$ and $B$ parameters automatically. Firstly, I went for searching the grid of different $A$ and $B$ values, i.e. fitting the model to different transformations of $x_{k+1}$ and returning list of models, sorted by $R^2$, $RMSE$ or $AIC$. However, this usually favours very sharp curves:

enter image description here

This doesn't make much sense in the eye test though. I know that eye-test may be biased but to me, more propable dependency would be captured by a shape as below:

enter image description here

One could already noticed that there are plenty of observations where $x_{k+1}$ is or is near to $0$, which is actually a characteristic of those variables that will be approximated by $arctan$. That may probably affect the fitting, but still something clearly doesn't work here.

Another approach that I considered is to fit the model using gradient descent search rather than using standard software (like lm from R). The problem is that with this $arctan$ transformation, we don't have a convex objective function anymore.

How can I then automatically test for the best parameters of the $arctan$ transformation? Am I on the right track with either idea or can I approach this problem differently?r

  • $\begingroup$ Did you check generalised additive models? They are extremely flexible in finding a good parameterisation more or less endogenously... so no need for manual feature engineering. $\endgroup$
    – Peter
    Jan 7, 2020 at 12:35

4 Answers 4


One approach would be to use an algorithm designed for non-convex problems like Bayesian optimization. However, if you have already evaluated a fine grid of parameters this is unlikely to offer significant improvement. Here is an example of how you could implement Bayesian optimization for this problem.

First, we need some data. Just for fun let’s extract the data from the images you posted (in brief since this is off topic).

In Mathematica:

img = Import[NotebookDirectory[]<>"LHrXQ.png"]
img2 = ImageResize[ImageTake[img, {40, 450}, {50, 1000}], 200]
pixels = {(#[[1]]-10)*3.8,#[[2]]*9+100}&/@PixelValuePositions[img2,Black, 0.4];

ListPlot[pixels, Frame->True, ImageSize->500]

enter image description here

Now in order to use Bayesian optimization we need to define the objective as a function of parameters A and B. Here we will maximize the R^2 value for the fitted model.

In python:

import pandas as pd
import numpy as np
import random
from sklearn.linear_model import LinearRegression

data = pd.read_csv('mathematica_data.csv')

def objective(params):
    """Whatever you want to do for your regression."""

    # So it works with GpyOpt
    A = params[0][0]
    B = params[0][1]

    temp = data.copy()

    # Transform variable
    xt = [np.arctan((x - A)/B) for x in data['x1'].tolist()]
    temp['x1'] = xt

    # Fit a linear model
    reg = LinearRegression().fit(temp.drop('y', axis=1), temp['y'])

    # Compute scores of interest
    r2 = reg.score(temp.drop('y', axis=1), temp['y'])

    # GPyOpt will minimize so we want -f
    return - r2

Now use GPyOpt to optimize the objective.

import GPyOpt

domain = [{'name': 'A', 'type': 'continuous', 'domain': (-300.0,300.0)},
          {'name': 'B', 'type': 'continuous', 'domain': (0.1,300.0)}]

bo = GPyOpt.methods.BayesianOptimization(f=objective,

# Run optimization

We can plot the optimizers convergence:


enter image description here

And information about how it sampled the parameters:


enter image description here

Of course using only the independent variable extracted from your plot even the best R^2 values suggests there is little to no relationship between arctan((x-A)/B).


You can build an optimization problem.

Your features will be the parameters, then you create a loss function and try to minimize it with gradient descent or brute force (depends on the search space of your problem)

  • $\begingroup$ I actually mentioned that at the end, but please note that this is not a convex problem anymore due to partial derivatives of $arctan$ $\endgroup$
    – jakes
    Jan 7, 2020 at 18:13
  • $\begingroup$ Do you care about optimality? Just do a brute force search across all values and if not possible, maybe genetic algorithm? $\endgroup$ Jan 7, 2020 at 18:27
  • $\begingroup$ It doesn't have to be perfect optimality. But I've been doing grid search as mentioned in the post and the result were clearly wrong. I am not sure, but I thing this is probably what you call brute force search? $\endgroup$
    – jakes
    Jan 7, 2020 at 19:42

I'm not sure if you are bound to the type of model presented in your question. However, an alternative would be to use generalised additive models (GAM), e.g. with regression splines or locale regression. These methods usually give a very good fit with non-linear patterns in $X$ and there is no need to provide parameterization of $X$ so that it is easy to find good models.

Here is an example based on simulated data (R code):

# Generate data
x <- -50:100
y <- 0.001*x^3
df = data.frame(y,x)

# Linear regression
pred_ols = predict(reg_ols, newdata=df)

# GAM with regression splined (df=3)
reg_gam = gam(y~s(x,3), data=df)
pred_gam = predict(reg_gam, newdata=df)

# Find opt. number of splines
for (sp in seq(1:50)){
  gamx=gam(y~s(x,sp), data=df)
  print(mse(y, predict(gamx, newdata=df)))

# Plot prediction and actual data
df2 = data.frame(x,y,pred_ols, pred_gam)
ggplot(df2, aes(x)) +                    
  geom_line(aes(y=y),size=1, colour="red") +  
  geom_line(aes(y=pred_ols),size=1, colour="blue") +
  geom_line(aes(y=pred_gam),size=1, colour="black", linetype = "dashed")

As you can see, the model gives a very good fit to my non-linear function while there is no need to provide parameterization (see figure). See ISL Chapter 7.7 for an applied introduction.

enter image description here


You could use Nonlinear Least Squares, in which one of the regressors is your arctan function with two more parameters to be estimated.

In R, for example:


df <- datasets::airquality

my_atan <- function(x, A, B){atan((x-A)/B)}

nlsLM(Ozone ~ a + b * Temp + c * my_atan(Temp, A, B),
      data = df,
      start = list(a = 0, b = 0, c = 0, A = 0, B = 1))

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.