-
Notifications
You must be signed in to change notification settings - Fork 0
/
gumbel.py
111 lines (96 loc) · 4.44 KB
/
gumbel.py
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
import math
import numpy as np
def randomDataGenerator(alpha, beta, n):
return np.random.gumbel(alpha, beta, n)
def firstDerivative(alpha, beta, all_points, n, flag):
if (flag == 'alpha'):
first_term = 0
for datapoint in all_points:
first_term += math.exp(-((datapoint - alpha) / beta))
final_result = n / beta - first_term / beta
if (flag == 'beta'):
first_term = 0
second_term = 0
for datapoint in all_points:
first_term += ((datapoint - alpha) / (beta * beta))
second_term += ((datapoint - alpha) / (beta * beta)) * (math.exp(-((datapoint - alpha) / beta)))
final_result = -n / beta + first_term - second_term
return final_result
def secondDerivative(alpha, beta, all_points, n, flag):
if (flag == 'alpha'):
first_term = 0
for datapoint in all_points:
first_term += math.exp(-((datapoint - alpha) / beta))
final_result = (-1 / (beta * beta)) * first_term
if (flag == 'beta'):
first_term = 0
second_term = 0
third_term = 0
for datapoint in all_points:
first_term += (datapoint - alpha)
second_term += ((datapoint - alpha) * math.exp(-((datapoint - alpha) / beta)))
third_term += ((datapoint - alpha) * (datapoint - alpha) * math.exp(-((datapoint - alpha) / beta)))
first_term = (-2 / (beta * beta * beta)) * first_term
second_term = (2 / beta * beta * beta) * second_term
third_term = -third_term / (beta * beta * beta * beta)
final_result = (n / beta * beta) + first_term + second_term + third_term
return final_result
def secondDerivativeAlphaBeta(alpha, beta, all_points, n):
first_term = 0
second_term = 0
for datapoint in all_points:
first_term += math.exp(-((datapoint - alpha) / beta))
second_term += ((datapoint - alpha) * (math.exp(-(datapoint - alpha) / beta)))
first_term = first_term / (beta * beta)
second_term = -second_term / (beta * beta * beta)
final_result = (-n / (beta * beta)) + first_term + second_term
return final_result
def estimateGumbelParameters(n):
alpha_list = []
beta_list = []
for i in range(0, 10):
dataset = randomDataGenerator(0, 0.1, n)
mean_dataset = np.mean(dataset)
std_dataset = np.std(dataset)
beta = 0.7797 * std_dataset
alpha = mean_dataset - 0.5772 * beta
itr=0
while True:
itr+=1
first_alpha_derivative = firstDerivative(alpha, beta, dataset, n, 'alpha')
first_beta_derivative = firstDerivative(alpha, beta, dataset, n, 'beta')
second_alpha_derivative = secondDerivative(alpha, beta, dataset, n, 'alpha')
second_beta_derivative = secondDerivative(alpha, beta, dataset, n, 'beta')
second_alpha_beta_derivative = secondDerivativeAlphaBeta(alpha, beta, dataset, n)
gumbel_function = np.array([[first_alpha_derivative], [first_beta_derivative]])
hessian = np.array([[second_alpha_derivative, second_alpha_beta_derivative],
[second_alpha_beta_derivative, second_beta_derivative]])
hessian_inv = np.linalg.pinv(hessian)
out_old = np.array([[alpha], [beta]])
out_new = np.subtract(out_old, (hessian_inv * gumbel_function))
new_alpha = out_new[0][0]
new_beta = out_new[1][0]
threshold = ((new_alpha - alpha) ** 2) + ((new_beta - beta) ** 2)
alpha = new_alpha
beta = new_beta
if threshold < 10 ** -6:
print(itr)
break
if itr > 10:
break
alpha_list.append(alpha)
beta_list.append(beta)
return alpha_list, beta_list
def reportValues(n):
alpha_list, beta_list = estimateGumbelParameters(n)
print("--------------ALPHA VALUES FOR N = ",n," -------------------------------------------------------------------")
print("List of Alphas", alpha_list)
print("Mean value of alpha ", np.mean(alpha_list))
print("Standard Deviation alpha ", np.std(alpha_list))
print("--------------BETA VALUES FOR N = ",n," --------------------------------------------------------------------")
print("List of Betas", beta_list)
print("Mean value of beta ", np.mean(beta_list))
print("Standard Deviation beta ", np.std(beta_list))
temp = [100, 1000, 10000]
for entry in temp:
reportValues(entry)