What distribution should I fit to this data I am trying to find a distribution that fits my data (3500+ data points) with satisfying goodness of fit (gof), I use the Kolmogorov-Smirnov test and its p-value as a gof measurement (p-value > 0.1).
I have tried the plfit.m and plpva.m program from Clauset et al. to fit a power-law distribution to my data, but the p-value is close to zero indicating it is not a good fit.
The log-normal and exponential distribution is also tested using the R package poweRlaw, but I still can't get a good enough p-value (> 0.1). However, I think the fitted curve is just close enough to the empirical data, as the picture shows (generated by the poweRlaw package).
I am completely new to this kind of job but I have to report the gof in my paper, so I am wondering:


*

*Did I do the fitting job in the correct way? (I didn't modify the Matlab or R programs)

*Is Kolmogorov-Smirnov test a proper approach to measure the gof? 

*If I remove the extreme cases (2 data points on the right tail), will it become better?

*What distribution should I fit my data?


Thanks!

 A: I've plotted log(data) vs abscissa and dropped two last points, graph below

data going below x=500 would be good fitted with just a linear function - which means exponential distribution because we plotted log(data).
But from 0 to 500 things looks different, maybe power law? Thus, I don't know if data could be fitted with one distribution function.
UPDATE
I played with truncated power law, and indeed this might work
Code (Python3.7, Anaconda x64 Win 10)
import numpy as np
from scipy.special import gamma
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

xs = []
ys = []

with open('my_data.csv') as f:
     for line in f:
         x, y = line.split(',')
         xs.append(float(x))
         ys.append(float(y))

xs = np.asarray(xs)
ys = np.asarray(np.log(ys))

x0 = xs[0]

def F(x, a, b):
    return np.power((x-x0)+1.0, -a) * np.exp(-b*(x-x0))

def logF(x, a, b):
    return -a*np.log((x-x0)+1.0) + (-b*(x-x0))

popt, pcov = curve_fit(logF, xs, ys)
print(popt)

plt.plot(xs, ys, 'b*', label='data')
plt.plot(xs, logF(xs, *popt), 'g-', label='fit')
plt.show()

produced reasonable fit graph. And if you print sqrt of pcov diagonal, errors looks small
print(popt)
print(np.sqrt(np.diag(pcov)))

produced
[0.38962489 0.00140291]
[2.51804031e-03 2.10496530e-05]


