diff --git a/Project.toml b/Project.toml
index 484756e..692090b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "BiweightStats"
 uuid = "5bf8a1e9-d5f8-4697-9608-80edd97af0ad"
 authors = ["Miles Lucas <mdlucas@hawaii.edu> and contributors"]
-version = "1.0.0"
+version = "1.0.1"
 
 [deps]
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
diff --git a/src/BiweightStats.jl b/src/BiweightStats.jl
index 007fa98..a3547f4 100644
--- a/src/BiweightStats.jl
+++ b/src/BiweightStats.jl
@@ -55,7 +55,7 @@ end
 Creates an iterator based on the biweight transform.[^1] This iterator will first filter all input data so that only finite values remain. Then, the iteration will progress using a custom state, which includes a flag to indicate whether the value is within the cutoff, which is `c` times the median-absolute-deviation (MAD). The MAD is based on the deviation from `M`, which will default to the median of `X` if `M` is `nothing`.
 
 !!! note "Advanced usage"
-    
+
     This transform iterator is used for the internal calculations in `BiweightStats.jl`, which is why it has a somewhat complicated iterator implementation.
 
 # Examples
@@ -106,6 +106,10 @@ julia> (d, u2, flag), _ = iterate(bt, 10)
 """
 function BiweightTransform(X; c=9, M=nothing)
     data = filter(isfinite, X)
+    # if no valid input return NaN
+    if isempty(data)
+        return BiweightTransform(NaN, NaN, NaN)
+    end
     if isnothing(M)
         med = median(data)
     else
@@ -312,7 +316,7 @@ Computes biweight midcovariance between the two vectors. If only one vector is p
 ```
 
 !!! warning
-    
+
     `NaN` and `Inf` cannot be removed in the covariance calculation, so if they are present the returned value will be `NaN`. To prevent this, consider imputing values for the non-finite data.
 
 # Examples
@@ -372,7 +376,7 @@ midcov(X; kwargs...) = midvar(X; kwargs...)
 Computes the variance-covariance matrix using the biweight midcovariance. By default, each column is a separate variable, so an `(M, N)` matrix with `dims=1` will create an `(N, N)` covariance matrix. If `dims=2`, though, each row will become a variable, leading to an `(M, M)` covariance matrix.
 
 !!! warning
-    
+
     `NaN` and `Inf` cannot be removed in the covariance calculation, so if they are present the returned value will be `NaN`. To prevent this, consider imputing values for the non-finite data.
 
 # Examples
diff --git a/test/runtests.jl b/test/runtests.jl
index a218933..ab0bf85 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -41,6 +41,10 @@ rng = StableRNG(1123)
             val = location([1, 2, 3, v, 2])
             @test val ≈ 2 atol = 1e-5
         end
+
+        # All NaN
+        val = location([NaN, NaN, NaN])
+        @test isnan(val)
     end
 
     @testset "scale" begin
@@ -81,6 +85,10 @@ rng = StableRNG(1123)
             val = midvar([1, 2, 3, v, 2])
             @test val ≈ 0.55472 atol = 1e-5
         end
+
+        # All NaN
+        val = midvar([NaN, NaN, NaN])
+        @test isnan(val)
     end
 
     @testset "midcov" begin