Sadly, there is no general solution to this problem.
You need to commit to a model, and tune the model parameters to suit the known points.
If there is no noise involved and the points are not shaped in an easily recognizable shape then I would use the lowest polynomial fit that passes through all the points and use the integral of that as an approximation of the wanted integral.
Your approach is also sensible, though it will only work reasonably if there is a high volume of points.
If the approximation is important then you can also try writing a program which generates fits using different combinations of typical functions (in the Wolfram Language this is somewhat easy to achieve) and pick the function which fits with the lowest complexity (where complexity can be measured in number and weirdness of terms or similar. Again, in Wolfram Language there are facilities for this).