diff --git a/parmorceauxMax-git/Dynamics_Model.m b/parmorceauxMax-git/Dynamics_Model.m index 9ba2a4d16eda2a7e0846af657202756b2758300b..57a8d2403918e226eab2fff17f50e4b172c8e2b7 100644 --- a/parmorceauxMax-git/Dynamics_Model.m +++ b/parmorceauxMax-git/Dynamics_Model.m @@ -54,97 +54,9 @@ classdef Dynamics_Model < matlab.mixin.Copyable %handle obj.sigma_estimation = sigma_estimation; obj.sigma_test = sigma_test; end - function process_new_frame(obj, Yt, Tt) - obj.T = [obj.T; Tt]; - obj.Y = [obj.Y; Yt(:,1:obj.num_dim_Y)]; - - % estimate alpha and adjust the window size - alpha = obj.estimate_alpha(); - sprintf('estimated alpha_variable: %f', alpha); - obj.current_alpha = alpha; - - % estimate X inside the window - obj.estimate_X(); - %sprintf('"normalised" log likelihood of X (or confidence): %f', obj.confidence(frame)); - if obj.display_screen ~= -1 - obj.plot_estimated_X(); - end - end - function alpha = estimate_alpha(obj) - N = size(obj.X,1); - - if N < 3 - alpha = obj.alpha_initial; - return - end - - first_index = 1; - nb_points = N; - - if nb_points < 3 % we need at least 2 intervals to compute the mean - alpha = obj.alpha_initial; - else - X_2_to_N = obj.X(first_index+1:N); - X_1_to_Nminus1 = obj.X(first_index:N-1); - T_2_to_N = obj.T(first_index+1:N); - T_1_to_Nminus1 = obj.T(first_index:N-1); - - delta_X = X_2_to_N - X_1_to_Nminus1; - delta_T = T_2_to_N - T_1_to_Nminus1; - - factors = delta_X ./ delta_T; - - alpha = mean(factors); - - if alpha <= 0.85 * obj.alpha_initial; %0.75 %0.015 0.75%0.85%0.9 - alpha = 0.85 * obj.alpha_initial; %0.75%0.015;%0.85%0.9 - end - if alpha >= 1.15 * obj.alpha_initial; %1.25%0.015 1.25%1.15%1.1 - alpha = 1.15 * obj.alpha_initial; %1.25%0.015;%1.15%1.1 - end - end - end - function estimate_X(obj) - %function to be minimised over new_X, with extra parameters Y, alpha and T - f = @(variable_X)(-obj.logP_x_knowing_y_paquet(variable_X,obj.current_alpha,obj.sigma_estimation)); - - N = size(obj.X,1); - - if N == 0 - X_ini = obj.initial_estimate_X; - else - X_ini = [obj.X; obj.X(N) + obj.current_alpha * (obj.T(N+1) - obj.T(N))]; - end - - lb = repmat(0,1,length(X_ini)); - ub = repmat(1,1,length(X_ini)); - A = []; - b = []; - Aeq = []; - beq = []; - [new_X_variable,loglikelihood] = fmincon(f,X_ini,A,b,Aeq,beq,lb,ub); - - obj.X = new_X_variable; - end function spread_X(obj) N = 25; -% if N ==1 -% obj.X = 0; -% obj.current_alpha = obj.alpha_initial; -% -% obj.plot_estimated_X() -% obj.compute_loglikelihood_dynamics(1:N) -% return -% elseif N == 2 -% obj.X = [0 1]; -% obj.current_alpha = 1; -% -% obj.plot_estimated_X() -% obj.compute_loglikelihood_dynamics(1:N) -% return -% end - %%% estimation initiale de X X_tmp = 0:1/(N-1):1; @@ -197,33 +109,6 @@ classdef Dynamics_Model < matlab.mixin.Copyable %handle if obj.current_alpha >= 1.15 * obj.alpha_initial; %1.25%0.015 1.25%1.15%1.1 obj.current_alpha = 1.15 * obj.alpha_initial; %1.25%0.015;%1.15%1.1 end - - - obj.plot_estimated_X() - obj.compute_loglikelihood_dynamics(1:N) - end - function loglikelihood = logP_x_knowing_y_paquet(obj, variable_X, alpha, sigma) - N_variable = size(variable_X,1); - - % if N == 1, we maximise P_cond_y(Y(1),X(1)) * P_x_Markov_simple(X(1),X(0)) avec X(0) virtuel - % if N == 2, we maximise P_cond_y(Y(2),X(2)) * P_x_Markov_simple(X(2),X(1)) * P_cond_y(Y(1),X(1)) with P_x_Markov_simple(X(2),X(1)) = 0 if delta_X is negative, and uniform proba otherwise - - loglikelihood = obj.logP_cond_y(obj.Y(1,:),variable_X(1)); - loglikelihood = loglikelihood + obj.logP_x_Markov(variable_X(1), variable_X(1)-alpha, 2, 1, alpha, sigma); - - if N_variable > 1 - Y_2_to_N = obj.Y(2:N_variable,:); - X_2_to_N = variable_X(2:N_variable); - X_1_to_Nminus1 = variable_X(1:N_variable-1); - T_2_to_N = obj.T(2:N_variable); - T_1_to_Nminus1 = obj.T(1:N_variable-1); - - ll_proba_cond_y = obj.logP_cond_y(Y_2_to_N, X_2_to_N); - ll_proba_x_Markov = obj.logP_x_Markov(X_2_to_N,X_1_to_Nminus1,T_2_to_N,T_1_to_Nminus1,alpha,sigma); - vect_tmp = ll_proba_cond_y + ll_proba_x_Markov; - - loglikelihood = loglikelihood + sum(vect_tmp); - end end function loglikelihood = logP_x_knowing_y_spread(obj, variable_X, alpha, sigma) N_variable = size(variable_X,1); @@ -281,7 +166,6 @@ classdef Dynamics_Model < matlab.mixin.Copyable %handle loglikelihood(index2) = -100;%a été supprimé end end - function llh_dynamics = compute_loglikelihood_dynamics(obj, frames) %llh_dynamics(frames <= 1) = nan; @@ -323,7 +207,6 @@ classdef Dynamics_Model < matlab.mixin.Copyable %handle end end - function [im]=final_plot_estimated_X(obj) if obj.display_screen ~= -1 N = size(obj.X,1); @@ -346,7 +229,6 @@ classdef Dynamics_Model < matlab.mixin.Copyable %handle %saveas(obj.display_screen, strcat('estimatedX_frame', int2str(N), '.png')) end end - function X = normalise_X(obj, X) X(X<0) = 0; X(X>1) = 1;