I have developed a theoretical distribution in the form of
$$ f(x) = \frac{\beta}{\alpha}\left(1+\frac{x}{\alpha}\right)^{-\beta - 1} $$
Where $\alpha$ and $\beta$ are parameters of the model with physical meanings. I have read that we can't simply estimate parameters from KS test and then evaluate its goodness of it. How do I go about finding the "best" parameters and then testing the goodness of fit?
What you have developed is a Lomax distribution.
Due to the heavy right tail, the (more-or-less) default approach of maximum likelihood doesn't work so well in small to medium samples when the sample coefficient of variation is less than $1$ (the population coefficient of variation, when it exists, is always $>1$.) (An analysis of the maximum likelihood estimates for the Lomax distribution.) If $\beta \leq 2$, you are also out of luck for the "usual" method of moments estimator, as the second moment (the variance) doesn't exist. Of course, we don't know that a priori, but if our estimate of $\beta$ is close to $2$, we should be concerned.
The method-of-moments estimators use the sample mean and variance and "reverse-engineer" the parameter values that would make the population mean and variance equal to the sample mean and variance. The relevant formulae can be found on the Wikipedia page linked to above:
$$\begin{eqnarray} \mu &=& {\alpha \over \beta - 1} \\ \sigma^2 &=& {\alpha^2\beta \over (\beta-1)^2(\beta-2)} = \mu^2{\beta \over \beta-2} \end{eqnarray}$$
A minor bit of algebra, substituting the sample mean $\bar{x}$ and sample variance $s^2$ in the above, gets us to:
$$\begin{eqnarray} \hat{\beta} &=& {2s^2 \over s^2 - \bar{x}^2} \\ \hat{\alpha} &=& \bar{x}(\hat{\beta}-1) \end{eqnarray}$$
Note that this will fail to provide meaningful estimates when the sample coefficient of variation ${s \over \bar{x}} \leq 1$.
The maximum likelihood approach does not admit a closed-form solution. Instead, we have the following optimization problem:
$$\max_{\alpha, \beta} \;\;n\ln\beta - n\ln\alpha-(\beta+1)\sum_{i=1}^n\ln\left(1 + {x_i\over \alpha}\right)$$
This has to be solved numerically, and any number of multivariate optimization algorithms will do so.
