estimate normal learning model and see how parameters update over time
2 views (last 30 days)
Show older comments
I want to estimate a normal learning model by simulating a variable via draws from a normal distribution and get the posterior, expected value, and variance of the expected value. A baseline I would like to run is this: https://www2.bcs.rochester.edu/sites/jacobslab/cheat_sheet/bayes_Normal_Normal.pdf
Here is my sample code I have from text. However, it starts with a linear regression. How can I adapt this so it just is a simple case where I am sampling from one variable? Do I just set this up so x=mu and xi=mu in the log likelihood?
Code:
n=831741; k=3; % set number of obs and n number of vars
y = data.Pi;
x = data.X; b = ones(k,1); % generate data set
y = x*b + randn(n,1);
r = [1.0 1.0 1.0]'; % prior means
R = eye(k); T = eye(k); % prior variance
Q = chol(inv(T)); q = Q*r;
b0 = (x'*x)\(x'*y); % use ols as initial values
sige = (y-x*b0)'*(y-x*b0)/(n-k);
xpx = x'*x; xpy = x'*y; % calculate x’x, x’y (just once)
qpq = Q'*Q; qpv = Q'*q; % calculate Q’Q, Q’q (just once)
ndraw = 5000;
bsave = zeros(ndraw,k); % allocate storage for results
ssave = zeros(ndraw,1);
for i=1:ndraw % Start the sampling
xpxi = inv(xpx + sige*qpq);
b = xpxi*(xpy + sige*qpv); % update b
b = norm_rnd(sige*xpxi) + b; % draw MV normal with mean(b), var(b)
bsave(i,:) = b'; % save b draws
e = y - x*b; ssr = e'*e; % update sige
chi = chis_rnd(1,n); % do chisquared(n) draw
sige = ssr/chi;
ssave(i,1) = sige; % save sige draws
end % End the sampling
1 Comment
Ayush Anand
on 13 Sep 2023
Hi Joshua, I don't exactly understand what do you mean by "sampling from just one variable". If you mean that you do not want the linear regression setup, you can use the "normrnd" function to generate Gaussian distributed numbers. You can also use the "normpdf" function to generate Gaussian distributed numbers from a chosen vector in accordance with the pdf. Let me know if I am misunderstanding your question or if you are trying to do something else
Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!