Added fit function with prior, plus bugfix

src/FitRoutines.jl

@@ -6,6 +6,7 @@ using ADerrors, LsqFit, LinearAlgebra, ForwardDiff
 include("types.jl")
 include("pvalue.jl")
 include("fit_functions.jl")
+include("fit_routine_prior.jl")
 include("find_fit.jl")
 include("fit_scan.jl")
 include("utils.jl")

src/fit_functions.jl

@@ -10,7 +10,7 @@ It assumes that `x`, `data` and `W` have compatible dimensions.
 """
 function get_chi(model,x,par,data,W::AbstractMatrix)
     ytrue = model(x,par)
-    return sum((data[i]-ytrue[i]) * W[i,j] * (data[j]-model(x[j],par)) for i in axes(W,1), j in axes(W,2))
+    return sum((data[i]-ytrue[i]) * W[i,j] * (data[j]-ytrue[j]) for i in axes(W,1), j in axes(W,2))
 end
 
 function get_chi(model,x,par,data,W::AbstractVector)
@@ -84,7 +84,7 @@ function fit_routine(model::Function,
 
     nobs = length(ydata)
     if length(xdata)!=nobs
-        error("xdata and ydata must be the of same size")
+        error("xdata and ydata must be of the same size")
     end
 
     [uwerr(y, wpm) for y in ydata]
@@ -270,8 +270,6 @@ function fit_routine(models::AbstractArray{<:Function},
     end
 
     fit = fit_routine(Model,X,Y,npar,guess = guess,W=WW,C=CC,corr=corr,logfile=NoPrint(),wpm=wpm)
-
-
     compute_error_for_logging(logfile, fit.par)
     flag = typeof(xdata) <: AbstractArray{<:AbstractArray{ADerrors.uwreal}}
     println(logfile, "fit_routine output:")

src/fit_routine_prior.jl

+@doc raw"""
+
+y = model(x,p,prior)
+
+"""
+function fit_routine_prior(model::Function,xdata::AbstractArray{<:Real},ydata::AbstractArray{uwreal},prior::AbstractArray{uwreal},npar::Int64;
+                           wpm::Union{Dict{Int64,Vector{Float64}}, Dict{String,Vector{Float64}}} = Dict{Int64,Vector{Float64}}(),
+                           corr::Bool = true,
+                           W::AbstractArray{Float64}=Vector{Float64}(),
+                           guess::AbstractVector{Float64} = fill(0.5,npar),
+                           logfile = NoPrint(),
+                           C::AbstractArray{Float64} = Matrix{Float64}(undef,0,0))
+
+    nobs = length(ydata)
+    nprior = length(prior)
+    length(xdata)==nobs || error("xdata and ydata musth be of the same size")
+
+    function F(i::Vector,par)
+        res = Vector{Any}(undef,nobs+nprior)
+        res[1:nobs] .= model(xdata,par[1:npar],par[npar+1:npar+nprior])
+        res[nobs+1:nobs+nprior] .= par[npar+1:npar+nprior]
+        return res
+    end
+    guess = vcat(guess,value.(prior))
+    y = vcat(ydata,prior)
+    fit = fit_routine(F, collect(eachindex(y)), y, length(guess),
+                      guess=guess, W=W, C=C, corr=corr,
+                      logfile=NoPrint(),
+                      wpm=wpm)
+    return fit
+end

test/prior.jl

+using Revise, ADerrors, FitRoutines, PyPlot
+
+xdata = collect(1:10)
+true_model(x) = 1.0 + 3exp(-x/5)
+ydata = let
+    y = true_model.(xdata)
+    [uwreal(randn(1000).+_y, "test") for _y in y]
+end
+
+Q = uwreal([5.,0.1], "prior")
+
+model(x,p,q) = @. p[1] + p[2]*exp(-x/q[1])
+model(x,p)= p[1] + p[2]*exp(-x/p[3])
+fit = FitRoutines.fit_routine_prior(model,xdata,ydata,[Q],2,corr=true)
+
+let
+    fig,ax =subplots(1)
+    y,dy = value.(ydata), ADerrors.err.(ydata)
+    ax.errorbar(xdata,y,dy,color="k",marker="s",fillstyle="none",ls="")
+    r = collect(1:0.1:10)
+    ax.plot(r,true_model.(r),color="r",ls="-")
+    y = [model(x,fit.par) for x in r]
+    uwerr.(y)
+    y,dy = value.(y), ADerrors.err.(y)
+    ax.fill_between(r,y.-dy, y.+dy,color="g",alpha=0.5)
+end