Decay curve best fit SciPy

I am having a problem when I try to find best fit to my data. Using scipy.optimize.curve_fit to create best fit. My data and code is:

EDIT You can download the data file from here. data is,

         a             b            b2
55478   1.07E+43    54395.93833 
56333   1.63E+43    54380.01385 
57540   2.57E+43    52393.31605 
61866   7.32E+43    52212.22838 52212.22838

code:

#!/usr/bin/env python
# -*- coding: utf-8 -*-


from __future__ import division

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import fit
import glob
import os
from scipy.optimize import curve_fit
import matplotlib.patches as patches

pf = pd.read_csv('/home/imhotep/Desktop/lala.csv', sep=',', encoding ='utf-8')



a1= pf['a'].max()
b1 = pf['b2'].max()
npoc=100

x = np.linspace((b1), (pf['b'].max()),npoc)
yy = np.linspace((pf['a'].min()), (pf['a'].max()), npoc)


fig = plt.figure()

ax4 = fig.add_subplot(111)

def h(x,k):
    return a1* (((x-(b1))/(k))**(-(5./3.)))


popt,pcov = curve_fit(h,x,yy)

print 'POPT,', popt,'PCOV',pcov
y_fi1 = h(x, *popt)

ax4.plot(x, y_fi1, label='fit', ls='-', color='blue')

ax4.plot(pf['b'], pf['a'], ls='None', color='blue', marker='o')

plt.show()

like that. When I run the code I'm getting that fit:

But, it should be roughly like that:

Can anyone tell me where I go wrong? I am beginner about curve fitting.

You want to fit a model to your 4 blue points described by a and b?

You should be doing fit in this direction then:

popt,pcov = curve_fit(h,b,a)

EDIT:

as mentioned in the comments of the question and this answer, you should be using the fit function only on your original data and then the newly created array using np.linspace to show the fit.

Here is what I got from your code:

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from __future__ import division

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

pf = pd.read_csv('lala.csv', sep=',', encoding ='utf-8')

a1 = pf['a'].max()
#b1 = pf['b2'].max()

x = pf["b"]
y = pf["a"]

def h(x,k,b1):
    return a1*((x-b1)/k)**(-5/3)

popt,pcov = curve_fit(h,x,y)

print 'POPT,', popt,'PCOV',pcov

xfit = np.linspace(x.min(),x.max(),100)
y_fi1 = h(xfit, *popt)

fig = plt.figure()
ax4 = fig.add_subplot(111)
ax4.plot(xfit, y_fi1, label='fit', ls='-', color='blue')
ax4.plot(x, y, ls='None', color='blue', marker='o')
plt.show()

Using the curve_fit to find only parameter k resulted in an error, therefore I included b1 as a search parameter. Then it does find a fit, however still not completelly satisfying. The output:

POPT, [   238.09666313  51973.04601693] 
PCOV [[ 21500.32886377 -22370.88448044] [-22370.88448044  23850.34961769]]