Кривая бистабильности в питоне

Question

У меня есть график плотности рассеяния между 2 переменными. Давайте предположим, что ось X является ответом, а ось Y - управляющей переменной. Я хочу нарисовать линию бистабильности на графике, чувствительную к разбросу. Линия и график выглядят примерно так: Код, используемый для построения этих графиков (кроме оранжевой линии, которая нарисована от руки) :

x = np.array(data['Response'])
y = np.array(data['Control'])

#histogram definition
bins = [200, 200] # number of bins

# histogram the data
hh, locx, locy = np.histogram2d(x, y, bins=bins)

# Sort the points by density, so that the densest points are plotted last
z = np.array([hh[np.argmax(a<=locx[1:]),np.argmax(b<=locy[1:])] for a,b in zip(x,y)])
idx = z.argsort()
x2, y2, z2 = x[idx], y[idx], z[idx]
plt.figure(1,figsize=(9,8)).clf()
plt.scatter(x2, y2, c=z2, cmap='jet', marker='.', s =50,vmax = 30, vmin = 0)  

#plt.plot([0,250],[0,250],'--',color = 'black', alpha = 0.6)
plt.xlabel('Response', fontsize = '15')
plt.ylabel('Control', fontsize = '15')  ##Change
plt.xticks(fontsize = '14')
plt.yticks(fontsize = '14')
#plt.xlim(7,20)
#plt.ylim(120,180)
#plt.text(60,220,'$R^2=0.78$',fontsize = 15)
cbar =plt.colorbar()
#cbar.ax.tick_params(labelsize=30) 

Кто-нибудь знает, как это сделать ?? Пожалуйста, обратитесь к этой странице для получения дополнительной информации ( здесь ). Я все еще новичок в Python, поэтому мне нужно некоторое время, чтобы понять это.

Этот код как-то связан, но мне все еще трудно расшифровать его.

#equilibria of x(t)

fig = plt.figure(figsize=(15,7))
mpl.rc('font', size = 14)

#define parameters and predefine array variables:
h = 1.; b = 1.; p = 12; r = 0.6; A=[]; XA=[]

EX1 = []; EX2 = []; EX3 = []

a = 0.5
while a <= 0.9:
    GX = []; X = []
    #calculate g(x) for parameter t in a range of (0,2):
    for x in np.linspace(0,2,1000):
        gx = a-b*x + (r*x**p)/(x**p + h**p)
        GX.append(gx); X.append(x)
        if -10**-4 < gx < 10**-4:
            #calculate g'(x) (= dg(x)/dx) to "check" stability:
            Dgx = -b+(r*p*x**(p-1)*h**p)/((x**p + h**p)**2)
            #write x-coordinate of extremum and depending value of parameter a in an array:
            if Dgx > 0: EX2.append([x,a])      #array for unstable extrema
            if Dgx < 0:                        
                if x < 1: EX1.append([x,a])    #array for stable extrema with x<1
                if x > 1: EX3.append([x,a])    #array for stable extrema with x>1
    #plot f(x) for specific values of a:
    if len(str(a)) == 3:        
        plt.subplot(2, 1, 1)
        plt.plot(X,GX,label = 'a = ' + str(a))
    a += 0.0001

#configure plot properties for subplot(2, 1, 1):
plt.subplot(2, 1, 1);
plt.tight_layout(h_pad=5)
plt.title('g(x) for specific values of a')
plt.xlabel('$x$', fontsize=18)
plt.ylabel('$g_{(x)}$', fontsize=18)
#plt.plot(0, 'k')
plt.grid()
plt.text(0.05, -0.9, 'fig. 4', fontsize=18)
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))

#sort equilibria containing arrays in respect to the increasing order of parameter a:
EX1 = list(zip(*sorted(EX1,key=lambda x: x[0]))); EX1x = EX1[1]; EX1y = EX1[0]
EX2 = list(zip(*sorted(EX2,key=lambda x: x[0]))); EX2x = EX2[1]; EX2y = EX2[0]
EX3 = list(zip(*sorted(EX3,key=lambda x: x[0]))); EX3x = EX3[1]; EX3y = EX3[0]

#plot curve representing equilibria of x depending on the variation of parameter a:
plt.subplot(2, 1, 2);
#lower part of the curve (stable extrema):
plt.plot(EX1x,EX1y,color = 'k', linewidth=2, label = 'stable equilibria');
#central part of the curve (unstable extrema):
plt.plot(EX2x,EX2y,color = 'k', linewidth=2, label = 'unstable equilibria', linestyle = 'dashdot');
#upper part of the curve (stable extrema):
plt.plot(EX3x,EX3y,color = 'k', linewidth=2);

#configure plot properties:
plt.title('equilibria of x(t) for specific values of a')
plt.grid();
plt.xlabel('$a$', fontsize=18)
plt.ylabel('$x$', fontsize=18)
plt.text(0.51, 1.4, 'fig. 5', fontsize=18)
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))