etalement des X

parmorceauxMax-git/Dynamics_Model.m

@@ -223,6 +223,68 @@ classdef Dynamics_Model < matlab.mixin.Copyable %handle
             obj.N_fixed = new_N_fixed;
             obj.X = new_X;
         end
+        function spread_X(obj)
+            N = size(obj.X,1);
+            
+            if N ==1
+                obj.X = 0;
+                obj.confidence = 1;
+                obj.current_alpha = obj.alpha_initial;
+            
+                obj.plot_estimated_X()
+                obj.compute_loglikelihood_dynamics_forced(1:N)
+                return
+            elseif N == 2
+                obj.X = [0 1];
+                obj.confidence = [0 1];
+                obj.current_alpha = 1;
+            
+                obj.plot_estimated_X()
+                obj.compute_loglikelihood_dynamics_forced(1:N)
+                return
+            end
+            
+            X_tmp = 0:1/(N-1):1;
+            X_2_to_N = X_tmp(2:N)';
+            X_1_to_Nminus1 = X_tmp(1:N-1)';
+            T_2_to_N = obj.T(2:N);
+            T_1_to_Nminus1 = obj.T(1: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);
+            
+
+            %function to be minimised over new_X, with extra parameters Y, alpha and T
+            f = @(variable_X)(-obj.logP_x_knowing_y_spread(variable_X,alpha,obj.sigma_estimation));
+
+            
+            X_ini = X_tmp(2:end-1)';
+            
+            [new_X,loglikelihood] = fminunc(f,X_ini);
+            
+            new_X = [0; obj.normalise_X(new_X); 1];
+            obj.X = new_X;
+            obj.confidence = [1, -loglikelihood / size(new_X,1), 1];
+            
+            
+            
+            X_2_to_N = obj.X(2:N);
+            X_1_to_Nminus1 = obj.X(1:N-1);
+            T_2_to_N = obj.T(2:N);
+            T_1_to_Nminus1 = obj.T(1: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;
+            obj.current_alpha = mean(factors);
+            
+            obj.plot_estimated_X()
+            obj.compute_loglikelihood_dynamics_forced(1:N)
+        end
         function loglikelihood = logP_x_knowing_y_partialUpdate(obj, last_fixed_X, variable_X, alpha, sigma)
             N_variable = size(variable_X,1);
 
@@ -249,6 +311,21 @@ classdef Dynamics_Model < matlab.mixin.Copyable %handle
                 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);
