Commit 5df8d28c authored by Michael Rudolf's avatar Michael Rudolf
Browse files

fixed VST fitting

Fixed weird behaviour of bootstrap fitting for VST tests
Switched to scipy.optimize.curve_fit for this case.
parent c7e0b4b3
...@@ -18,7 +18,8 @@ import matplotlib.gridspec as gridspec ...@@ -18,7 +18,8 @@ import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import nptdms import nptdms
import numpy as np import numpy as np
from scipy import optimize, stats from scipy.optimize import least_squares, curve_fit
from scipy.stats import norm
from scipy.signal import medfilt from scipy.signal import medfilt
...@@ -395,12 +396,12 @@ def fitfric(x, a, b): ...@@ -395,12 +396,12 @@ def fitfric(x, a, b):
# %%=====================VELOCITY STEPPING ANALYSIS============================ # %%=====================VELOCITY STEPPING ANALYSIS============================
def vst_analysis(data): def vst_analysis(data):
mu_app = data['shearstress']/data['normalstress'] # Apparent friction mu_app = data['shearstress']/data['normalstress'] # Apparent friction
pfit, perr = fit_bootstrap( pfit, pcov = curve_fit(
[-.1, .5], fitfric,
np.log10(data['velocity']), np.log10(data['velocity']),
mu_app, mu_app
function=fitfric
) )
perr = 2*np.sqrt(np.diag(pcov))
name_fit = ( name_fit = (
r'y = ({:.4f}$\pm${:.5f})log10(x) +' + r'y = ({:.4f}$\pm${:.5f})log10(x) +' +
'\n\t\t' + '\n\t\t' +
...@@ -434,7 +435,7 @@ def fit_bootstrap(p0, datax, datay, function=None, ...@@ -434,7 +435,7 @@ def fit_bootstrap(p0, datax, datay, function=None,
return function(x, *p) - y return function(x, *p) - y
# Fit first time # Fit first time
pfit, perr = optimize.leastsq(errfunc, pfit, perr = leastsq(errfunc,
p0, p0,
args=(datax, datay), args=(datax, datay),
full_output=0) full_output=0)
...@@ -453,7 +454,7 @@ def fit_bootstrap(p0, datax, datay, function=None, ...@@ -453,7 +454,7 @@ def fit_bootstrap(p0, datax, datay, function=None,
randomdataY = datay + randomDelta randomdataY = datay + randomDelta
randomfit, randomcov = \ randomfit, randomcov = \
optimize.leastsq(errfunc, p0, args=(datax, randomdataY), leastsq(errfunc, p0, args=(datax, randomdataY),
full_output=0) full_output=0)
ps.append(randomfit) ps.append(randomfit)
...@@ -545,10 +546,10 @@ def plothist(path, name, strength, data_mut): ...@@ -545,10 +546,10 @@ def plothist(path, name, strength, data_mut):
color='royalblue', color='royalblue',
edgecolor='black') edgecolor='black')
lnspc = np.linspace(np.nanmin(coef), np.nanmax(coef), len(coef)) lnspc = np.linspace(np.nanmin(coef), np.nanmax(coef), len(coef))
histfit_coef = stats.norm.pdf( histfit_coef = norm.pdf(
lnspc, lnspc,
stats.norm.fit(coef[~np.isnan(coef)])[0], norm.fit(coef[~np.isnan(coef)])[0],
stats.norm.fit(coef[~np.isnan(coef)])[1]) norm.fit(coef[~np.isnan(coef)])[1])
axrow[0].plot(lnspc, histfit_coef, axrow[0].plot(lnspc, histfit_coef,
'r--', 'r--',
linewidth=2, linewidth=2,
...@@ -557,10 +558,10 @@ def plothist(path, name, strength, data_mut): ...@@ -557,10 +558,10 @@ def plothist(path, name, strength, data_mut):
axrow[0].set_xlabel('Friction coefficient $\\mu$') axrow[0].set_xlabel('Friction coefficient $\\mu$')
axrow[0].set_ylabel('Counts') axrow[0].set_ylabel('Counts')
text = 'Mean: ' + \ text = 'Mean: ' + \
str(round(stats.norm.fit(coef[~np.isnan(coef)])[0], 3)) + \ str(round(norm.fit(coef[~np.isnan(coef)])[0], 3)) + \
'\n' + \ '\n' + \
'Std.: ' + \ 'Std.: ' + \
str(round(stats.norm.fit(coef[~np.isnan(coef)])[1], 3)) + \ str(round(norm.fit(coef[~np.isnan(coef)])[1], 3)) + \
'\n' + \ '\n' + \
' (' + \ ' (' + \
str(coef[~np.isnan(coef)].size) + \ str(coef[~np.isnan(coef)].size) + \
...@@ -575,21 +576,21 @@ def plothist(path, name, strength, data_mut): ...@@ -575,21 +576,21 @@ def plothist(path, name, strength, data_mut):
density=True, density=True,
color='royalblue', color='royalblue',
edgecolor='black') edgecolor='black')
statscoh = stats.norm.fit(coh[~np.isnan(coh)]) statscoh = norm.fit(coh[~np.isnan(coh)])
lnspc = np.linspace(np.nanmin(coh), np.nanmax(coh), len(coh)) lnspc = np.linspace(np.nanmin(coh), np.nanmax(coh), len(coh))
histfit_coh = stats.norm.pdf( histfit_coh = norm.pdf(
lnspc, lnspc,
stats.norm.fit(coh[~np.isnan(coh)])[0], norm.fit(coh[~np.isnan(coh)])[0],
stats.norm.fit(coh[~np.isnan(coh)])[1]) norm.fit(coh[~np.isnan(coh)])[1])
axrow[1].plot(lnspc, histfit_coh, 'r--', linewidth=2) axrow[1].plot(lnspc, histfit_coh, 'r--', linewidth=2)
axrow[1].set_title(tit_C) axrow[1].set_title(tit_C)
axrow[1].set_xlabel('Cohesion $C$ [Pa]') axrow[1].set_xlabel('Cohesion $C$ [Pa]')
axrow[1].set_ylabel('Counts') axrow[1].set_ylabel('Counts')
text = 'Mean: ' + \ text = 'Mean: ' + \
str(round(stats.norm.fit(coh[~np.isnan(coh)])[0], 2)) + \ str(round(norm.fit(coh[~np.isnan(coh)])[0], 2)) + \
'\n' + \ '\n' + \
'Std.: ' + \ 'Std.: ' + \
str(round(stats.norm.fit(coh[~np.isnan(coh)])[1], 2)) + \ str(round(norm.fit(coh[~np.isnan(coh)])[1], 2)) + \
'\n' + \ '\n' + \
' (' + \ ' (' + \
str(coh[~np.isnan(coh)].size) + \ str(coh[~np.isnan(coh)].size) + \
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment