forked from latz-io/bayesian_inversion
-
Notifications
You must be signed in to change notification settings - Fork 0
/
IS_N2t.m
51 lines (41 loc) · 1.04 KB
/
IS_N2t.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
% Importance Sampling: Breaks.
% Simulate the 2nd Moment of a student's t-distribution given a normal
% proposals
% Jonas Latz, M.Sc.
% Lehrstuhl f?r Numerische Mathematik
% Fakult?t f?r Mathematik
% Technische Universit?t M?nchen
% 2017 -
%% Configurations
% Degrees of Freedom
%k = 3.5;
k = 200000;
% Exact 2nd Moment
true = k/(k-2);
%Number of Simulations
N = 100000;
%% Sample
X = normrnd(0,1,[N,1]);
Y = X.^2;
w = tpdf(X,k)./normpdf(X,0,1);
Est = sum(Y.*w)/sum(w);
%% Normalise the weights (just used to plot the the Dist)
Norm = sum(w);
w = w/Norm;
%% Plot
% Plot requires the function histwv.m
% (Copyright (c) 2016, Brent All rights reserved.)
% that can be downloaded:
% https://de.mathworks.com/matlabcentral/fileexchange/58450-histwv-v--w--min--max--bins-/content/histwv.m
figure(5)
X_vals = min(X):0.5:max(X);
[histw, intervals] = histwv(X',w',min(X),max(X),length(X_vals));
bar(X_vals,histw)
hold on
x = min(X):0.05:max(X);
y = tpdf(x,k);
plot(x,y);
hold off
%% Return the relative error
rel_error = abs(Est-true)/true