search

sbegueria / SPEI

An R package for computing the Standardized Precipitation-Evapotranspiration Index (SPEI), Penman-Monteith and other reference evapotranspiration, SPI, etc.
78 stars 34 forks source link

Using max-lik fit returns all NA #67

Open ThomasOtt314 opened 1 year ago

ThomasOtt314 commented 1 year ago

Excecuting the code below returns a time series of NA for the SPI calulation using a fit of 'max-lik' The code works fine with the 'ub-pwm' fit. Any help would be geat! 

library(SPEI)

date1 <- '1991-01-01'
date2 <- '2022-12-31'

# Vector of precip data
prcp <- c(12.385, 14.56, 84.177, 82.407, 140.11, 71.865, 122.404, 19.689, 74.63, 106.985, 104.498, 29.507,
          25.848, 23.135, 32.936, 52.319, 61.574, 51.701, 70.253, 82.455, 102.266, 38.425, 118.912, 56.103,
          35.346, 1.337, 8.652, 76.22, 124.201, 137.701, 55.681, 87.971, 85.293, 48.03, 43.44, 11.535,
          24.986, 15.565, 38.001, 55.143, 51.189, 67.051, 127.859, 62.913, 147.208, 52.577, 45.874, 3.914,
          17.281, 8.121, 49.453, 57.111, 95.307, 41.837, 83.787, 173.801, 69.909, 124.574, 40.477, 56.383,
          92.18, 31.096, 45.262, 100.924, 30.166, 129.566, 151.449, 48.891, 101.847, 78.021, 25.922, 84.927,
          70.358, 6.09, 51.512, 18.66, 67.187, 94.862, 57.436, 100.444, 58.134, 92.501, 9.297, 17.738, 63.015,
          33.73, 83.186, 40.893, 31.306, 85.689, 43.181, 68.625, 38.336, 56.033, 56.388, 25.188, 67.26,
          46.581, 5.931, 42.391, 157.64, 76.645, 274.669, 73.131, 54.014, 52.27, 25.362, 23.102, 41.704,
          34.995, 53.772, 59.42, 50.702, 114.356, 145.441, 87.555, 138.506, 19.834, 56.441, 28.874, 31.283,
          33.95, 15.962, 72.31, 86.051, 96.976, 55.149, 42.773, 86.452, 63.163, 55.593, 34.201, 6.104, 44.355,
          72.547, 140.646, 73.256, 80.686, 87.629, 91.156, 78.412, 182.947, 11.309, 10.609, 7.404, 37.726,
          78.151, 74.319, 72.562, 61.255, 26.069, 134.898, 97.566, 31.052, 48.441, 39.8, 24.741, 48.004,
          96.378, 53.094, 100.726, 81.675, 103.615, 70.379, 41.732, 114.762, 29.506, 47.318, 28.767, 27.808,
          18.898, 34.903, 77.017, 78.357, 76.733, 55.887, 69.255, 149.134, 71.69, 29.588, 34.341, 34.929,
          35.232, 19.747, 145.687, 31.257, 82.48, 119.92, 65.942, 48.313, 45.942, 48.417, 23.48, 18.238,
          42.067, 42.48, 40.272, 65.081, 53.738, 38.337, 76.987, 147.274, 15.899, 61.78, 39.436, 30.22, 
          21.472, 123.231, 73.043, 106.254, 62.813, 15.926, 57.42, 49.503, 30.557, 54.914, 15.794, 29.734, 
          18.37, 111.956, 84.803, 55.557, 46.68, 76.849, 31.811, 116.134, 16.584, 55.444, 13.689, 10.121,
          15.546, 10.445, 34.91, 178.911, 196.371, 68.862, 190.498, 64.618, 54.982, 32.544, 14.395, 3.747,
          42.236, 86.608, 99.329, 94.798, 46.932, 33.443, 89.661, 40.399, 30.978, 36.376, 31.049, 27.903,
          35.629, 31.926, 115.1832, 49.418, 79.482, 34.024, 53.833, 85.008, 34.768, 35.674, 32.502, 36.513,
          42.602, 99.468, 79.037, 90.192, 101.719, 77.357, 51.073, 66.202, 69.336, 33.474, 34.997, 34.114,
          31.174, 87.496, 69.826, 88.95, 38.837, 108.239, 164.872, 84.869, 71.43, 53.455, 15.693, 6.145,
          18.74, 79.687, 125.433, 89.968, 65.675, 91.385, 78.546, 74.637, 67.385, 102.195, 22.46, 22.227,
          73.878, 65.999, 55.27, 88.588, 93.457, 84.378, 109.766, 115.427, 43.971, 50.174, 36.989, 37.836,
          27.074, 158.36, 120.466, 177.778, 68.125, 102.783, 58.884, 84.327, 34.88, 47.297, 28.806, 30.719,
          16.388, 61.368, 59.829, 104.492, 70.354, 88.753, 120.826, 173.388, 55.717, 32.44, 30.041, 108.874,
          28.476, 78.3, 128.971, 40.769, 100.428, 67.155, 124.736, 108.245, 31.703, 107.648, 47.778, 11.688,
          75.845, 68.2, 55.504, 130.53, 195.374, 105.965, 94.044, 67.995, 58.563, 15.757, 15.1248, 30.168,
          45.9477, 76.3019, 41.0862, 106.5615, 84.8749, 64.9139, 74.4176, 18.4552, 24.8467, 65.7223, 19.3086,
          20.6611, 69.5475, 79.5091, 117.4526, 87.4791, 39.3405, 110.511, 37.1871, 28.1281, 68.0051, 39.1535)

# Build the date frame of the data
dates <- seq(as.Date(date1, "%Y-%m-%d"), as.Date(date2, "%Y-%m-%d"), by = "month")
months <- as.numeric(format(dates, "%m"))
years <- as.numeric(format(dates, "%Y"))
data <- data.frame(years, months, prcp)

# Compute 1-month SPI
spi1_ub_pwm <- spi(prcp, 1, fit='ub-pwm')
spi1_max_lik <- spi(prcp, 1, fit="max-lik")

print(spi1_ub_pwm)
print(spi1_max_lik)
ThomasOtt314 commented 1 year ago

Forgot to add versions. R version 4.3.1 (2023-06-16 ucrt) Package SPEI (1.8.1)

phdsdc commented 1 year ago

i found the same problem

xiaoqi1010 commented 1 year ago

anyone may solve this problem??