That's an okay toy project, but you shouldn't expect too much from it. Dozens (or hundreds!) of functionals have been fit to experimental data over several decades with varying success. You didn't specify what property you would be looking at; however, one of the key problems in the fitting effort is that depending on the property, the experiments may be quite far from the calculation.
Single-point energies are the easiest target, since they are comeasurable between wave function theory and density functional theory. But, this problem has already been solved, starting out with Truhlar's functionals and finishing with the tour de force by Mardirossian and Head-Gordon, in which billions of functional forms were fit to a huge dataset of high-level (i.e. coupled-cluster) data, and the most predictive functional forms were chosen; even their most sophisticated functionals contain few parameters. (They also do a good job of testing the functionals with molecules that don't appear in their training set.) You should definitely have a look at their work on the wB97X-V Phys. Chem. Chem. Phys. 16, 9904 (2014), B97M-V J. Chem. Phys. 142, 074111 (2015), and wB97M-V
J. Chem. Phys. 144, 214110 (2016) functionals, as well as the double hybrid wB97M(2) functional J. Chem. Phys. 148, 241736 (2018).
Anyway, the fitting procedure itself is quite simple. What you do is set up a functional form and make an initial guess for the parameters. You converge the wave functions with this functional, and calculate the energy and its derivatives with respect to the parameters in your functional for all of your molecules (this will probably require some heavy modifications to the quantum chemistry code). Next, you update the parameters by requiring that e.g. the root-mean-squared error is minimized for your molecular test set; if your parameters are linear then you get a matrix equation, if you have non-linear parameters you may end up having to do direct optimization with line searches and so on. Then, you recompute the wave functions, the energy and its derivatives for all your molecules, and re-fit. Once the parameters have converged, you're done.
The real problem is ensuring that what you're doing actually makes sense. You have to have a large training set of molecules, ensure that the functional is numerically stable i.e. doesn't exhibit pathological grid dependence (e.g. SCAN is terrible in this respect even though it's "from first principles"), and make sure that the accuracy is transferable. These are all things that have been screwed up in one way or another in many recent functionals, except the aforementioned ones.