0

我正在使用 5 个模型，用于 5 个数据集，它们之间共享一些参数以进行全局拟合。我不确定我是否真的在这里进行了全局调整。我正在尝试为剂量反应类型的数据实施 Schild / Global fit。 

def gauss(params,xdata,ydata):
        ## share parameters: hillSlope,schildSlope,top,bottom,pA2
        hillSlope = params['hillSlope'].value
        schildSlope = params['SchildSlope'].value
        top = params['top'].value
        bottom = params['bottom'].value
        pA2 = params['pA2'].value


        EC50_1 = params['ec50_1'].value
        B_1 = params['B_1'].value

        Antag_1 = 1+(B_1/(10**(-1*pA2)))**schildSlope
        LogEC_1=np.log10(EC50_1*Antag_1)
        y_model_1 = bottom + (top-bottom)/(1+10**((LogEC_1-xdata)*hillSlope))
        resid1 = ydata[0] - y_model_1

        #print ydata[0]



        EC50_2 = params['ec50_2'].value
        B_2 = params['B_2'].value
        Antag_2 = 1+(B_2/(10**(-1*pA2)))**schildSlope
        LogEC_2 = np.log10(EC50_2*Antag_2)
        y_model_2 = bottom + (top-bottom)/(1+10**((LogEC_2-xdata)*hillSlope))
        resid2 = ydata[1] - y_model_2


        EC50_3 = params['ec50_3'].value
        B_3 = params['B_3'].value
        Antag_3 = 1+(B_3/(10**(-1*pA2)))**schildSlope
        LogEC_3=np.log10(EC50_3*Antag_3)
        y_model_3 = bottom + (top-bottom)/(1+10**((LogEC_3-xdata)*hillSlope))
        resid3 = ydata[2] - y_model_3


        EC50_4 = params['ec50_4'].value
        B_4 = params['B_4'].value
        Antag_4 = 1+(B_4/(10**(-1*pA2)))**schildSlope
        LogEC_4=np.log10(EC50_4*Antag_4)
        y_model_4 = bottom + (top-bottom)/(1+10**((LogEC_4-xdata)*hillSlope))
        resid4 = ydata[3] - y_model_4


        EC50_5 = params['ec50_5'].value
        B_5 = params['B_5'].value
        Antag_5 = 1+(B_5/(10**(-1*pA2)))**schildSlope
        LogEC_5=np.log10(EC50_5*Antag_5)
        y_model_5 = bottom + (top-bottom)/(1+10**((LogEC_5-xdata)*hillSlope))
        resid5 = ydata[4] - y_model_5

        return np.concatenate((resid1,resid2,resid3,resid4,resid5))




    params = Parameters()
    B = [1.25e-7,2.5e-8,1e-8,5e-9,0]

    #data
    dataComplete = [[]]*5
    dataComplete[0]=conc1
    dataComplete[1]=conc2
    dataComplete[2]=conc3
    dataComplete[3]=conc4
    dataComplete[4]=conc5

    ydata = np.array(dataComplete)
    xdata = np.array(dose)

    ## parameter top
    topValues = []
    bottomValues = []
    for i in range(5):
        topValues.append(max(dataComplete[i]))
        bottomValues.append(min(dataComplete[i]))

    print "bottomValues", bottomValues
    print "topValues", topValues

    ec50Values = []

     #calculate ec50 values using some data
    def calculateEC50Values(data):
        Yvalue_50 = (max(data)-min(data)) /2
        Distance_To_Central_Y_Value = np.abs(data - Yvalue_50)
        LocationOfNearest = np.argmin(Distance_To_Central_Y_Value)
        Xvalue_50 = xdata[LocationOfNearest]
        return 10**Xvalue_50


    for i in range(5):
        ec50Values.append(calculateEC50Values(np.array(dataComplete[i])))



    ## Add parameters
    ## Fixed
    params.add('bottom',value=-0.67,vary=True)
    params.add('top',value=116.96,vary=True)
    params.add('hillSlope',value=1.,vary=True)
    params.add('SchildSlope', value=1.,vary=True)
    params.add('pA2',value=7., vary=True)

    ## concentration values
    params.add('B_1',value=B[0], min=0, vary=False)
    params.add('B_2',value=B[1], min=0, vary=False)
    params.add('B_3',value=B[2], min=0, vary=False)
    params.add('B_4',value=B[3], min=0, vary=False)
    params.add('B_5',value=B[4], min=0, vary=False)

    ## estimated ec50 values
    params.add('ec50_1',value=ec50Values[0], min=0, vary=True)
    params.add('ec50_2',value=ec50Values[1], min=0, vary=True)
    params.add('ec50_3',value=ec50Values[2], min=0, vary=True)
    params.add('ec50_4',value=ec50Values[3], min=0, vary=True)
    params.add('ec50_5',value=ec50Values[4], min=0, vary=True)


    ## curve for 0M, agonist only.
    result = minimize(gauss,params,args=(xdata,dataComplete))
    #final = dataComplete[dataIndex] + result.residual[-22:]
    lmfit.printfuncs.report_fit(result.params)

    print result.residual
4

0 回答 0