+
+            % we maximise P_cond_y(Y(1),0) * P_x_Markov_simple(X(2),0) * P_cond_y(Y(i),X(i)) * P_x_Markov_simple(X(i+1),X(i)) * P_x_Markov_simple(1,X(N-1)) * P_cond_y(Y(N),1)
+
+            loglikelihood = obj.logP_cond_y(obj.Y(1,:),0) + obj.logP_x_Markov(variable_X(1),0,obj.T(2),obj.T(1),alpha,sigma);
+            loglikelihood = loglikelihood + sum(obj.logP_cond_y(obj.Y(2:end-1,:), variable_X));
+            loglikelihood = loglikelihood + obj.logP_x_Markov(1, variable_X(end),obj.T(end),obj.T(end-1),alpha,sigma);
+            loglikelihood = loglikelihood + obj.logP_cond_y(obj.Y(end,:),1);
+            
+            if N_variable > 1
+                ll_proba_x_Markov = obj.logP_x_Markov(variable_X(2:end),variable_X(1:end-1),obj.T(3:end-1),obj.T(2:end-2),alpha,sigma);
+                loglikelihood = loglikelihood + sum(ll_proba_x_Markov);
+            end
+        end
         function loglikelihood = logP_cond_y(obj, Yi, Xi)
             s = size(Yi,1);
             comp_small = repmat(obj.limits_Yi(:,1)',s,1);
@@ -287,27 +364,6 @@ classdef Dynamics_Model < matlab.mixin.Copyable %handle
                 loglikelihood = log(normpdf(difference,0,sigma));
             end
         end
-        function llh_dynamics = compute_loglikelihood_dynamics_inwindow(obj)
-            N = size(obj.X,1);
-            start_frame = obj.N_fixed + 1;
-            nb_points = N - start_frame + 1;
-            
-            terme1 = sum( obj.logP_cond_y(obj.Y(start_frame:N,:), obj.X(start_frame:N)) );
-            
-            if nb_points == 1
-                llh_dynamics = terme1 + obj.logP_x_Markov(obj.alpha_initial, 0, 2, 1, obj.alpha_initial, obj.sigma_test);
-                return
-            end
-            
-            X_2_to_N = obj.X(start_frame+1:N);
-            X_1_to_Nminus1 = obj.X(start_frame:N-1);
-            T_2_to_N = obj.T(start_frame+1:N);
-            T_1_to_Nminus1 = obj.T(start_frame:N-1);
-            
-            
-            terme2 = sum( obj.logP_x_Markov(X_2_to_N,X_1_to_Nminus1,T_2_to_N,T_1_to_Nminus1,obj.current_alpha,obj.sigma_test) );
-            llh_dynamics = (terme1 + terme2) / nb_points;
-        end
         function list_frames = get_converged_frames(obj)
             if obj.previous_N_fixed < obj.N_fixed
                 list_frames = [obj.previous_N_fixed+1:obj.N_fixed];
@@ -315,27 +371,7 @@ classdef Dynamics_Model < matlab.mixin.Copyable %handle
                 list_frames = [];
             end
         end
-        function llh_dynamics = compute_loglikelihood_dynamics(obj, frames)
-            indexes_float = frames > obj.N_fixed;
-            if sum(indexes_float) ~= 0
-                sprintf('compute_loglikelihood_dynamics() warning: the following frames have not converged yet:')
-                frames(indexes_float)
-                
-                llh_dynamics(indexes_float) = nan;
-            end
-            
-            %llh_dynamics(frames <= 1) = nan;
-            
-            indexes = frames > 1 & frames <= obj.N_fixed;
-            
-            terme1 = obj.logP_cond_y(obj.Y(frames(indexes),:), obj.X(frames(indexes)));
-            terme2 = obj.logP_x_Markov(obj.X(frames(indexes)),obj.X(frames(indexes)-1),obj.T(frames(indexes)),obj.T(frames(indexes)-1),obj.alpha_variable(frames(indexes))',obj.sigma_test);
-            llh_dynamics(indexes) = terme1 + terme2;
-            
-            indexes = frames <= 1;
-            terme1 = obj.logP_cond_y(obj.Y(frames(indexes),:), obj.X(frames(indexes)));
-            llh_dynamics(frames <= 1) = terme1 + obj.logP_x_Markov(obj.alpha_initial, 0, 2, 1, obj.alpha_initial, obj.sigma_test);
-        end
+        
         function llh_dynamics = compute_loglikelihood_dynamics_forced(obj, frames)
             %llh_dynamics(frames <= 1) = nan;
             
@@ -355,7 +391,8 @@ classdef Dynamics_Model < matlab.mixin.Copyable %handle
             indexes = frames <= 1
             frames(indexes)
             terme1 = obj.logP_cond_y(obj.Y(frames(indexes),:), obj.X(frames(indexes)));
-            llh_dynamics(frames <= 1) = terme1 + obj.logP_x_Markov(obj.alpha_initial, 0, 2, 1, obj.alpha_initial, obj.sigma_test);
+            %llh_dynamics(frames <= 1) = terme1 + obj.logP_x_Markov(obj.alpha_initial, 0, 2, 1, obj.alpha_initial, obj.sigma_test);
+            llh_dynamics(frames <= 1) = terme1 + obj.logP_x_Markov(obj.X(1), -alpha, 2, 1, alpha, obj.sigma_test);
         end
         function plot_estimated_X(obj)
             if obj.display_screen ~= -1

parmorceauxMax-git/Interface.m

@@ -104,6 +104,8 @@ classdef Interface < handle
             obj.meilleure_sequence_models = [obj.morceaux_finis{indice} obj.dynamics_models{indice}];
             
             %% calcul les scores pour les frames du dernier morceau
+            % mise à jour des x pour bien les étaler
+            obj.meilleure_sequence_models(end).spread_X();
             new_scores = obj.meilleure_sequence_models(end).compute_loglikelihood_dynamics_forced([1:obj.N-obj.first_frame_of_models(indice)+1]);
             obj.loglikelihood = [obj.scores_finis{indice} new_scores];
             
@@ -165,7 +167,7 @@ classdef Interface < handle
                     new_first_frame_of_models(morceau) = obj.first_frame_of_models(morceau);
                     
                     new_models{morceau}.process_new_frame(Yt, Tt);
-                    proba_obs = new_models{morceau}.compute_loglikelihood_dynamics_inwindow();
+                    proba_obs = new_models{morceau}.compute_loglikelihood_dynamics_forced(1);
                     
                     new_proba_morceaux(morceau) = log(obj.proba_morceaux(morceau)) + proba_obs;
                     new_morceaux_finis{morceau} = [];
@@ -213,7 +215,7 @@ classdef Interface < handle
                         else
                             nouveau_morceau_fini_branche{branche} = copy(obj.dynamics_models{branche}); % ce modèle n'a pas encore vu la frame courante, on peut donc le sauvegader dans son état à la frame précédente
 %%%%                        %%%%%% mettre à jour la copie ici
-                            
+                            nouveau_morceau_fini_branche{branche}.spread_X();
                         end
                         
                         %% on récupère les scores des frames vues par le modèle