clear
clc
u=[];
MSE=[];
times_k=50;
for i=1:times_k
Data= load('data_dui.txt');
S=sort(Data,'descend');
SS=S(10:4000);
X_k=SS(randi([1 length(SS)]));
S_k=S(S>=X_k);
k=find(S==X_k);
k_max=max(k);
times=200;
k_GP_200=[];
for i=1:times
S_k_i=[];
for ii=1:k_max
X_k_k=S_k(randi(numel(S_k),1,1));
S_k_i=[S_k_i;X_k_k];
end
y = S_k_i(S_k_i >X_k) - X_k;
paramEsts = gpfit(y);
k_GP = paramEsts(1);
k_GP_200=[k_GP_200;k_GP];
end
y_0 = S_k(S_k >X_k) - X_k;
paramEsts_0 = gpfit(y_0);
K_GP_O= paramEsts_0(1);
SUM_ki=sum(k_GP_200);
MM=k_GP_200-K_GP_O;
MM_2=MM.*MM;
SUM_KIK=sum(MM_2);
bias=0.005.*SUM_ki-K_GP_O;
MSE_K=bias.*bias+1./199.*SUM_KIK;
u=[u;X_k];
MSE=[MSE;MSE_K];
end