-
Notifications
You must be signed in to change notification settings - Fork 0
/
MBHTM.m
83 lines (60 loc) · 2.19 KB
/
MBHTM.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
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
function [y] = MBHTM(u,m3,m4)
% GOAL
% Transform a gaussian-distributed variable into a non-gaussian one,
% given a target skewness (m3) and kurtosis (m4).
% INPUT
% u: random gaussian vector : size : [ 1 x N]
% m3: Target Skewness; is [1x1]
% m4: Target Kurtosis; is [1x1]
% OUTPUT
% y : Non gaussian variable with target m1,m2,m3 and m4. size : [ 1 x N]
% example: [y] = MBHTM(randn(1,1e5),0.5,4)
%% Reduction of the input
% variable used in input is reduced into a gaussian variable of variance equal to
% 1 and mean equal to 0:
m1 = nanmean(u);
m2 = nanvar(u);
u = (u-m1)/sqrt(m2) ;
% if the target skewness and kurtosis are the same as for a gaussian
% distribution, the output is the same as the input
if and(m3==0,m4 ==3),
y = m1+sqrt(m2).*u;
return
end
%% Test of applicability of the transformation
Ye = m4-3;% excess coefficient
test = (1/3.*(1-1.5.*Ye-sqrt(1+1.5*Ye))).^2.*(6./(sqrt(1+1.5*Ye)-1)-1);
if (m3^2> test||isnan(test)),
fprintf(' applicability condition is not verified \n');
fprintf([num2str(test,3),' < m3^2 = ',num2str(m3^2,2),' \n']);
% the input is NaN;
y = NaN(1,length(u));
return
else
fprintf('applicability condition is verified \n');
fprintf([num2str(test,3),' >= m3^2 = ',num2str(m3^2,2),' \n']);
end
%% Main calculation
% guess of h3 and h4 before solving the system of equations_
h3_guess = m3./(4+2.*sqrt(1+1.5*Ye));
h4_guess = (sqrt(1+1.5*Ye)-1)./18;
guess = [h3_guess,h4_guess];
% Call solver
options = optimoptions('fsolve','TolX',1e-10); % Option to display output
[x,~] = fsolve(@(x) myPolyFun(x,m3,m4),guess,options);
% get h3 and h4
h3=x(1);
h4=x(2);
% get a,b and alpha
a = h3/(3*h4);
b = 1/(3*h4);
alpha = 1./(sqrt(1+2*h3^2+6*h4.^2));
% get the non gaussian variable y
y = alpha./b.*(u.^3/3+a.*u.^2+(b-1).*u-a);
y = m1+sqrt(m2).*y;
function F = myPolyFun(x,Y3,Y4)
F= [-Y3+(2.*x(1).*(3+4.*x(1).^2+18.*x(2)+54.*x(2).^2))./(sqrt(1+2*x(1).^2+6*x(2).^2)).^3;
-Y4+3.*(1+20.*x(1).^4+8.*x(2)+84.*x(2).^2+432.*x(2).^3 +1116.*x(2).^4+(4.*x(1).^2.*(5+48.*x(2)+168.*x(2).^2)))./...
(1+2*x(1).^2+6*x(2).^2).^2];
end
end