Merging piecewise polynomial (pp) structures

Is there any method how to merge two or more "pp" (peicewise polynomial) structures in a case of consecutive "pp" structures with some range overlap to get one "pp" structure describing merged function?

Answer 1

[Bruno Luong]

Using BSFK from this FEX https://www.mathworks.com/matlabcentral/answers/1810010-free-knot-spline-approximation-bsfk-problem?s_tid=ans_lp_feed_leaf

% Example as in John's answer
breaks = [4 10]; coefs = [ 0 0 1 -2 53];
PP1 = mkpp(breaks,coefs);
breaks = [8 15]; coefs = [ -1 6 1 4 77];
PP2 = mkpp(breaks,coefs);

% overlap interval
oI = [8 10];
li = diff(oI); % its length
f = @(x) ppval(PP1,x);
g = @(x) ppval(PP2,x);
% define smooth transition weight function
sigfun = @(x) x.^2.*(3-2*x);
w2 = @(x) sigfun((x-oI(1))/li);
w1 = @(x) 1-w2(x);
% weighted sum of f and g, applied only in oI
h = @(x) w1(x).*f(x) + w2(x).*g(x);

%% Fit on 3 intervals
arg = {5,[],[],struct('KnotRemoval', 'yes','Animation',0)}; % EDIT first argument 4 => 5

xo = linspace(oI(1),oI(2),100);
ppo = BSFK(xo,h(xo),arg{:});

x1 = linspace(PP1.breaks(1),oI(1));
ppf = BSFK(x1,f(x1),arg{:}); % we can also use PP1 instead

x2 = linspace(oI(end),PP2.breaks(end));
ppg = BSFK(x2,g(x2),arg{:});

% Construct ppmerge structure
ppmerge = ppf;
ppmerge.breaks = [ppf.breaks(1:end-1),ppo.breaks(1:end-1),ppg.breaks(1:end)];
ppmerge.coefs = [ppf.coefs;
    ppo.coefs;
    ppg.coefs];
ppmerge.pieces = size(ppmerge.coefs,1);

% Graphical check
figure
hold on
xf = linspace(PP1.breaks(1),PP1.breaks(end));
xg = linspace(PP2.breaks(1),PP2.breaks(end));
plot(xf,f(xf),'b+');
plot(xg,g(xg),'rx');
xm = linspace(ppmerge.breaks(1),ppmerge.breaks(end));
plot(xm,ppval(ppmerge,xm),'g','LineWidth',2);
legend('f','g','merge')

A closer look in the overlap interval, it seems the transition is smooth as expected

[Michal]

I am just looking into it ... :)

Just a few preliminary remarks:

1. ppf and ppg re-fit (by BSFK) are not required, I hope. I already have these regions well covered by PP1 and PP2 (produced for example by BSFK during free-knot spline approximation of two consecutive processing window approximation - fit). Instead of ppf is possible to use directly PP1 (witt some minor modifications ... order field ... = 5 at definition, and order = 4 after BSFK re-fit), as you already mentioned in the code. Problem is with ppg, where should be first PP2 transformed by John's "shiftPolySeg" function instead of BSFK re-fit.
2. I think that your method is simply applicable on a case of multiple braks at ovelap region, too.

Am I right?

[Bruno Luong]

If you don't want to fit to rebuild ppg, you can recycle PP2 keep only tail intervals that are matter, i.e., intersecting with (10,pp2.breaks(end)) . The first kept interval you must rebuild with shifted polynomial using John's function, the rest of intervals can be concatenated as it is.

None of the steps are difficult, just do it in sorted order of the break points that are either recycled or rebuilt.

[Bruno Luong]

Same technique applied to J'ohn first example.

PP1 = mkpp([0 2],[1 0]);
PP2 = mkpp([1 3],[-1 2]);

