(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



(* ::Input:: *)

(*********************************************************************************)
(*                              BAYESIAN ESTIMATION                              *)
(* Mathematica 6.0 Code by Tugrul Temel, 4 November 2008, IVO-Tilburg University *)
(* See the Magnus-Danilov papers for the mathematical formulation of the problem *)
(*********************************************************************************)

SetOptions["stdout",PageWidth->79] ;                                               (* Setting page width  *)

BeginPackage["BayesianEstimation`"]  ;
Context[]  ;
$ContextPath  ; 
BayesianEstimation::usage=
"BayesianEstimation Package evaluates the following Modules: rank[data,ratio,identity,vecs] tests the rank condition;IndependentQ[eqns,vecs,method] constructs a system of linearly independent equations;LinearF[eqns,vecs] constructs a function of linearly dependent equations; VarCov[ratios1,identity,vecs,method] calculates variance-covariance matrices; bayes[data,ratio,identity,vecs,method] solves the Bayesian estimation problem." ;



(* ::Input:: *)

(*********************************************************************************)
(* Data, Parameters and Definitions                                              *)
(*********************************************************************************)

X                        = ToExpression[Flatten[Flatten[data1]]] ;
n                        = Length[X] ; 

benchmark       = Partition[ToExpression[Flatten[Flatten[data2]]], 2] ;  
Nbm                   = Length[benchmark] ;
Xbm                   = Table[benchmark[[i]][[2]], {i, 1, Nbm}] ;

basicdata       = Partition[ToExpression[Flatten[Flatten[data3]]], 2] ;
p                       = Length[basicdata] ;
data                 = Rationalize[Table[basicdata[[i]][[1]] + basicdata[[i]][[2]] == 0, {i, 1, p}]] ;

ratios             = Partition[ToExpression[Flatten[Flatten[data4]]], 2] ;
m1                     = Length[ratios] ;
ratio               = Rationalize[Table[ratios[[i,2]], {i,1,m1}]] ;

ratios1           = Drop[Partition[Flatten[Flatten[data5]],4], None, 1] ;

otherpriors   = Partition[ToExpression[Flatten[Flatten[data6]]], 2] ;                         
m2                     = Length[otherpriors] ; 
other               = Rationalize[Table[otherpriors[[i]][[1]] + otherpriors[[i]][[2]] == 0, {i, 1, m2}]] ;

identities    = Partition[ToExpression[Flatten[Flatten[data7]]], 2] ;
m3                    = Length[identities] ;
identity        = Rationalize[Table[identities[[i,2]], {i,1,m3}]] ;

stddata          = Table[Partition[ToExpression[Flatten[Flatten[data8]]], 2]  [[i,2]], {i,1,p}]  ; 
stdratio        = Table[Partition[ToExpression[Flatten[Flatten[data9]]], 2]  [[i,2]], {i,1,m1}] ;
stdother        = Table[Partition[ToExpression[Flatten[Flatten[data10]]], 2]  [[i,2]], {i,1,m2}] ;

originalVarNames=Flatten[data1/.Table[ data11[[i]][[1]]-> data11[[i]][[2]], {i,1,n}  ]]  ;

methods        ={"GramSchmidt", "ModifiedGramSchmidt", "Householder", "Reorthogonalization"} ;



(* ::Input:: *)

(*********************************************************************************)
(* STEP 1: Testing the rank condition                                            *)
(*********************************************************************************)

rank[data_,ratio_,identity_,vecs_]:=
   Module[ {rkDRI,rkDR,rkRI,rkDa,rkR,rkId},
{dri,DRI}=CoefficientArrays[Flatten[{data,ratio,identity}],vecs];
{dr,DR}    =CoefficientArrays[Flatten[{data,ratio}],vecs];
{ri,RI}    =CoefficientArrays[Flatten[{ratio,identity}],vecs];
{da,Da}    =CoefficientArrays[Flatten[{data}],vecs];
{r,R}        =CoefficientArrays[Flatten[{ratio}],vecs];
{id,Id}   =CoefficientArrays[Flatten[{identity}],vecs];
 RankCondition={
                                    "rank(data,ratio,identity)=n",
                                    "rank(data,ratio)         =(p+m1)",
                                    "rank(ratio,identity)     =(m1+m3)",
                                    "rank(data)               =p",
                                    "rank(ratio)              =m1",
                                    "rank(identity)           =m3"
                               };
rkDRI=If[MatrixRank[DRI]==Length[vecs],"satisfied","NOT satisfied"];
rkDR  =If[MatrixRank[DR]  ==Length[DR],    "satisfied","NOT satisfied"];
rkRI  =If[MatrixRank[RI]  ==Length[RI],     "satisfied","NOT satisfied"];
rkDa  =If[MatrixRank[Da]  ==Length[Da],     "satisfied","NOT satisfied"];
rkR    =If[MatrixRank[R]    ==Length[R],       "satisfied","NOT satisfied"];
rkId =If[MatrixRank[Id]  ==Length[Id],     "satisfied","NOT satisfied"];
Print[
             TableForm[
                                    {{"=p+m1+m3"Length[DRI],RankCondition[[1]],MatrixRank[DRI],rkDRI},
                                    {"=p+m1"Length[DR],RankCondition[[2]],MatrixRank[DR],rkDR},
                                    {"=m1+m3"Length[RI],RankCondition[[3]],MatrixRank[RI],rkRI},
                                    {"=p"Length[Da],RankCondition[[4]],MatrixRank[Da],rkDa},
                                    {"=m1"Length[R],RankCondition[[5]],MatrixRank[R],rkR},
                                    {"=m3"Length[Id],RankCondition[[6]],MatrixRank[Id],rkId}},
                                       TableHeadings->{
                                                                          { "Full-System" ,"Data-Ratio","Ratio-Identity","Data","Ratio","Identity"},
                                                                          {"# of eqns","Rank condition tested","Rank","Test result"}
                                                                      }
                             ]
        ]
];

CellPrint[Cell["Step 1. Testing Rank Condition","Section",Background->GrayLevel[.7]]] ;
rank[data,ratio,identity,X] ;



(* ::Input:: *)

(*********************************************************************************)
(* STEP 2: Constructing a system of linearly independent equations               *)
(*********************************************************************************)

independentQ[eqns_,vecs_,method_]:=
Module[{h,V,row,col,zero,orthogonalV,selectionV},
                   {h,V}=CoefficientArrays[eqns,vecs];
                    {row,col}=Dimensions[V];
                    zero=ConstantArray[0,col];
                    orthogonalV=Orthogonalize[V,Method->method]//Chop; 
                    selectionV=Map[If[#!=zero,"keep","discard"]&,orthogonalV];    
                    reducedV={};               
                    reducedh={};
                    Do[
                            If[orthogonalV[[i]]!=zero,AppendTo[reducedV,V[[i]]]];
                            If[orthogonalV[[i]]!=zero,AppendTo[reducedh,h[[i]]]],{i,1,row}
                        ];
                    Print[
                                  TableForm[Table[{equ[i],selectionV[[i]]},{i,1,row}],TableHeadings->{{},{"Equations","ModifiedGramSchmidt"}}],
                                  V//Normal//MatrixForm, reducedV//Normal//MatrixForm,reducedh//Normal//MatrixForm
                              ]
];

systemRI=Flatten[{ratio,identity}] ;
CellPrint[Cell["Step 2. A System of Linearly Independent Equations","Section",Background->GrayLevel[.7]]] ;
independentQ[systemRI,X,methods[[2]]] ;



(* ::Input:: *)

(*********************************************************************************)
(* STEP 3: Constructing a variance-covariance matrix                             *)
(*********************************************************************************)

VarCov[ratios1_,identity_,vecs_,method_]:= 
Module[{rMat,na,Xi,Xj,Bmat1,B1,Bmat2,B2,Bmat3,B3,vardata, varratio,H3}, 
rMat        = Normal[SparseArray[Table[{IntegerPart[ratios1[[i,1]]], IntegerPart[ratios1[[i,2]]]}-> ratios1[[i,3]], {i,1,m1}], Nbm ]] ;
na            = 0 ;
Xi            = Xbm ;
Xj            = Xbm ;
Bmat1     = ConstantArray[na,{Nbm,Nbm}] ;
B1            = Table[
                                  If[Xi[[i]]!=na && Xj[[j]]!=na && rMat[[i,j]]!=na,Bmat1[[i,j]]=(rMat[[i,j]]/(1+rMat[[i,j]]^2))*Xi[[i]]+(1/(1+rMat[[i,j]]^2))*Xj[[j]],na],
                                  {i,1,Nbm},{j,1,Nbm}
                           ] ;
Bmat2   = ConstantArray[na,{Nbm,Nbm}] ;
B2         = Table[
                              If[Xi[[i]]!=na && Xj[[j]]==na && rMat[[i,j]]!=na,Bmat2[[i,j]]=(Xi[[i]]/rMat[[i,j]]),na],{i,1,Nbm},{j,1,Nbm}
                         ] ;
Bmat3  = ConstantArray[na,{Nbm,Nbm}] ;
B3        = Table[
                             If[Xi[[i]]==na && Xj[[j]]!=na && rMat[[i,j]]!=na,Bmat3[[i,j]]=Xj[[j]],na],{i,1,Nbm},{j,1,Nbm}
                         ] ;
vardata    =  stddata*stddata ; 
varratio  = stdratio*stdratio ; 
S1                = DiagonalMatrix[vardata] ;
H1                = DiagonalMatrix[Table[varratio[[i]]*B1[[IntegerPart[ratios1[[i,1]]], IntegerPart[ratios1[[i,2]]]]]^2,{i,1, m1}]]//N ;
S                  = ArrayFlatten[{{S1,0},{0,H1}}]//N ;
independentQ[identity, vecs, method]//N ;
H3              = DiagonalMatrix[ConstantArray[0,Length[reducedh]]]//N ;
H                = ArrayFlatten[{{H1,0},{0,H3}}]//N ;
Print[S1//MatrixForm,H1//MatrixForm,S//MatrixForm,H3//MatrixForm, H//MatrixForm] 
] ; 

CellPrint[Cell["Step 3. Variance-Covariance Matrices","Section",Background->GrayLevel[.7]]] ;
VarCov[ratios1,identity,X,methods[[2]]] ;



(* ::Input:: *)

(*********************************************************************************)
(*                           BAYESIAN ESTIMATION                                 *)
(*          Var-covar matrices (S1,H1,S,H3,H) were calculated above              *)
(*********************************************************************************)

bayes[data_,ratio_,identity_,vecs_,method_]:= 
 Module[{d1,D1,hh,AA,row,col,zero,orthogonalAA,selectionAA,A,h,mm,nn,Ap,S0,L,CC,K,Mean,Variance,postDist},
    {d1,D1}=CoefficientArrays[Flatten[{data}],vecs];
    {hh,AA}=CoefficientArrays[Flatten[{ratio,identity}],vecs];
    {row,col}=Dimensions[AA];
    zero=ConstantArray[0,col];
    orthogonalAA=Orthogonalize[AA,Method->method]//Chop; 
    selectionAA=Map[If[#!=zero,"keep","discard"]&,orthogonalAA];    
    A={};
    h={};
    Do[
           If[orthogonalAA[[i]]!=zero,AppendTo[A,AA[[i]]]];
           If[orthogonalAA[[i]]!=zero,AppendTo[h,hh[[i]]]],{i,1,row}
       ];
{mm,nn}  =  Dimensions[A];
Ap               =  PseudoInverse[A]//Chop//N;
S0               =  S1 + D1.Ap.H.Transpose[Ap].Transpose[D1]//Chop//N;
Which[
             mm==nn, {
                              Mean            = Ap.h- (Ap.H.Transpose[Ap]).Transpose[D1].PseudoInverse[S0].(D1.Ap.h-d1)//Chop//N,
                              Variance   = Ap.H.Transpose[Ap]-Ap.H.Transpose[Ap].Transpose[D1].PseudoInverse[S0].D1.Ap.H.Transpose[Ap]//Chop//N
                            },
             mm<nn, {
                             L                    = Transpose[Orthogonalize[NullSpace[A]]]//Chop//N ,
                             Print[{Chop[A.L]==ConstantArray[0,{mm,nn-mm}]}, {Chop[Transpose[L].L]==IdentityMatrix[nn-mm]}],
                             K                   = L.PseudoInverse[Transpose[L].Transpose[D1].PseudoInverse[S0].D1.L].Transpose[L]//Chop//N, 
                             CC                 = IdentityMatrix[n]-Ap.H.Transpose[Ap].Transpose[D1].PseudoInverse[S0].D1//Chop//N,
                             Mean            = Ap.h-(Ap.H.Transpose[Ap]+ CC.K).Transpose[D1].PseudoInverse[S0].(D1.Ap.h-d1)//Chop//N,
                            Variance = Ap.H.Transpose[Ap]-Ap.H.Transpose[Ap].Transpose[D1].PseudoInverse[S0].D1.Ap.H.Transpose[Ap]+ CC.K.Transpose[CC]//Chop//N
                         }
           ];
postDist   =  Table[{originalVarNames[[i]], Mean[[i]],Sqrt[Variance[[i,i]]]}, {i, 1,n}];
Print[ TableForm[postDist,TableHeadings->{{},{"Variable", "Post-mean","Post-se"}}]]
]

CellPrint[Cell["Bayesian Estimation","Section",Background->GrayLevel[.7]]] ;
bayes[data,ratio,identity,X,methods[[2]]] ;



(* ::Input:: *)

EndPackage[] ;
