-
Notifications
You must be signed in to change notification settings - Fork 7
/
Copy pathExtrVal2.m
121 lines (103 loc) · 4.05 KB
/
ExtrVal2.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
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
classdef ExtrVal2 < dContinuous
% Extreme Value Type II distribution with parameters mu, scale, shapek
% Reference: JKN, Vol 2, page 3. (What they call epsilon I call mu, what they call theta I call scale)
% shapek is the shape parameter
% Also mentioned by Luce (1986), page 508.
properties(SetAccess = protected)
mu, scale, shapek
end
methods (Static)
function Reals = ParmsToReals(Parms,~)
Reals = [Parms(1) NumTrans.GT2Real(eps,Parms(2)) NumTrans.GT2Real(eps,Parms(3))];
end
function Parms = RealsToParms(Reals,~)
Parms = [Reals(1) NumTrans.Real2GT(eps,Reals(2)) NumTrans.Real2GT(eps,Reals(3))];
end
end
methods
function obj=ExtrVal2(varargin)
obj=obj@dContinuous('ExtrVal2');
obj.ParmTypes = 'rrr';
obj.DefaultParmCodes = 'rrr';
obj.NDistParms = 3;
obj.SearchOptions.MaxFunEvals = 2000;
obj.SearchOptions.MaxIter = 2000;
switch nargin
case 0
case 3
ResetParms(obj,[varargin{:}]);
otherwise
ME = MException('ExtrVal2:Constructor', ...
'ExtrVal2 constructor needs 0 or 3 arguments.');
throw(ME);
end
end
function []=ResetParms(obj,newparmvalues)
ClearBeforeResetParmsC(obj);
obj.mu = newparmvalues(1);
obj.scale = newparmvalues(2);
obj.shapek = newparmvalues(3);
ReInit(obj);
end
function PerturbParms(obj,ParmCodes)
newmu = ifelse(ParmCodes(1)=='f', obj.mu,obj.mu+obj.scale/100);
newscale = ifelse(ParmCodes(2)=='f', obj.scale, 0.95*obj.scale);
newshapek = ifelse(ParmCodes(3)=='f', obj.shapek, 1.05*obj.shapek);
obj.ResetParms([newmu newscale newshapek]);
end
function []=ReInit(obj)
assert(obj.scale>0,'ExtrVal2 scale must be > 0.');
assert(obj.shapek>0,'ExtrVal2 shapek must be > 0.');
obj.Initialized = true;
obj.LowerBound = InverseCDF(obj,obj.CDFNearlyZero);
obj.UpperBound = InverseCDF(obj,obj.CDFNearlyOne);
if (obj.NameBuilding)
BuildMyName(obj);
end
end
function thispdf=PDF(obj,X)
[thispdf, InBounds, Done] = MaybeSplinePDF(obj,X);
if Done
return;
end
XX = (X(InBounds) - obj.mu) / obj.scale;
XXP = XX.^(-obj.shapek);
thispdf(InBounds) = obj.shapek*XXP.*exp(-XXP)./(XX*obj.scale);
% for i=1:numel(X)
% if (X(i) > obj.LowerBound) && (X(i) < obj.UpperBound)
% XX = (X(i) - obj.mu) / obj.scale;
% XXP = XX^(-obj.shapek);
% thispdf(i) = obj.shapek*XXP*exp(-XXP)/(XX*obj.scale);
% end
% end
end
function thiscdf=CDF(obj,X)
[thiscdf, InBounds, Done] = MaybeSplineCDF(obj,X);
if Done
return;
end
XX = (X(InBounds) - obj.mu) / obj.scale;
XX = XX.^(-obj.shapek);
thiscdf(InBounds) = exp(-XX);
% for i=1:numel(X)
% if X(i) <= obj.LowerBound
% elseif X(i) >= obj.UpperBound
% thiscdf(i) = 1;
% else
% XX = (X(i) - obj.mu) / obj.scale;
% XX = XX^(-obj.shapek);
% thiscdf(i) = exp(-XX);
% end
% end
end
function thisval=InverseCDF(obj,P)
[thisval, ~, Done] = MaybeSplineInvCDF(obj,P);
if Done
return;
end
XX = -log(P);
XX = XX.^(-1/obj.shapek);
thisval = XX * obj.scale + obj.mu;
end
end % methods
end % class ExtrVal2