# Data stitching
**URL:** https://discourse.julialang.org/t/data-stitching/89662
**Category:** Probabilistic Programming
**Tags:** turing
**Created:** [November 2, 2022, 2:21pm UTC](https://discourse.julialang.org/t/data-stitching/89662 "2022-11-02T14:21:00Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@gustaphe](https://discourse.julialang.org/u/gustaphe)
#### Post date: [November 2, 2022, 2:21pm UTC](https://discourse.julialang.org/t/data-stitching/89662/1 "2022-11-02T14:21:00Z")
I’m looking for a package or an algorithm, alternatively an explanation of why what I want is impossible or bad:
I have a set of data `x1, y1 ≈ f.(x1)` and one `x2, y2 ≈ a*f.(x2) .+ b` for some natural black box function `f` and constants `a, b`. I would like to merge these two datasets by scaling and shifting `y2`, essentially finding `a` and `b`. The physical motivation is, for instance, spectrometry of the same features using different integration times.
If enough of `x2` has “close neighbors” in `x1`, this is easy to do decently manually, but I find surprisingly little information on automating the process. My current approach is to minimize wrt `p` the sum over each possible pair of points `(x1[i], y1[i]), (x2[j], y2[j])` of the error function
```julia
(p[1]*y2+p[2] - y1)^2/(1 + (x2 - x1)^2/w^2)
```
where the points are weighted by the lorentzian `1/(1 + Δx^2/w^2)` of width `w` chosen as some multiple of `mean(diff(x))`. Results are not great so far.
Has anyone here seen anything similar to this? Image stitching seems adjacent, but not a perfect match.
---
### Author:  [@gustaphe](https://discourse.julialang.org/u/gustaphe)
#### Post date: [November 2, 2022, 2:39pm UTC](https://discourse.julialang.org/t/data-stitching/89662/2 "2022-11-02T14:39:02Z")
MWE:
```julia
using Plots, Optim
f(x) = exp(-x^2)
x1 = randn(50) .- 1.0
x2 = randn(30) .+ 1.0
y1 = f.(x1)
y2 = 1.5*f.(x2) .+ 0.5
p1 = plot([x1, x2], [y1, y2]; seriestype=:scatter) # don't line up
w = 0.5*sum(diff(sort(vcat(x1, x2))))/(length(x1)+length(x2))
error_fun(p) = sum(
(p[1]*y2 .+ p[2] - y1)^2/(1 + (x2 - x1)^2/w^2) for (x1, y1) in zip(x1, y1), (x2, y2) in zip(x2, y2)
)
result = optimize(error_fun, [1.0, 0.0]).minimizer
p2 = plot(vcat(x1, x2), vcat(y1, result[1]*y2 .+ result[2]); seriestype=:scatter) # should line up
plot!(p2, f; seriestype=:path)
plot(p1, p2)
```
---
### Author:  [@dlakelan](https://discourse.julialang.org/u/dlakelan)
#### Post date: [November 2, 2022, 3:33pm UTC](https://discourse.julialang.org/t/data-stitching/89662/4 "2022-11-02T15:33:11Z")
My usual thing when faced with unknowns is to formulate a Bayesian inference problem.
I’m going to assume that you don’t have an explicit form for the `f()` function. In that case, at least what you have is some noisy measures of `f()` in your plot maybe that’s the blue data. And you have some noisy measures of `a*f()+b` the red data.
I would choose some small ranges of data, and say for each range that a line fits through the data… (let’s call this a set of of local taylor series of f)
```julia
blue ~ Normal((m*(x-xcenter) + c),scale)
red ~ Normal(a * (m*(x-xcenter) + c) + b, scale*a)
```
where there’s a different m and c for each small range and xcenter is the center of each range. You probably don’t need to deal with all the ranges, but several ranges where m varies is probably preferred because you will get better inference for a I think.
The a and b would be shared across all the ranges, and the m and c would be different for each range.
Sample from this posterior distribution using Turing, and then you’ll see not only inference for a and b but also their uncertainty, and the correlation between them.
---
### Author:  [@Eben60](https://discourse.julialang.org/u/Eben60)
#### Post date: [November 2, 2022, 7:55pm UTC](https://discourse.julialang.org/t/data-stitching/89662/5 "2022-11-02T19:55:31Z")
I can think here of two approaches. One - maybe the data can be described by some relatively simple function with a few (3…6) parameter - a polynom, rational function, gauss, some combination of these, whatever (“With four parameters I can fit an elephant, and with five I can make him wiggle his trunk.” 🙂). It is better if there were some physical rational behind, but not necessary. Then, first fit the first dataset, then add the `a` and `b` to the parameters and fit the combined dataset.
Another approach - local one. Interpolate both data separately onto a regular grid, then you have your “new” data all with the same x-values, then just do LSQ fit for `a` and `b` .
---
### Author:  [@SteffenPL](https://discourse.julialang.org/u/SteffenPL)
#### Post date: [November 2, 2022, 8:53pm UTC](https://discourse.julialang.org/t/data-stitching/89662/6 "2022-11-02T20:53:48Z")
Here is my attempt (please note that statistics is something I would really like to be good at; but I’m not…)
If the sample points are dense enough, a simple interpolation could be sufficient to approximate `f.(x2)`. Once we have an approximation of `f.(x2)` the rest is linear regression.
```julia
using CurveFit, Interpolations
f(x) = exp(-x^2)
x1 = randn(50) .- 1.0
x2 = randn(30) .+ 1.0
y1 = f.(x1)
y2 = 1.5*f.(x2) .+ 0.5
p = sortperm(x1)
permute!(x1, p)
permute!(y1, p)
f_fit = linear_interpolation(x1, y1)
# restrict to range of x1 (otherwise one would need a better extrapolation model)
ind = minimum(x1) .< x2 .< maximum(x1)
b, a = linear_fit(f_fit( x2[ind] ), y2[ind])
# a = 1.5092573052474938
# b = 0.4960624235320712
```
---
### Author:  [@dlakelan](https://discourse.julialang.org/u/dlakelan)
#### Post date: [November 2, 2022, 9:17pm UTC](https://discourse.julialang.org/t/data-stitching/89662/7 "2022-11-02T21:17:40Z")
> [@SteffenPL](#):
>
> If the sample points are dense enough, a simple interpolation could be sufficient to approximate `f`
Not if there is significant noise in the measurements, interpolation will be fitting a lot of noise
---
### Author:  [@Dan](https://discourse.julialang.org/u/Dan)
#### Post date: [November 2, 2022, 10:48pm UTC](https://discourse.julialang.org/t/data-stitching/89662/8 "2022-11-02T22:48:47Z")
@SteffenPL basically did what I wanted. So this is just another version, incorporating intuition from original and other posts that close points are key. Luckily we have to estimate just two parameters, and we may get over some noise.
Now the code:
```julia
using IterTools
f(x) = exp(-x^2)
x1 = randn(50) .- 1.0
x2 = randn(30) .+ 1.0
y1 = f.(x1)
y2 = 1.5*f.(x2) .+ 0.5
function findclosepairs(x1, x2, thresh)
res = Tuple{Float64,Int,Int}[]
for i in eachindex(x1), j in eachindex(x2)
abs(x1[i]-x2[j])abs(xs[i][1]-xs[i-1][1]), filter(i->xs[i][2] != xs[i-1][2], 2:length(xs))))
return xdeltas[clamp(length(xdeltas)÷5, 1, length(xdeltas))]
end
function findaffine(y1, y2, cp)
xx = y1[nth.(cp, 2)]
yy = y2[nth.(cp, 3)]
M = [ones(size(xx)) xx]
v = M \ yy
end
```
After this prep, we can:
```julia
julia> cp = findclosepairs(x1, x2, threshsuggest(x1, x2))
4-element Vector{Tuple{Float64, Int64, Int64}}:
(0.0018581113413391748, 46, 13)
(0.0034200499950651775, 18, 2)
(0.005578601033654262, 47, 15)
(0.006237570357774835, 17, 2)
julia> findaffine(y1,y2,cp)
2-element Vector{Float64}:
0.507535289631714
1.4827932614125852
```
The actual estimation is more accurate with just 3 pairs.
But I guess this code is a good part of the automation OP discussed (and some more optimizations possible for sure).
---
### Author:  [@dlakelan](https://discourse.julialang.org/u/dlakelan)
#### Post date: [November 3, 2022, 4:18am UTC](https://discourse.julialang.org/t/data-stitching/89662/9 "2022-11-03T04:18:44Z")
I’m working up a Turing example based on your MWE and including it as one of my Vignettes for a series on “how to use Julia for data analysis”
---
### Author:  [@dlakelan](https://discourse.julialang.org/u/dlakelan)
#### Post date: [November 3, 2022, 4:46am UTC](https://discourse.julialang.org/t/data-stitching/89662/10 "2022-11-03T04:46:09Z")
Ok, here’s my example. The Turing sampling took 2-3 seconds.
```julia
using StatsPlots, Turing, LinearAlgebra,Random
Random.set_global_seed!(1)
f(x) = exp(-x^2)
x1 = randn(50) .- 1.0
x2 = randn(30) .+ 1.0
y1 = f.(x1)
y2 = 1.5*f.(x2) .+ 0.5 ## true values of a = 1.5 and b = 0.5
p1 = plot([x1, x2], [y1, y2]; seriestype=:scatter) # don't line up
display(p1)
## how to determine a,b if we don't know f()?
@model function stitchdata(xbase,ybase,xtform,ytform)
a ~ Gamma(5,1.0/4) # a is of order 1
b ~ Gamma(5,1.0/4) # b is of order 1
m ~ MvNormal(repeat([0.0],4),100.0^2*I(4))
c ~ MvNormal(repeat([0.0],4),100.0^2*I(4))
xbneighs = Vector{eltype(xbase)}[]
ybneighs = Vector{eltype(ybase)}[]
xtneighs = Vector{eltype(xbase)}[]
ytneighs = Vector{eltype(ybase)}[]
centers = (-1.0,-0.5,.5,1.0)
for (i,dx) in enumerate(centers)
f = x -> x > dx - .25 && x < dx + 0.25
xneighindex = findall(f, xbase)
#@show xneighindex
push!(xbneighs, xbase[xneighindex])
push!(ybneighs, ybase[xneighindex])
tneighindex = findall(f,xtform)
push!(xtneighs, xtform[tneighindex])
push!(ytneighs, ytform[tneighindex])
end
for (i,(x,y)) in enumerate(zip(xbneighs,ybneighs))
Turing.@addlogprob!(logpdf(MvNormal(m[i].*(x .- centers[i]) .+ c[i], 0.05^2*I(length(y))), y))
end
for (i,(x,y)) in enumerate(zip(xtneighs,ytneighs))
Turing.@addlogprob!(logpdf(MvNormal(a.*(m[i].*(x .- centers[i]) .+ c[i]) .+ b, 0.05^2*I(length(y))), y))
end
end
mm = stitchdata(x1,y1,x2,y2)
s = sample(mm,NUTS(400,0.8),200)
plot(s[:,[:a,:b],1])
```
This produces the data:
And the inference for a and b are:
I just ad-hoc decided on the “neighborhoods” to use.
I had to use the `@addlogprob!` macro because I’m subsetting the data and Turing is confused about that if I don’t.
Also, it might make more sense to actually have some measurement noise and/or fit the scale of the measurement noise…
when I added Normal(0.0,0.05) noise to the data, I get:
and
---
### Author:  [@gustaphe](https://discourse.julialang.org/u/gustaphe)
#### Post date: [November 3, 2022, 9:28am UTC](https://discourse.julialang.org/t/data-stitching/89662/11 "2022-11-03T09:28:48Z")
Nice! Guess I’m learning Turing… As an addition, one could run a final
```julia
y3 = (y2 .- mean(s[:b]))/mean(s[:a])
plot([x1, x2], [y1, y3]; seriestype=:scatter)
plot!(f; seriestype=:path)
```
to get it back in the domain of data stitching.
---
### Author:  [@Dan](https://discourse.julialang.org/u/Dan)
#### Post date: [November 3, 2022, 11:23am UTC](https://discourse.julialang.org/t/data-stitching/89662/12 "2022-11-03T11:23:38Z")
`Turing` certainly gives the prettiest solution. Here is another Turing model for this problem. I’ve tried to simplify a bit the model for people getting anxious from `@addlogprob!`:
```julia
using StatsPlots, Turing, LinearAlgebra,Random, IterTools
Random.set_global_seed!(1)
# prep example data
f(x) = exp(-x^2)
x1 = randn(50) .- 1.0
x2 = randn(30) .+ 1.0
y1 = f.(x1)
y2 = 1.5*f.(x2) .+ 0.5 ## true values of a = 1.5 and b = 0.5
# helper functions for close pairs
function findclosepairs(x1, x2, thresh)
res = Tuple{Float64,Int,Int}[]
for i in eachindex(x1), j in eachindex(x2)
abs(x1[i]-x2[j])<=thresh && push!(res, (abs(x1[i]-x2[j]),i,j))
end
sort!(res)
end
function threshsuggest(x1,x2)
xs = vcat([(x,0) for x in x1], [(x,1) for x in x2])
sort!(xs)
xdeltas = sort(map(i->abs(xs[i][1]-xs[i-1][1]), filter(i->xs[i][2] != xs[i-1][2], 2:length(xs))))
return xdeltas[clamp(length(xdeltas)÷5, 1, length(xdeltas))]
end
# find close pairs for Bayesian model
cp = findclosepairs(x1, x2, 4*threshsuggest(x1, x2))
xbase = x1[nth.(cp, 2)]
ybase = y1[nth.(cp, 2)]
xtform = x2[nth.(cp, 3)]
ytform = y2[nth.(cp, 3)]
# Bayesian model or "how likely is the data given some parameters"
@model function bayes_stitch(xbase,ybase,xtform,ytform)
a ~ Gamma(5,1.0/4) # a is of order 1
b ~ Gamma(5,1.0/4) # b is of order 1
σ ~ InverseGamma(1.0,1.0)
# using Exponential(0.1) also good, especially for wackier f(.)
L = length(xbase)
for i=1:L
ytform[i] ~ Normal(a*ybase[i]+b, abs(xtform[i]-xbase[i])*σ)
# instead of Normal this could be Cauchy for noisier data
end
end
# create and sample model from closepoint pairs data
bayesmm = bayes_stitch(xbase, ybase, xtform, ytform);
s = sample(bayesmm,NUTS(400,0.8),200)
plot(s[:,[:a,:b,:σ],1])
```
This gives:
Notes:
- The resulting variance posterior distribution is not too large, meaning the rest of the data fits the model decently (otherwise model would claim variance is high and everything is noise).
- The goal of this version is to remove some manual parameters from other model and simplify the model description.
- Having said that, this model loses a bit of power that @dlakelan 's model gets from large clusters (more than 2 points).
- There are still a few parameters which should depend on nature of samples and especially nature of `f(.)`:
- choosing close points threshold (or list size)
- choosing prior distributions for variance and error
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [November 3, 2022, 11:55am UTC](https://discourse.julialang.org/t/data-stitching/89662/13 "2022-11-03T11:55:06Z")
A beginner’s question for probability experts:
For the above models that are free of noise and using a trivial case of dense overlapping input data (`x1 = -4:0.1:4; x2 = -4:0.1:4`), why are the Turing solutions (highest probability peak, mean or median) always shifted from the true solution (`a = 1.5; b = 0.5`), as if there is a bias…
Any reference or link would be great.
---
### Author:  [@Dan](https://discourse.julialang.org/u/Dan)
#### Post date: [November 3, 2022, 12:15pm UTC](https://discourse.julialang.org/t/data-stitching/89662/14 "2022-11-03T12:15:16Z")
> [@rafael.guerra](#):
>
> For the above models that are free of noise and using a trivial case of dense overlapping input data (`x1 = -4:0.1:4; x2 = -4:0.1:4`),
Did you try this? When there are 2 `x` points in both `x1` and `x2`, the distributions should very concentrated (and your ranges have tons of `x`’s in common).
For example, in the model I’ve just ran, there is a varince of `abs(x1-x2)*sigma` which would be zero in this case (forcing the approx equality to be equality).
This comment applies to the model I’ve shown. dlakelan’s model is a bit more lineant.
The priors introduce a bit of bias (which decreases as more data is given).
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [November 3, 2022, 12:23pm UTC](https://discourse.julialang.org/t/data-stitching/89662/15 "2022-11-03T12:23:18Z")
Plugging-in such x1 and x2 in the first dlakelan’s model, the results plot like:
---
### Author:  [@Dan](https://discourse.julialang.org/u/Dan)
#### Post date: [November 3, 2022, 12:27pm UTC](https://discourse.julialang.org/t/data-stitching/89662/16 "2022-11-03T12:27:57Z")
Because essentially the sampler is choosing a best-fit line, there is some ambiguity as to the slope, and when you under-estimate the slope (`a`), you need to compensate with over-estimating the intercept (`b`).
Thus the uncertainty is in an ellipse-ish shape in `a`-`b` space but not an axis-parallel ellipse, and the distributions we are getting are the marginals of this rotated ellipse-ish.
Non-axis ellipses projection are still centered, but in this case, there is a non-linearity in `f(.)` which introduces heteroscedasticity (changing variance of `b` as function of `a`), and this distorts the ellipses and produces this bias.
This can be tested using a linear `f` such as `f(x) = kx+h`.
Another issue is the convergence quality / chain correlation and length of the Bayesian sampler chains.
Perhaps someone will give a good link in a reply.
UPDATE:
As suggested, I’ve tested dlakelan model with -4:0.1:4 and linear function f(x). The resulting distributions:
---
### Author:  [@dlakelan](https://discourse.julialang.org/u/dlakelan)
#### Post date: [November 3, 2022, 2:08pm UTC](https://discourse.julialang.org/t/data-stitching/89662/17 "2022-11-03T14:08:51Z")
> [@Dan](#):
>
> The priors introduce a bit of bias (which decreases as more data is given).
Indeed the priors basically say that both a and b are most likely to be 1.0 so it is pulling the slope down a little and the offset up a little. It would be possible to make a prior which was a little flatter through this region if desired.
Also the line is fitting over a region of width 0.5 but over this region there is some mild nonlinearity. It would be possible to introduce an extra Taylor series term and handle this nonlinearity a little better.
---
### Author:  [@dlakelan](https://discourse.julialang.org/u/dlakelan)
#### Post date: [November 3, 2022, 7:34pm UTC](https://discourse.julialang.org/t/data-stitching/89662/18 "2022-11-03T19:34:52Z")
> [@gustaphe](#):
>
> Nice! Guess I’m learning Turing…
I’ve got a plan to produce a bunch of vignettes about doing common tasks in Julia and put them on YouTube. In addition to that I have some plans for larger multi-part projects. One of my plans is to talk about Bayesian inference in general, and Turing in particular. The thing is that my shortest coverage of these ideas was at least 1 hr long, which I think is too much to plop into my YouTube stream before people have come up to speed.
in any case, I hope you’ll tune in when I get around to doing the Bayesian + Turing version. I’ll post the YouTube playlist probably next week. Right now I have only longer form content and I decided the vignettes are really needed.
---
### Author:  [@bertschi](https://discourse.julialang.org/u/bertschi)
#### Post date: [November 3, 2022, 8:26pm UTC](https://discourse.julialang.org/t/data-stitching/89662/19 "2022-11-03T20:26:37Z")
This problem appears to be a nice fit for Gaussian processes, i.e., using the following generative model:
\begin{align\*} f &\sim GP(\cdot) \\ y\_1 &\sim \mathcal{N}(f(x\_1), \sigma) \\ y\_2 &\sim \mathcal{N}(a \cdot f(x\_1) + b, \sigma) \end{align\*}
Thus, the Gaussian process describes the latent/unobserved function f and depending on the kernel takes automatically care of smoothing, i.e., drawing information from neighbouring data points. I’m not quite sure which kernel would correspond to your Lorentzian weighting though. Where did you get that formula from? Just a modelling assumptions or is there a more fundamental reason for choosing 1 / \left(1 + \left(\frac{\Delta x}{w}\right)^2\right) to weight neighbouring data points?
---
### Author:  [@gustaphe](https://discourse.julialang.org/u/gustaphe)
#### Post date: [November 3, 2022, 8:35pm UTC](https://discourse.julialang.org/t/data-stitching/89662/20 "2022-11-03T20:35:20Z")
Just intuition. Didn’t want stuff far away to factor in, lorentzians are fairly simple. Could just as well have been a gaussian or a tent.
---
### Author:  [@bertschi](https://discourse.julialang.org/u/bertschi)
#### Post date: [November 3, 2022, 8:38pm UTC](https://discourse.julialang.org/t/data-stitching/89662/21 "2022-11-03T20:38:26Z")
Ok, then you might want to try Gaussian processes which allow you to compare different smoothness assumptions on f as well as learn the effective neighbourhood range via kernel parameters.
When I find some time, I might try to put together a MWE Turing example …
[Next page](https://discourse.julialang.org/t/data-stitching/89662.md?page=2)