fix pow for negative value

src/Dist/logLogistic.jl

@@ -26,7 +26,7 @@ end
 
 function cdf_logLogistic(x::Real, params)
   β, α, γ = params[1:3]
-  1 / (1 + ((α / (x - β))^γ))
+  1 / (1 + (pow(α / (x - β), γ)))
 end
 
 

src/SPEI.jl

@@ -8,7 +8,14 @@ using SpecialFunctions: gamma_inc, lgamma, loggamma
 
 
 qnorm(p::Real) = quantile(Normal(), p)
-pow = ^
+
+export pow
+# pow = ^
+# import Base:^
+pow(x, y) = x^y
+function pow(x::T, y::AbstractFloat) where {T<:AbstractFloat}
+  x < 0 ? T(NaN) : x^y
+end
 
 
 include("PWM.jl")

src/main_spei.jl

@@ -1,6 +1,8 @@
 function param_spei(x_ref::AbstractVector)
   x2 = x_ref[.!isnan.(x_ref)] |> sort # have to sort
   beta = pwm(x2, 0.0, 0.0, 0)
+  @show beta
+
   params = _fit_logLogistic(beta)
   params
 end
@@ -21,6 +23,8 @@ function param_spi(x_ref::AbstractVector; fit="lmom")
   params(D), q
 end
 
+export param_spei, param_spi
+
 
 """
     spei(x::AbstractVector)

test/test-spei.jl

@@ -2,24 +2,27 @@ using Test
 using RCall
 using SPEI
 
+R"""
+library(SPEI2)
+"""
+
 function nanmax(x)
   x2 = x[.!isnan.(x)]
   maximum(x2)
 end
 
-R"""
-library(SPEI2)
-"""
-wb = R"wb" |> rcopy |> x -> Float32.(x)
+# wb = R"wb" |> rcopy |> x -> Float32.(x)
+wb = deserialize("data/wb")
+x, x_ref = deserialize("data/x_ref")
 
 @testset "spei" begin
-  @time r_R = R"cal_spei(wb)" |> rcopy
-  @time r_jl = spei(wb)
+  @time r_R = R"cal_spei($x, $x_ref)$z" |> rcopy
+  @time r_jl = spei(x, x_ref).z
 
-  z_r = r_R[:z]
-  z_jl = r_jl.z
-  e = z_r - z_jl
-  @test maximum(abs.(e)) <= 1e-3
+  # z_r = r_R[:z]
+  # z_jl = r_jl.z
+  e = r_R - r_jl
+  @test nanmax(e) <= 1e-3
 end
 
 @testset "spi" begin
@@ -33,11 +36,8 @@ end
 end
 
 # > 计算结果完全一致
-
-# ## 误差还挺大的
 # begin
 #   using Plots
-
 #   plot(z_r, label="R")
 #   plot!(z_jl, label="Julia")
 #   plot!(e, label="R - Julia")