From 86e86dc115a2d743ce99dc7cbc6731f39303567a Mon Sep 17 00:00:00 2001
From: Sourish <sourish@cmi.ac.in>
Date: Tue, 11 Oct 2022 21:26:59 +0530
Subject: [PATCH] add user specific alpha and beta

---
 src/CRRao.jl                      |   4 +-
 src/bayesian/linear_regression.jl | 103 ++++++++++++++++++++++++++++++
 2 files changed, 105 insertions(+), 2 deletions(-)

diff --git a/src/CRRao.jl b/src/CRRao.jl
index 2769b04..330fede 100644
--- a/src/CRRao.jl
+++ b/src/CRRao.jl
@@ -41,7 +41,7 @@ Type representing the Poisson Regression model class.
 struct PoissonRegression end
 
 
-
+struct Prior_Gauss end
 struct Prior_Ridge end
 struct Prior_Laplace end
 struct Prior_Cauchy end
@@ -135,7 +135,7 @@ end
 Cauchit() = Cauchit(Cauchit_Link)
 
 export LinearRegression, LogisticRegression, PoissonRegression, NegBinomRegression
-export Prior_Ridge, Prior_Laplace, Prior_Cauchy, Prior_TDist, Prior_Uniform, Prior_HorseShoe
+export Prior_Ridge, Prior_Laplace, Prior_Cauchy, Prior_TDist, Prior_Uniform, Prior_HorseShoe, Prior_Gauss
 export CRRaoLink, Logit, Probit, Cloglog, Cauchit, fit
 export coeftable, r2, adjr2, loglikelihood, aic, bic, sigma, predict, residuals, cooksdistance
 export FrequentistRegression, BayesianRegression
diff --git a/src/bayesian/linear_regression.jl b/src/bayesian/linear_regression.jl
index d8a241d..9d49cb3 100644
--- a/src/bayesian/linear_regression.jl
+++ b/src/bayesian/linear_regression.jl
@@ -497,5 +497,108 @@ function fit(
         y ~ MvNormal(α .+ X * β, σ)
     end
 
+    return linear_reg(formula, data, LinearRegression, sim_size)
+end
+
+"""
+```julia
+fit(formula::FormulaTerm, data::DataFrame, modelClass::LinearRegression, prior::Prior_Gauss,alpha_prior_mean::Float64 = 0.0, beta_prior_mean::Float64, sim_size::Int64 = 1000, h::Float64 = 0.1)
+```
+
+Fit a Bayesian Linear Regression model on the input data with a Gaussian prior with user specific prior mean for α and β. User doesnot have
+    idea about the prior sd of α and β. User ignore the specification for sd of α and β.
+
+# Example
+```julia-repl
+julia> using CRRao, RDatasets, StableRNGs, StatsModels
+julia> CRRao.set_rng(StableRNG(123));
+julia> df = dataset("datasets", "mtcars");                                                                                                
+julia> container = fit(@formula(MPG ~ HP + WT + Gear), df, LinearRegression(), Prior_Gauss(),0.0,[0.0,-3.0,1.0],1000)
+```
+"""
+function fit(
+    formula::FormulaTerm
+    , data::DataFrame
+    , modelClass::LinearRegression
+    , prior::Prior_Gauss
+    , alpha_prior_mean::Float64
+    , beta_prior_mean::Vector{Float64}
+    , sim_size::Int64 = 1000
+    , h::Float64 = 0.1
+)
+    @model LinearRegression(X, y) = begin
+        p = size(X, 2)
+        α0 = alpha_prior_mean
+        β0 = beta_prior_mean
+
+        #priors
+        a0 = 0.1
+        b0 = 0.1
+
+        Ip = 1*Matrix(I,p,p)
+
+        S = cov(X)+Ip
+
+        v ~ InverseGamma(h, h)
+        σ ~ InverseGamma(a0, b0)
+        α ~ Normal(α0, v * σ)
+        β ~ MvNormal(β0, S)
+
+        #likelihood
+        y ~ MvNormal(α .+ X * β, σ)
+    end
+
+    return linear_reg(formula, data, LinearRegression, sim_size)
+end
+
+
+"""
+```julia
+fit(formula::FormulaTerm, data::DataFrame, modelClass::LinearRegression, prior::Prior_Gauss, alpha_prior_mean::Float64, alpha_prior_sd::Float64, beta_prior_mean::Vector{Float64}, beta_prior_sd::Vector{Float64}, sim_size::Int64 = 1000)
+```
+
+Fit a Bayesian Linear Regression model on the input data with a Gaussian prior with user specific prior mean and sd for α and β. 
+
+# Example
+```julia-repl
+julia> using CRRao, RDatasets, StableRNGs, StatsModels
+julia> CRRao.set_rng(StableRNG(123));
+julia> df = dataset("datasets", "mtcars");                                                                                                
+julia> container = fit(@formula(MPG ~ HP + WT + Gear), df, LinearRegression(), Prior_Gauss(),30.0,10.0,[0.0,-3.0,1.0],[0.1,1.0,1.0],1000)
+```
+"""
+function fit(
+    formula::FormulaTerm
+    , data::DataFrame
+    , modelClass::LinearRegression
+    , prior::Prior_Gauss
+    , alpha_prior_mean::Float64
+    , alpha_prior_sd::Float64
+    , beta_prior_mean::Vector{Float64}
+    , beta_prior_sd::Vector{Float64}
+    , sim_size::Int64 = 1000
+)
+    @model LinearRegression(X, y) = begin
+        p = size(X, 2)
+        α0 = alpha_prior_mean
+        σ_α0 = alpha_prior_sd
+        β0 = beta_prior_mean
+        σ_β0 = beta_prior_sd
+
+        S = Matrix(Diagonal(σ_β0))
+        S = S*S
+
+        #priors
+        a0 = 0.1
+        b0 = 0.1
+
+        σ ~ InverseGamma(a0, b0)
+        α ~ Normal(α0, σ_α0)
+        β ~ MvNormal(β0, S)
+
+        #likelihood
+        y ~ MvNormal(α .+ X * β, σ)
+    end
+
     return linear_reg(formula, data, LinearRegression, sim_size)
 end
\ No newline at end of file