Skip to content

Commit

Permalink
Merge pull request #452 from mrclary/extend-uncertainty
Browse files Browse the repository at this point in the history
Evaluate uncertainty for arbitrary independent variable value
  • Loading branch information
newville authored Feb 8, 2018
2 parents fc5e38d + 198adb5 commit 1affbb9
Showing 1 changed file with 10 additions and 8 deletions.
18 changes: 10 additions & 8 deletions lmfit/model.py
Original file line number Diff line number Diff line change
Expand Up @@ -1298,28 +1298,30 @@ def eval_uncertainty(self, params=None, sigma=1, **kwargs):
give the same results, within precision errors.
"""
self.userkws.update(kwargs)
userkws = self.userkws.copy()
userkws.update(kwargs)
if params is None:
params = self.params

nvarys = self.nvarys
ndata = self.ndata
covar = self.covar / self.redchi
fjac = np.zeros(ndata*nvarys).reshape((nvarys, ndata))
# ensure fjac and df2 are correct size if independent var updated by kwargs
ndata = self.model.eval(self.params, **userkws).size
covar = self.covar
fjac = np.zeros((nvarys, ndata))
df2 = np.zeros(ndata)

# find derivative by hand!
pars = self.params.copy()
for i in range(nvarys):
pname = self.var_names[i]
pars = self.params
val0 = pars[pname].value
dval = pars[pname].stderr/3.0

pars[pname].value = val0 + dval
res1 = self.model.eval(pars, **self.userkws)
res1 = self.model.eval(pars, **userkws)

pars[pname].value = val0 - dval
res2 = self.model.eval(pars, **self.userkws)
res2 = self.model.eval(pars, **userkws)

pars[pname].value = val0
fjac[i] = (res1 - res2) / (2*dval)
Expand All @@ -1332,7 +1334,7 @@ def eval_uncertainty(self, params=None, sigma=1, **kwargs):
prob = sigma
else:
prob = erf(sigma/np.sqrt(2))
return np.sqrt(df2*self.redchi) * t.ppf((prob+1)/2.0, ndata-nvarys)
return np.sqrt(df2) * t.ppf((prob+1)/2.0, self.ndata-nvarys)

def conf_interval(self, **kwargs):
"""Calculate the confidence intervals for the variable parameters.
Expand Down

0 comments on commit 1affbb9

Please sign in to comment.