Another sigmoidal regression equation question

I posted an earlier version of this yesterday, but I cannot seem to add this version to that posting because someone seems to have closed that posting for editing, so here is the new version in a new posting.

I have the script below that does the following things:

1.) plot a best fit curve to sigmoidal data.

2.) re-size the data based on new max and min coordinates for x and y.

3.) calculate and plot a new best fit curve for the resized data.

Steps 1 and 2 seem to work fine, but step 3 does not. If you run the script, you will see that it plots a totally invalid curve for the re-sized data.

Can anyone show me how to revise the code below so that it creates and plots a true best fit sigmoidal curve for the re-sized data? This needs to be reproducible when re-sized across a spectrum of possible max and min values.

I seem to be able to track the problem to New_p, which is defined in the following line of code:

New_p, New_cov, New_infodict, New_mesg, New_ier = scipy.optimize.leastsq( 
    residuals,New_p_guess,args=(NewX,NewY),full_output=1,warning=True)   

But I cannot figure out how to get deeper into the problem than that. I think the problem may have something to do with the difference between global and local variables, but perhaps it is something else.

Here is the current draft of my complete code:

import numpy as np 
import matplotlib.pyplot as plt 
import scipy.optimize 

def GetMinRR(age):
    MaxHR = 208-(0.7*age)
    MinRR = (60/MaxHR)*1000
    return MinRR

def sigmoid(p,x):
    x0,y0,c,k=p 
    y = c / (1 + np.exp(-k*(x-x0))) + y0 
    return y 

def residuals(p,x,y): 
    return y - sigmoid(p,x) 

def resize(x,y,xmin=0.0,xmax=1.0,ymin=0.0,ymax=1.0):
    # Create local variables
    NewX = [t for t in x]
    NewY = [t for t in y]
    # If the mins are greater than the maxs, then flip them.
    if xmin>xmax: xmin,xmax=xmax,xmin 
    if ymin>ymax: ymin,ymax=ymax,ymin
    #----------------------------------------------------------------------------------------------    
    # The rest of the code below re-calculates all the values in x and then in y with these steps:
    #       1.) Subtract the actual minimum of the input x-vector from each value of x
    #       2.) Multiply each resulting value of x by the result of dividing the difference
    #           between the new xmin and xmax by the actual maximum of the input x-vector
    #       3.) Add the new minimum to each value of x
    # Note: I wrote in x-notation, but the identical process is also repeated for y
    #----------------------------------------------------------------------------------------------    
    # Subtracts right operand from the left operand and assigns the result to the left operand.
    # Note: c -= a is equivalent to c = c - a
    NewX -= x.min()

    # Multiplies right operand with the left operand and assigns the result to the left operand.
    # Note: c *= a is equivalent to c = c * a
    NewX *= (xmax-xmin)/NewX.max()

    # Adds right operand to the left operand and assigns the result to the left operand.
    # Note: c += a is equivalent to c = c + a
    NewX += xmin

    # Subtracts right operand from the left operand and assigns the result to the left operand.
    # Note: c -= a is equivalent to c = c - a
    NewY -= y.min()

    # Multiplies right operand with the left operand and assigns the result to the left operand.
    # Note: c *= a is equivalent to c = c * a
    NewY *= (ymax-ymin)/NewY.max()

    # Adds right operand to the left operand and assigns the result to the left operand.
    # Note: c += a is equivalent to c = c + a
    NewY += ymin

    return (NewX,NewY)

# Declare raw data for use in creating logistic regression equation
x = np.array([821,576,473,377,326],dtype='float') 
y = np.array([255,235,208,166,157],dtype='float') 

# Call resize() function to re-calculate coordinates that will be used for equation
MinRR=GetMinRR(50)
MaxRR=1200
minLVET=(y[4]/x[4])*MinRR
maxLVET=(y[0]/x[0])*MaxRR

#x,y=resize(x,y,xmin=0.3, ymin=0.3) 
NewX,NewY=resize(x,y,xmin=MinRR,xmax=MaxRR,ymin=minLVET,ymax=maxLVET) 
print 'x is:  ',x 
print 'y is:  ',y
print 'NewX is:  ',NewX
print 'NewY is:  ',NewY

# p_guess is the starting estimate for the minimization
p_guess=(np.median(x),np.median(y),1.0,1.0) 
New_p_guess=(np.median(NewX),np.median(NewY),1.0,1.0) 

# Calls the leastsq() function, which calls the residuals function with an initial
# guess for the parameters and with the x and y vectors.  The full_output means that
# the function returns all optional outputs.  Note that the residuals function also
# calls the sigmoid function.  This will re开发者_JS百科turn the parameters p that minimize the
# least squares error of the sigmoid function with respect to the original x and y
# coordinate vectors that are sent to it.
p, cov, infodict, mesg, ier = scipy.optimize.leastsq( 
    residuals,p_guess,args=(x,y),full_output=1,warning=True)   

New_p, New_cov, New_infodict, New_mesg, New_ier = scipy.optimize.leastsq( 
    residuals,New_p_guess,args=(NewX,NewY),full_output=1,warning=True)   

# Define the optimal values for each element of p that were returned by the leastsq() function.
x0,y0,c,k=p 
print('''Reference data:\ 
x0 = {x0} 
y0 = {y0} 
c = {c} 
k = {k} 
'''.format(x0=x0,y0=y0,c=c,k=k)) 

New_x0,New_y0,New_c,New_k=New_p 
print('''New data:\ 
New_x0 = {New_x0} 
New_y0 = {New_y0} 
New_c = {New_c} 
New_k = {New_k} 
'''.format(New_x0=New_x0,New_y0=New_y0,New_c=New_c,New_k=New_k))

# Create a numpy array of x-values
xp = np.linspace(x.min(), x.max(), x.max()-x.min())
New_xp = np.linspace(NewX.min(), NewX.max(), NewX.max()-NewX.min())
# Return a vector pxp containing all the y values corresponding with the x-values in xp
pxp=sigmoid(p,xp)
New_pxp=sigmoid(New_p,New_xp)

# Plot the results 
plt.plot(x, y, '>', xp, pxp, 'g-')
plt.plot(NewX, NewY, '^',New_xp, New_pxp, 'r-')
plt.xlabel('x')
plt.ylabel('y',rotation='horizontal')
plt.grid(True)
plt.show()

Try this:

import numpy as np 
import matplotlib.pyplot as plt 
import scipy.optimize 

def GetMinRR(age):
    MaxHR = 208-(0.7*age)
    MinRR = (60/MaxHR)*1000
    return MinRR

def sigmoid(p,x):
    x0,y0,c,k=p 
    y = c / (1 + np.exp(-k*(x-x0))) + y0 
    return y 

def residuals(p,x,y): 
    return y - sigmoid(p,x) 

def resize(arr,lower=0.0,upper=1.0):
    # Create local copy
    result=arr.copy()
    # If the mins are greater than the maxs, then flip them.
    if lower>upper: lower,upper=upper,lower
    #----------------------------------------------------------------------------------------------    
    # The rest of the code below re-calculates all the values in x and then in y with these steps:
    #       1.) Subtract the actual minimum of the input x-vector from each value of x
    #       2.) Multiply each resulting value of x by the result of dividing the difference
    #           between the new xmin and xmax by the actual maximum of the input x-vector
    #       3.) Add the new minimum to each value of x
    #----------------------------------------------------------------------------------------------    
    # Subtracts right operand from the left operand and assigns the result to the left operand.
    # Note: c -= a is equivalent to c = c - a
    result -= result.min()

    # Multiplies right operand with the left operand and assigns the result to the left operand.
    # Note: c *= a is equivalent to c = c * a
    result *= (upper-lower)/result.max()

    # Adds right operand to the left operand and assigns the result to the left operand.
    # Note: c += a is equivalent to c = c + a
    result += lower
    return result


# Declare raw data for use in creating logistic regression equation
x = np.array([821,576,473,377,326],dtype='float') 
y = np.array([255,235,208,166,157],dtype='float') 

# Call resize() function to re-calculate coordinates that will be used for equation
MinRR=GetMinRR(50)
MaxRR=1200
# x[-1] returns the last value in x
minLVET=(y[-1]/x[-1])*MinRR
maxLVET=(y[0]/x[0])*MaxRR

print(MinRR, MaxRR)
#x,y=resize(x,y,xmin=0.3, ymin=0.3) 
NewX=resize(x,lower=MinRR,upper=MaxRR)
NewY=resize(y,lower=minLVET,upper=maxLVET) 
print 'x is:  ',x 
print 'y is:  ',y
print 'NewX is:  ',NewX
print 'NewY is:  ',NewY

# p_guess is the starting estimate for the minimization
p_guess=(np.median(x),np.min(y),np.max(y),0.01) 
New_p_guess=(np.median(NewX),np.min(NewY),np.max(NewY),0.01) 

# Calls the leastsq() function, which calls the residuals function with an initial
# guess for the parameters and with the x and y vectors.  The full_output means that
# the function returns all optional outputs.  Note that the residuals function also
# calls the sigmoid function.  This will return the parameters p that minimize the
# least squares error of the sigmoid function with respect to the original x and y
# coordinate vectors that are sent to it.
p, cov, infodict, mesg, ier = scipy.optimize.leastsq( 
    residuals,p_guess,args=(x,y),full_output=1,warning=True)   

New_p, New_cov, New_infodict, New_mesg, New_ier = scipy.optimize.leastsq( 
    residuals,New_p_guess,args=(NewX,NewY),full_output=1,warning=True)   

# Define the optimal values for each element of p that were returned by the leastsq() function.
x0,y0,c,k=p 
print('''Reference data:\ 
x0 = {x0} 
y0 = {y0} 
c = {c} 
k = {k} 
'''.format(x0=x0,y0=y0,c=c,k=k)) 

New_x0,New_y0,New_c,New_k=New_p 
print('''New data:\ 
New_x0 = {New_x0} 
New_y0 = {New_y0} 
New_c = {New_c} 
New_k = {New_k} 
'''.format(New_x0=New_x0,New_y0=New_y0,New_c=New_c,New_k=New_k))

# Create a numpy array of x-values
xp = np.linspace(x.min(), x.max(), x.max()-x.min())
New_xp = np.linspace(NewX.min(), NewX.max(), NewX.max()-NewX.min())
# Return a vector pxp containing all the y values corresponding with the x-values in xp
pxp=sigmoid(p,xp)
New_pxp=sigmoid(New_p,New_xp)

# Plot the results 
plt.plot(x, y, '>', xp, pxp, 'g-')
plt.plot(NewX, NewY, '^',New_xp, New_pxp, 'r-')
plt.xlabel('x')
plt.ylabel('y',rotation='horizontal')
plt.grid(True)
plt.show()

Your other related question is not closed, it appears you've registered twice, and stackoverflow is not letting you edit the other question because it doesn't recognize this user to be the same as this user.

Mainly all I've done in the code above is alter the New_p_guess. Finding the right values for the initial guess is kind of an art. If it could be done algorithmically, scipy wouldn't be asking you to do it. A little analysis can help, as well as "a feel" for your data. Knowing in advance what the solution should roughly look like, and therefore what values are reasonable within the context of the problem also helps. (Which is all just a long way of saying I guessed my way to choosing k=0.01.)