diff --git a/Project.toml b/Project.toml
index 777f3c7..b2f6135 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,11 +8,13 @@ BoundaryValueDiffEq = "764a87c0-6b3e-53db-9096-fe964310641d"
 Convex = "f65535da-76fb-5f13-bab9-19810c17039a"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ReinforcementLearning = "158674fc-8238-5cab-b5ba-03dfc80d1318"
 SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13"
diff --git a/src/ObjectiveFunc/AsymptoticBound/CramerRao.jl b/src/ObjectiveFunc/AsymptoticBound/CramerRao.jl
index 39ca0a1..4b33c82 100644
--- a/src/ObjectiveFunc/AsymptoticBound/CramerRao.jl
+++ b/src/ObjectiveFunc/AsymptoticBound/CramerRao.jl
@@ -565,6 +565,49 @@ function FIM(p::Vector{R}, dp::Vector{Vector{R}}; eps=GLOBAL_EPS) where {R<:Real
     
 end
 
+"""
+
+    FI_Expt(y1, y2, dx; ftype=:norm)
+
+Calculate the classical Fisher information (CFI) based on the experiment data.
+- `y1`: Experimental data obtained at the truth value (x).
+- `y1`: Experimental data obtained at x+dx.
+- `dx`: A known small drift of the parameter.
+- `ftype`: The distribution the data follows. Options are: norm, gamma, rayleigh, and poisson.
+"""
+function FI_Expt(y1, y2, dx; ftype=:norm)
+    Fc = 0.0
+    if ftype == :norm
+        p1_norm = fit(Normal, y1)
+        p2_norm = fit(Normal, y2)
+        f_norm(x) = sqrt(pdf(p1_norm, x)*pdf(p2_norm, x))
+        fidelity, err = quadgk(f_norm, -Inf, Inf)
+        Fc = 8*(1-fidelity)/dx^2
+    elseif ftype == :gamma
+        p1_gamma = fit(Gamma, y1)
+        p2_gamma = fit(Gamma, y2)
+        f_gamma(x) = sqrt(pdf(p1_gamma, x)*pdf(p2_gamma, x))
+        fidelity, err = quadgk(f_gamma, 0.0, Inf)
+        Fc = 8*(1-fidelity)/dx^2
+    elseif ftype == :rayleigh
+        p1_rayl = fit(Rayleigh, y1)
+        p2_rayl = fit(Rayleigh, y2)
+        f_rayl(x) = sqrt(pdf(p1_rayl, x)*pdf(p2_rayl, x))
+        fidelity, err = quadgk(f_rayl, 0.0, Inf)
+        Fc = 8*(1-fidelity)/dx^2
+    elseif ftype == :poisson
+        p1_pois = pdf.(fit(Poisson, y1), range(0, maximum(y1), step=1))
+        p2_pois = pdf.(fit(Poisson, y2), range(0, maximum(y2), step=1))
+        p1_pois, p2_pois = p1_pois/sum(p1_pois), p2_pois/sum(p2_pois)
+        fidelity = sum([sqrt(p1_pois[i]*p2_pois[i]) for i in 1:length(p1_pois)])
+        Fc = 8*(1-fidelity)/dx^2
+    else
+        println("supported values for ftype are 'norm', 'poisson', 'gamma' and 'rayleigh'")
+    end
+    return Fc
+end
+
+
 #======================================================#
 ################# Gaussian States QFIM #################
 function Williamson_form(A::AbstractMatrix)
diff --git a/src/QuanEstimation.jl b/src/QuanEstimation.jl
index 625330f..0df94e7 100644
--- a/src/QuanEstimation.jl
+++ b/src/QuanEstimation.jl
@@ -12,6 +12,8 @@ using SCS
 using BoundaryValueDiffEq
 using Trapz
 using Interpolations
+using Distributions
+using QuadGK
 
 const pkgpath = @__DIR__