forked from latz-io/bayesian_inversion
-
Notifications
You must be signed in to change notification settings - Fork 0
/
IS_t2N.m
46 lines (36 loc) · 937 Bytes
/
IS_t2N.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
% Importance Sampling: Works.
% Simulate the 2nd Moment of a normal distribution given a t-distributed
% samples
% Jonas Latz, M.Sc.
% Lehrstuhl f?r Numerische Mathematik
% Fakult?t f?r Mathematik
% Technische Universit?t M?nchen
% 2017 -
%% Configuration
k = 3.5;
true = 1;
%Number of Simulations
N = 100000
%% Sample
X = trnd(k,[N,1]);
Y= X.^2;
w = normpdf(X,0,1)./tpdf(X,3);
w = w/sum(w);
Est = sum(Y.*w);
%% 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(4)
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,'Linewidth',1.2);
hold off
%% Return the relative error
rel_error = abs(Est-true)/true