diff --git a/core/MarkovProcesses.jl b/core/MarkovProcesses.jl
index d971b8013d6ccaa4d7352ec195d0066b7737d3b7..f6e5d042bd38c50c2fb07a648d68b059825aa53b 100644
--- a/core/MarkovProcesses.jl
+++ b/core/MarkovProcesses.jl
@@ -28,7 +28,7 @@ export get_index, get_value, length_var, isaccepted
 
 # Model related methods
 export simulate, volatile_simulate
-export distribute_mean_value_lha, mean_value_lha, distribute_prob_accept_lha
+export distribute_mean_value_lha, mean_value_lha, distribute_prob_accept_lha, probability_var_value_lha
 export set_param!, set_x0!, set_time_bound!, set_observed_var!, observe_all!
 export get_param, get_x0, getproperty, get_proba_model, get_observed_var
 export isbounded, isaccepted, check_consistency
diff --git a/core/model.jl b/core/model.jl
index 7857f9387ad28a9ed1f6e8f2da35f3f0b684626c..431e528671fe1c4ee06c094cf39eaf572bc7a2e6 100644
--- a/core/model.jl
+++ b/core/model.jl
@@ -382,6 +382,21 @@ function mean_value_lha(sm::SynchronizedModel, sym_var::VariableAutomaton, nbr_s
     end
     return sum_val / nbr_sim
 end
+"""
+    `distribute_var_value_lha(sm::SynchronizedModel, nbr_sim::Int, value = 0, sym_var = :d)
+
+Compute the probability that the variable `sym_var` is equal to  `value`
+of a LHA over `nbr_sim` simulations of the model.
+"""
+function probability_var_value_lha(sm::SynchronizedModel, nbr_sim::Int; 
+                                   value::Float64 = 0.0, sym_var::VariableAutomaton = :d)
+    sum_val = @distributed (+) for i = 1:nbr_sim
+        S = volatile_simulate(sm)
+        S[sym_var] == value
+    end
+    return sum_val / nbr_sim
+end
+
 
 function distribute_prob_accept_lha(sm::SynchronizedModel, nbr_sim::Int)
     sum_val = @distributed (+) for i = 1:nbr_sim