How to automatically test for the best parameters for transformed independent variable in linear model
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[]&lt;&gt;"LHrXQ.png"]
img2 = ImageResize[ImageTake[img, {40, 450}, {50, 1000}], 200]
pixels = {(#[[1]]-10)*3.8,#[[2]]*9+100}&amp;/@PixelValuePositions[img2,Black, 0.4];

ListPlot[pixels, Frame-&gt;True, ImageSize-&gt;500]




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,
                                         domain=domain,
                                         model_type='GP',
                                         acquisition_type='EI',
                                         initial_design_numdata=10,
                                         initial_design_type='random',
                                         acquisition_jitter=0.01,
                                         num_cores=-1,
                                         de_duplication=True,
                                         exact_feval=True)

# Run optimization
bo.run_optimization(max_iter=100)


We can plot the optimizers convergence:

bo.plot_convergence()




And information about how it sampled the parameters:

bo.plot_acquisition() 




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).