% overlap interval
oI = [PP2.breaks(1) PP1.breaks(end)]; % [1 2;
...

[Bruno Luong]

Here is a modified version of the code that recycle PP1 and PP2 to build ppf and ppg rather than fit.

% Example as in John's answer
breaks = [4 10]; coefs = [ 0 0 1 -2 53];
PP1 = mkpp(breaks,coefs);
breaks = [8 15]; coefs = [ -1 6 1 4 77];
PP2 = mkpp(breaks,coefs);

% overlap interval
oI = [PP2.breaks(1) PP1.breaks(end)]; % [8 10];
li = diff(oI); % its length
f = @(x) ppval(PP1,x);
g = @(x) ppval(PP2,x);
% define smooth transition weight function
sigfun = @(x) x.^2.*(3-2*x); % sigmoid (0,1) -> (0,1)
w2 = @(x) sigfun((x-oI(1))/li);
w1 = @(x) 1-w2(x);
% weighted sum of f and g, applied only in oI
h = @(x) w1(x).*f(x) + w2(x).*g(x);

%% Fit on overlap interval
n = size(PP1.coefs,2);
arg = {n,[],[],struct('KnotRemoval', 'yes','Animation',0)};
xo = linspace(oI(1),oI(2),100);
ppo = BSFK(xo,h(xo),arg{:});

% Recycle PP1 and PP2 on smaller intervals distinct from overlap interval
ppf = pptruncate(PP1, -Inf, oI(1));
ppg = pptruncate(PP2, oI(2), +Inf);

% Construct ppmerge structure
breaks = [ppf.breaks(1:end-1),ppo.breaks(1:end-1),ppg.breaks(1:end)];
coefs = [ppf.coefs;
         ppo.coefs;
         ppg.coefs];
ppmerge = mkpp(breaks,coefs);

% Graphical check
figure()
hold on
xf = linspace(PP1.breaks(1),PP1.breaks(end));
xg = linspace(PP2.breaks(1),PP2.breaks(end));
plot(xf,f(xf),'b+');
plot(xg,g(xg),'rx');
xm = linspace(ppmerge.breaks(1),ppmerge.breaks(end));
plot(xm,ppval(ppmerge,xm),'g','LineWidth',2);
legend('f','g','merge');

%% Truncate a pp on a new interval [left,right]
function pp = pptruncate(pp, left, right)
breaks = pp.breaks;
coefs = pp.coefs;
b = breaks < left;
if b(1)
    n = size(coefs,2);
    i = find(b,1,'last');
    breaks(1:i-1) = [];
    coefs(1:i-1,:) = [];
    if length(breaks) >= 2
        if breaks(2) > left % check if interval not reduced to a point
            % shift the first polynomial to the new most left break point == left
            % using polynomial fit (rather than conv, see John d'Errico's shiftPolySeg())
            x = linspace(left-breaks(1), breaks(2)-breaks(1), n);
            y = polyval(coefs(1,:), x);
            coefs(1,:) = polyfit(x-x(1),y,n-1);
            breaks(1) = left;
        else
            % Interval reduced t a point, remove it
            breaks(1) = [];
            coefs(1,:) = [];
        end
    end
end
b = breaks > right;
if ~isempty(b) && b(end)
    i = find(b,1,'first');
    breaks(i+1:end) = [];
    coefs(i+1:end,:) = [];
    breaks(end) = right;
end
pp = mkpp(breaks,coefs);
end

EDIT: modification to make pptruncate more robust

[Michal]

@Bruno Luong Thanks for help! Now I would like try to create more general version of your method with multi-breaks overlap. I still wait if John put some more general solution here ... but anyway problem of PP's merging is now much more doable.

[Bruno Luong]

My method works even if the overlap intervals contains multi-break points.

You only need to compute the overlap interval

oI = [8 10]; % which is [PP2.breaks(1) PP1.breaks(end)]

It also handles smooth transition (C^2), no increases polynomial order. The tradeoff is that the polynomial in the overlap region is from fitting.

[Michal]

@Bruno Luong Nice method, really!!! The only objection is the fact that weighted sum sampling need number of samples, which is not simple to optimally estimate for the final implementation. Re-fit over overlap region is then not-optimal and produce additional computing overhead.

[Bruno Luong]

I suggest you select the sample points as the a subdivisions of f/g subintervals inside overlap region rather than 100 points as I did in my code.

[Michal]

What does it mean: "subdivison f/g subintervals inside overlap regions"?

[Michal]

OK, you probably mean the primary sampling of the signals. Am I right?

[Bruno Luong]

What I means in "subdivisions of f/g subintervals inside overlap region" is that

in the region of [8,10] in the example is the "overlap region".

In this region f has multiple break points, they define a set of "subintervals of f inside overlap region". Each of the subinterval of f is a polynomial of order k-1, since f is pp. So if you select k sample points in each sub intervals (eg. equidistance including the two end break points, thus the word "subdivision") then f is fully characterized by the values computed at this subdivision samples.

That define the sample points for f. Then do the same procedure for g. Take the union of two sets of sampling points, this gives us a way to select the points in the overlap region, and I'm sure it is large enough to recover accurately the approximation of

h = w1*f  + w2*g.

on the overlap region by fitting. It is better than selecting emprirical 100 points as I did in my code and it can capture the sharp corner where it needs, etc...

Answer 2

[Matthieu]

Are you trying to 'concatenate' your polynomial structures ?

I am not an expert, although I know those are defined by 1 vector and a matrix, breaks & coefs. break defines the range of action of the coefficients of your polynomial, stored in coefs.

You can get those vectors with:

breaks = [0 4 10 15];
coefs = [0 1 -1 1 1; 0 0 1 -2 53; -1 6 1 4 77];
pp = mkpp(breaks,coefs) ;
pp.breaks

ans = 1×4

0 4 10 15

pp.coefs

ans = 3×5

0 1 -1 1 1 0 0 1 -2 53 -1 6 1 4 77

A way of merging those would then be to concatenate the breaks and coefs of all piecewise polynomials you want to merge, as below :

breaks2 = [15 20 22 28];
coefs2  = [0 0 1 2 3; 0 0 0 0 5; -0.1 1 0 0 8] ;
breaks_merged = horzcat(breaks,breaks2(2:end));
coefs_merged  = vertcat(coefs,coefs2);
pp_merged = mkpp(breaks_merged,coefs_merged) ;

t = -5:0.1:33 ;
plot(t,ppval(pp_merged,t))

Is this what you needed ?

[Michal]

The problem is not so simple. Main problem is the fact that separate pp structures can be overlapped. So first or several break points of second pp2 structure may be located before last break points of the pp1 structure. Open question is how should be treated the merged function at the overlapped domain.

[Bruno Luong]

"Open question is how should be treated the merged function at the overlapped domain."

Uh if it is not clear for you then no one can answer at your place.

I would say you can take a weighted sum with the weight changing in the overlap region to make the transition.

Or you can fit all functions discretized by a single splibne to create a single pp But no this is still an arbitrary choice

Or ...

But no this are still an arbitrary choices

What is your intention to do with the merged function?

[John D'Errico]

Note that a weighted sum will implicitly increase the order of the polynomial segment in the merged region.

A linear weight will result in a function that is one degree higher in the overlapped interval, but it will also result in a derivative discontinuity at the break. So continuous, but not differentiable. So it should also be possible to employ a nonlinear (actually just a higher polynomial order) weighting scheme, if it were important to acheive differentiability on the overlapped interval too.

[Michal]

@Bruno Luong @John D'Errico

My main intention is to perform running (moving window) approximation of continuously measured signals by piecewise polynomials. The crucial task is suitable merging separate "pp" structures corresponding to the separate consecutive processing windows with overlap. Transition domain (overlap) between two consecutive windows can be processed by "weighted sum" method.

The final merged function (pp structure) could be one of possible way how to archive all separately measured windows at once in compressed form for next additional processing.

What I am looking for is the method how to perform this task effectively in MATLAB.

Answer 3

[John D'Errico]

Is there any method to "merge" two PP functions? Well, no. Should there be one? I've been using and writing splines tools for 40+ years, in MATLAB for 35 years or so, and answering questions about them for as long. And I've never seen that request, certainly not as you are asking. You don't write code to do something that nobody will ever want to use.

Could there be a way to do so,? Well, yes. The obvious one is to just expand the list of breaks, then append the pp segments as suggested already.

However, your question is even less, let me say, "expected", because you have an overlap. What would you do in the overlap region? Arbitrarily choose one function of the other? Take the mean? Flip a coin to decide?

Sigh. COULD you do something? Well, I suppose you could do something. I could think of several choices even here. The problem is in the overlap region, you need to resolve the conflct. That is, suppose we have intervals at the ends of domain 1 as [A,B], and at the beginning of domain 2 as [C,D]. I'll assume from what you have said is that A < C < B < D. Otherwise there is no problem. Actually, I can see another issue. What if B < C? That is, the problem is simple, if B == C. Now the direct merger already suggested works. But if B < C, then there is an undefined region. I suppose you could extrapolate the curves, over the hole in the middle between domains. But extrapolating splines is about the most obscene thing you can do. Avoid it at all costs.

Anyway, suppose we have the sceneario where A < C < B < D? We could now adjust the breaks to be exactly that sequence. Where on the interval [A,B), we use PP1 from the lower domain. On the interval [B,D), we use PP2, so the upper segment.

But on the interval [C,B) What choice do we have? Perhaps we might just use the average of the two functions on that interval. But that would now almost certainly introduce a point of discontinuity at the breaks, a bad thing.

One idea might be to use an average that varies along the interval. Essentially, this would force the curve to be a higher order segment. So instead of a cubic segment there, it would now be, a degree 4 segment at least. That would allow the resulting function to be continuous, though not necessarily differentiable. If we are willing to make it a 5th degree segment, that should allow one to construct a "merged" segment that would be both continous AND differentiable. (Just thinking off the cuff there.)

So, doable. But, for example, what would you do here?

PP1 = mkpp([0 2],[1 0]);
PP2 = mkpp([1 3],[-1 2]);

fnplt(PP1)
hold on
fnplt(PP2)
hold off

Lets see. The new pp form would have 4 breaks.

allbreaks = sort([PP1.breaks,PP2.breaks]);
S1 = [0,PP1.coefs];
S3 = [0,PP2.coefs(1),fnval(PP2,allbreaks(3))];
S2a = [PP1.coefs(1),fnval(PP1,allbreaks(2))];
S2b = PP2.coefs;
S2 = conv([-1,1],S2a) + conv([1 0],S2b);
PPmerged = mkpp(allbreaks,[S1;S2;S3]);
fnplt(PPmerged,'g')
hold on
fnplt(PP1,'r')
fnplt(PP2,'b')
hold off

So the "merged function is the same as the old ones, in the intervals [0,1], and [2,3]. nd it is made to be continuous at the overlap points, but differentiability was not enforced. On the sub-interval [1 2], the segment is now quadratic, so one degree higher order. That was necessary to enforce continuity.

As I said, doable.

[Michal]

@John D'Errico

Thanks for your comprehensive comment to my question! But I still do not fully understand your claim regarding fact that: "weighted sum will implicitly increase the order of the polynomial segment in the merged region".

Yes, standard overlap introducing the scenario: A < C < B < D

[Bruno Luong]

When you do weighted sum you must multiply the weight function by the original polynomial, so if the weight function is linear the order of the produt will increased by 1 (compared to the original polynomial), etc...

[Michal]

@Bruno Luong ... well OK ... I understand this behaviour now :) Thanks!

[Michal]

@John D'Errico Does it mean, that the followig commands realize wighted sum method over transition region (overlap)?

S2a = [PP1.coefs(1),fnval(PP1,allbreaks(2))];
S2b = PP2.coefs;
S2 = conv([-1,1],S2a) + conv([1 0],S2b);

I am a bit confused by using of "conv". From where comes vectors [-1,1] and [1 0]? These vectors represents behavior of linear weight functions over transition region?

[Bruno Luong]

Read about polyval, pofitfit, ppval, etc... to learn more about how polynomials are represented by coefficients in an array.

[-1,1] represent wthe polynomial w1(x) = -x+1

[1,0 ] represente polynomial w2(x) = x = 1 - w1(x)

conv(P,Q) represent the polynomial P*Q

[Michal]

Yes, of course I know how polynomials are represented by coefficients, but in general case variable "x" should be scaled to the interval [0,1] by x_scaled = (x - xmin)/(xmax - xmin), where [xmin,xmax] represents transition region range, Am I right?

Moreover, @John D'Errico example is able to handle only special case of 2 break points transtion region. More general case (N break points overlap) requires some more sophisticated approach. I think @John D'Errico just shows basic idea of a "doable" method.

But, may I kindly ask for a more elaborate answer?

[Bruno Luong]

" but in general case variable "x" should be scaled to the interval [0,1] by x_scaled = (x - xmin)/(xmax - xmin), where [xmin,xmax] represents transition region range, Am I right? "

No. The polynomial in pp form are polynomials of variable xi = x-x(breaks(i)

There is no scaling.

[Michal]

@John D'Errico Regarding your comment: "And I've never seen that request, certainly not as you are asking. You don't write code to do something that nobody will ever want to use." see my motivation description here.

Of course, if you know better way how to solve my problem in a case of near-real-time regime of data creation (via moving window piecewise polynomial approximator) I will be really happy to learn something new.

[Michal]

@Bruno Luong So weigths sum is always w1(x) + w2(x) = 1, but w1(x) and w2(x) are not strictly non-negative? This is a strange weights behavior...?! Or I miss something again?

[Bruno Luong]

The weights w1(x) and w2(x) have values in [0,1] in the overlap interval.

[Michal]

OK, I just slightly changed @John D'Errico example:

breaks = [4 10]; coefs = [ 0 0 1 -2 53];
PP1 = mkpp(breaks,coefs)
breaks = [8 15]; coefs = [ -1 6 1 4 77];
PP2 = mkpp(breaks,coefs)
fnplt(PP1)
hold on
fnplt(PP2)
hold off

allbreaks = sort([PP1.breaks,PP2.breaks]);
S1 = [0,PP1.coefs];
S3 = [0,PP2.coefs(1:end-1),fnval(PP2,allbreaks(3))];
S2a = [PP1.coefs(1:end-1),fnval(PP1,allbreaks(2))];
S2b = PP2.coefs;
S2 = conv([-1,1],S2a) + conv([1 0],S2b);
PPmerged = mkpp(allbreaks,[S1;S2;S3]);
fnplt(PPmerged,'g')
hold on
fnplt(PP1,'r')
fnplt(PP2,'b')
hold off

Something does not work ... But previous example works well?!

Obviously, the PPmerged.coefs are completely wrong in general, due to the change of ranges for separate regions. The example with linear functions is too simplified for further development.

[Bruno Luong]

As I told earlier, in the overlap region and all intervals of pp there is no scaling so you should not use

w1 = 1 - x
w2 = x

beside John's example.

But

w1 = 1-x/2
w2 = 1 - w1 = x/2

Meaning

S2 = conv([-1/2,1],S2a) + conv([1/2 0],S2b);

2 is the length of the overlap interval (8,10) and you have to divide the first order term with the length.

[John D'Errico]

Sorry. What I wrote works for linear segments. But it needs to be fixed for higher degree segments.

breaks = [4 10]; coefs = [ 0 0 1 -2 53];
PP1 = mkpp(breaks,coefs);
breaks = [8 15]; coefs = [ -1 6 1 4 77];
PP2 = mkpp(breaks,coefs);
fnplt(PP1)
hold on
fnplt(PP2)

allbreaks = sort([PP1.breaks,PP2.breaks]);
S1 = [0,PP1.coefs];

But S3 needs to be corrected for a higher degree segment. It was fine for a linear segment as I did it.

S3 = [0,shiftPolySeg(PP2.coefs,diff(allbreaks(2:3)))];

S2a = shiftPolySeg(PP1.coefs,diff(allbreaks(1:2)));
S2b = PP2.coefs;

The above is fine. But there was also another problem in my code, necessary when the overlap interval is not exactly ONE unit in length. The linear combination between the two pieces needs to be:

intervallen = diff(allbreaks([2,3]));
S2 = conv([-1/intervallen,1],S2a) + conv([1/intervallen 0],S2b);
PPmerged = mkpp(allbreaks,[S1;S2;S3])

PPmerged = struct with fields:
      form: 'pp'
    breaks: [4 8 10 15]
     coefs: [3×6 double]
    pieces: 3
     order: 6
       dim: 1

fnplt(PPmerged)

function Ps = shiftPolySeg(P,b)
% Given a PP segment on the interval [0,a], find the coefficients of the
% same poly on the interval [b,a], then adjusted so the new polynomial
% lives on the interval [0,a-b].
Pdeg = numel(P)-1;
Ps = [zeros(1,Pdeg),P(end)];
T = 1;
for n = 1:Pdeg
    T = conv(T,[1,b]);
    Ps = Ps + [zeros(1,Pdeg - n),P(Pdeg - n + 1)*T];
end
end

[Michal]

@Bruno Luong Sorry, but this was originally the main reason why I asked to clarify scaling properties. Your answer was ... no scaling here! And now the proper scaling is a crucial part of the method, and depends on length of overlap interval. Maybe my questions are too confusing, so sorry!

But the result is still wrong ... see attached figure:

[Bruno Luong]

The result is right in the overlap interval.

I'll let you deal with John help.

[Michal]

@John D'Errico This is a what I looking for ... Thank you very much!!!

[John D'Errico]

Yes. I fixed my code. There were two problemsi nthat I wrote originally. One was that the scaling needs to be fixed for when the overlapped interval was not 1. Secondly, the polynom,ials needed to be shifted appropriately. I wrote a little helper function that did that. Now the code works for the more general case, with any degree pp segments.

[Michal]

@John D'Errico Will this method also work well for an overlapping region composed of multiple breakpoints? I am afraid not ... could you please help me to solve this final problem?

[Bruno Luong]

Manipulating directly the breaks points and coefs as John code in the overlap interval from weighted function is combersome (though quite doable).

An esier way is

• (i) to compute a sampling of weighted sum in the overlap region, then
• (ii) do a fitting on the overlap region to get a pp form. You can even enforce the function value or even function derivative at the two ends to make the transition continuous or smooth. Also you can select differennt weigting scheme such as sigmoid, etc...
• Then (iii) concatenate pp on a partition (formed by alternate non-overlap and overlapped interval)s of the overall interval to form a global pp.

I already give a hint of such method on my earlier comment.

[John D'Errico]

Still doable.

x1 = sort(rand(1,10));
y1 = exp(x1);
spl1 = spline(x1,y1);
x2 = sort(rand(1,10)+0.25);
y2= (x2+1).^2;
spl2 = spline(x2,y2);
fnplt(spl1,'r')
hold on
fnplt(spl2,'b')
hold off

So there is some overlap, of multiple breaks.

xline(spl1.breaks,'r')
xline(spl2.breaks,'b')

Now we would need to build higher order polynomials in each segment, in much the same way as I did it before. Still hardly seems worth it, for a problem I've never seen asked about before. I'll take a look at it. But, before I spend the time to write something, can you explain what purpose you have for it?

[Michal]

John OK, the main reason why I asked for this specific algorithm is based on overall processing scenario:

1. Data are created in the fixed time processing window (typically N ~ 1000 - 10000 equidistant samples length)
2. The data are smoothed (de-noised) by application of free-knot splines approximation procedure. Most of free-knot spline algorithms has typical numerical complexity O(N^2) or O(N^3) (depends on used optimization algorithm) ... so, for increasing processing window length (N) is free-knot spline algorithm progresivelly more time consuming. So, practically is optimal processing window length N ~ 1000-3000 samples
3. Typical overlap of the two consecutive processing windows is about 50-300 samples, to be able properly process any irregular behavior of signal near to the processing windows bounds. The transiotion (overlap) region containts, in general, several break-points!!!
4. Smooted data from the step 2 are represented as PPwindow (piecewise polynomial) structure.
5. Current PPwindow is then incrementally merged with previous archived PP structure PPmerge = mergepp(PPmerge,PPwindow), where function "mergepp" represents the algorithm I am looking for.
6. Overlap - Transition region is processed via weighted sum algorithm (preserving continuity or (optionally) smoothness at the end points)
7. Finally, the structure PPmerge containts smoothed signal for many processing steps in very suitable (compressed) form to very effectivelly evaluate any part of smoothed signal.

This scenario could be possibly processed by Bruno's steps, too. But I am still not sure, how to exactly implement this method. What is your opinion about that?

[Michal]

@Bruno Luong Could you elaborate the approach you propose in more details? I think you still need to apply functuion "shiftPolySeg" proposed by John.

[Bruno Luong]

I can reform the polynomial coefs from the sampling of the function, where as in the overlap or outside it, if the left break point change. The same function that do that will be used for both cases.

[Michal]

@Bruno Luong Sorry Bruno, I do not understand what do you propose exactly. Could you add some code to be your recomendation more clear?

[Bruno Luong]

@Michal reply in my answer

[Michal]

@John D'Errico John, any progress with your method based on "shiftPolySeg" function for multi-break overlap region? Is my purpose description clear enough for you? From my point of view, is the effective and robust overlapped PP structures merging method very important for anybody who needs to apply incremental version of piecewise polynomial approximation on running data stream.

I still think, that your approach could be more general and robust than Bruno's method based on re-fit of overlap region. The Bruno's re-fit based method requires additional user defined control parameter (number of weighted sum samples ~ 100), which is not simple to "optimally" estimate. Moreover, the re-fit introduce additional computing overhead and needs to use proper fitting algorithm. Your method works only with PP structures, which is an interesting way.