Discrepancy analysis

[makeabhishek]

Is there a way to perform Discrepancy analysis to chose best regularization parameter for refraction manager?

[halbmy]

If you know the measurement error very well, the discrepancy principle, i.e. setting the regularization parameter so that chi-square of 1 is achieved, makes most sense, i.e. the smoothest model is selected that is able to fit the data.

[makeabhishek]

I have created the forward model, where I know the error values (as I specified it). Now I want to perform the discrepancy analysis. But I am not sure how to proceed for that while running inversion. I tried to do L-curve analysis but cannot obtain the curvature when plotting mgr.inv.phiModel() versus mgr.inv.phiData().
There fore want to look the discrepancy, which was suggested in #257. But I don't have much knowledge on the implementation of it.

Is there any example available for that?

[halbmy]

In the old (core) inversion, there was a function invChi1 but now there is None. In the upcoming inversions, there will be such a function again. Up to now you can only decrease lambda step by step until you reach chi^2=1.

[makeabhishek]

Just read an issue related to that. So thought to add that for reference for the community. #440

Want to doublecheck to obtain $chi^2$ value; I can obtain it using mgr.inv.chi2History. Is that the correct way or there is a method to obtain $chi^2$ at each iteration to set the lambda value?
I tired like that

phiM =[]; 
phiD =[];
models = [];
chiHistory = [];
velModel=[];
lam = 0

while mgr.inv.chi2History ==1:
    lam= lam+2
    pg.info("Traveltime Inversion")
    velInv = mgr.invert(data, mesh=pMesh, 
                        lam=lam,      
                        secNodes=3, 
                        showProgress= True, 
                       )
    mgr.inv.echoStatus()
    phiM.append(mgr.inv.phiModel())
    phiD.append(mgr.inv.phiData())
    chiHistory.append(mgr.inv.chi2History)

[halbmy]

1. inv.chi2History is a vector of all chi2 of the individual iterations. Use inv.chi2() instead
2. while chi2 == 1 does not make sense. I would rather suggest to use a high lambda and decrease it like while inv.chi2() > 1: self.lam = self.lam * 0.9 ...

[makeabhishek]

I tried to do that but looks, I'm missing something.
Also, how can I get the lamda value after each iteration. I tried mgr.inv.lam, but thats giving me some strange values.

chi2  = []
lam = 100
mgr = tt.TravelTimeManager(data)
while mgr.inv.chi2() >1:
    velInv = mgr.invert(data, mesh=pMesh,
                        lam = lam * 0.9,
                        secNodes=3, 
                        showProgress= True,
                        verbose=True,
                       )
    chi2.append(mgr.inv.chi2())
    print(chi2)
    print(lam)

How can I initialise the mgr.inv.chi2(). Do I need to run an inversion first before the loop?

[halbmy]

Of course you have to run an inversion in order to have a data fit. It could be done with maxIter=0 (infinite lambda).

[makeabhishek]

But if I just initialize the inversion it will give the chi2 value less than 1 in the first inversion and stop the while loop after the first iteration. Isn't it?

chi2 = []
phiM =[]; 
phiD =[];
lam = 5000
print('initial Lmbda = ', lam)
mgr = tt.TravelTimeManager(data)
velInv = mgr.invert(data, mesh=pMesh,lam =1000,secNodes=3,maxIter=0)
while mgr.inv.chi2() >1:
    lam = lam * 0.9
    print('Lambda = ', lam)
    model= mgr.invert(data, mesh=pMesh,
                        lam = lam,
                        secNodes=3, 
                        showProgress= True,
                       )
    chi2.append(mgr.inv.chi2())
    phiM.append(mgr.inv.phiModel())
    phiD.append(mgr.inv.phiData())
    print(chi2)
    print(lam)

plt.plot(phiM,phiD)

[halbmy]

Then you should start with a higher lambda value going towards a homogeneous halfspace that should not fit the data.

[makeabhishek]

Not sure why, but it doesn't matter what lambda parameter I'm choosing chi² is always converging to less than 1 in the first iteration.
For example the output from the inversion at very small lambda

--------------------------------------------------------------------------------
inv.iter 1 ... chi² = 602.77 (dPhi = 76.78%) lam: 9.75
--------------------------------------------------------------------------------
inv.iter 2 ... chi² = 18.16 (dPhi = 96.97%) lam: 9.75
--------------------------------------------------------------------------------
inv.iter 3 ... chi² = 4.96 (dPhi = 72.53%) lam: 9.75
--------------------------------------------------------------------------------
inv.iter 4 ... chi² = 3.16 (dPhi = 36.05%) lam: 9.75
--------------------------------------------------------------------------------
inv.iter 5 ... chi² = 1.55 (dPhi = 50.97%) lam: 9.75
--------------------------------------------------------------------------------
inv.iter 6 ... chi² = 1.15 (dPhi = 25.38%) lam: 9.75
--------------------------------------------------------------------------------
inv.iter 7 ... chi² = 0.89 (dPhi = 22.73%) lam: 9.75

[halbmy]

If chi goes below 1 in the first iteration for a high lambda value, your anomalous data probably Do Not exceed the noise Level.
Sorry, but nobody can help you unless you share the whole Story with the community.

[makeabhishek]

I created a synthetic forward model as shown below and generated the data using traveltime manger. Than I'm trying to perform discrepancy analysis on the data. The model and data is shown below.

Forward modelling

mgr =  tt.TravelTimeManager()
ttData = mgr.simulate(mesh_fwd,
                      scheme=ttScheme, 
                      slowness=1/vel,
                      noiseLevel=0.01, 
                      secNodes=4,
                      seed=1337)

Discrepancy analysis

chi2 = []
phiM =[]; 
phiD =[];
lam = 10000
tolerance = 1.05
print('initial Lmbda = ', lam)
mgr = tt.TravelTimeManager(data)

velInv = mgr.invert(data, mesh=pMesh,lam =50000,secNodes=3,maxIter=0)

while mgr.inv.chi2() > tolerance:
    lam = lam / 0.95
    print(f'Running inversion with Lambda = {lam}')
    velInv = mgr.invert(data, mesh=pMesh,
                        lam = lam,
                        secNodes=3, 
                        showProgress= True,
                        verbose=True,
                       )
    chi2_value = mgr.inv.chi2()
    chi2.append(chi2_value)
    phiM.append(mgr.inv.phiModel())
    phiD.append(mgr.inv.phiData())
    print('Chi2 :', chi2_value)
    if len(chi2)>1 and abs(chi2[-1]-chi2[-2]) <1e-3:
        print("Chi2 not changing significantly, stopping!")
        break
        
plt.plot(phiM,phiD,'-o')
plt.xlabel('Model Regularization (PhiM)')
plt.ylabel('Data Misfit (PhiD)')
plt.title('Discrepancy Analysis')
plt.show()

Result of discrepancy analysis is shown in the last post.

[halbmy]

Your parameter contrast is so low that everything falls below the noise level as you can clearly see in the data image. No chance to resolve it.