Fitting in log space and weighted least squares

[Beforerr]

Hi! First off, I really like many of the design choices in this package — great work!

I just have a couple of user questions:
1. How can I perform a fit in log space using this package?
2. If my data have uncertainties, is there a way to perform a weighted least-squares fit?

Thanks in advance for any guidance!

[fjebaker]

Hi @Beforerr thanks for your questions! I'll address them in order of ease:

• 2.) You can assign uncertainties in either X or Y. For example, if you have simple data (i.e. a single X value for every Y value), you can construct an InjectiveData dataset with, with variances in one or both axes. This will then use weighted least squares if using the LevenbergMarquadt optimizer, and otherwise uses a Gaussian or Poisson likelihood function that is sensitive to the variances.

  Also apologies I thought this was in the docs, but I see we only have InjectiveData without any variances at the moment.

• 1.) Probably the best way to do that is to transform your data into log space and define a custom transformer function that transforms the output of the model during fitting to the log space. If you can tell me a little more about what you're trying to do I can maybe put an example together :)

[Beforerr]

Hi @fjebaker . Thanks for the quick reply!

First, just a small question about the name InjectiveData — we don’t imply/require injective in the pure mathematical sense here, right?

Also, to clarify what I’m trying to do: I’m fitting data for a function like power law that spans many orders of magnitude, so it’s not appropriate to fit the data in linear scale since it performs poorly there. The best approach seems to be transforming the data into log scale and fitting in log space to optimize the fractional deviation.

[fjebaker]

First, just a small question about the name InjectiveData — we don’t imply/require injective in the pure mathematical sense here, right?

Yes, that's right. As they say, the two hardest problems in computer science are naming things and cache invalidations. At some point in the past I settled on InjectiveData, but it isn't really a very good descriptor for it.

"Injective" in this case is solely to mean there is a one X to one Y correspondence (they need not be unique pairs), as opposed to binned data, where there may be a bin low and bin high in X for each Y (I welcome better names if you have a suggestion!).

Regarding logs:

I've never been completely sure of what happens to the likelihood functions when you log all your data. I'm sure there's a good derivation somewhere, but I haven't really had the time to look for it.

In SpectralFitting.jl you're probably best off applying whatever transformation you want to your data before wrapping it into one of the dataset types. You can then either fit a straight line model to get the powerlaw, or create a model wrapper that applies the log function to whatever model you have. That doesn't really avoid the huge dynamical range problems in floating point representation in the models, but might help you in the fit?

The latter is something I have had on the TODO list for a while (apply arbitrary functions as model wrappers), so if it's something you would be interested in trying, I can put that together tomorrow.

[Beforerr]

Regarding the name, I’d suggest PairedData or PointwiseData — both capture the idea of a one-to-one correspondence between X and Y values more clearly.

As for the likelihood functions, it depends on the noise assumption. When the flux spans several decades in energy and we do not know the uncertainty, it’s often reasonable to assume roughly constant fractional uncertainty. In that case, a log-normal noise model applies:

$$ \ln F_i \sim \mathcal{N}\left(\ln F_{\text {model }}\left(E_i ; \theta\right), \sigma^2\right) $$

which is equivalent to fitting in log space by minimizing squared residuals in $ln 𝐹$.

I previously thought Log10Flux does something like that, but after checking it seems that it serves a different purpose.

[fjebaker]

I like PointwiseData, I may adopt that in the next release and deprecate InjectiveData. I agree it captures the intention better.

If you assume constant fractional uncertainty (or have log-normal data), then yes it is okay to assume the data in log space is normally distributed. To do this in SpectralFitting.jl then you would do the equivalent of e.g. InjectiveData(log10.(x), log10.(y)) and define your models similarly -- for which it would be most convenient if I add the function wrapper, which I will try to address tomorrow.

The Log10Flux is an XSPEC model that measures the decades of flux of a given model. It does not, as you say, do any transformations.