.

Author: Quafadas

vecxt/src/intarray.scala

@@ -5,6 +5,7 @@ import scala.util.control.Breaks.*
 
 import vecxt.BooleanArrays.trues
 import vecxt.BoundsCheck.BoundsCheck
+import vecxt.arrays.sumSIMD
 
 object IntArrays:
 
@@ -35,6 +36,25 @@ object IntArrays:
       out
     end select
 
+    inline def countsToIdx: Array[Int] =
+      var total = arr.sumSIMD
+      var i = 0
+      val out = new Array[Int](total)
+      var j = 0
+      while i < arr.length do
+        val count = arr(i)
+        val idx = i + 1
+        var k = 0
+        while k < count do
+          out(j) = idx
+          j += 1
+          k += 1
+        end while
+        i += 1
+      end while
+      out
+    end countsToIdx
+
     inline def contiguous: Boolean =
       var i = 1
       var out = true

vecxt/test/src/intarray.test.scala

@@ -133,4 +133,31 @@ class IntArrayExtensionSuite extends munit.FunSuite:
     assertVecEquals(base.select(idx), Array.emptyIntArray)
   }
 
+  test("countsToIdx basic") {
+    val counts = Array(2, 3, 1)
+    // 2 ones, 3 twos, 1 three => [1, 1, 2, 2, 2, 3]
+    assertVecEquals(counts.countsToIdx, Array(1, 1, 2, 2, 2, 3))
+  }
+
+  test("countsToIdx with zeros") {
+    val counts = Array(1, 0, 2, 0, 1)
+    // 1 one, 0 twos, 2 threes, 0 fours, 1 five => [1, 3, 3, 5]
+    assertVecEquals(counts.countsToIdx, Array(1, 3, 3, 5))
+  }
+
+  test("countsToIdx empty array") {
+    val counts = Array.emptyIntArray
+    assertVecEquals(counts.countsToIdx, Array.emptyIntArray)
+  }
+
+  test("countsToIdx all zeros") {
+    val counts = Array(0, 0, 0)
+    assertVecEquals(counts.countsToIdx, Array.emptyIntArray)
+  }
+
+  test("countsToIdx single element") {
+    val counts = Array(5)
+    assertVecEquals(counts.countsToIdx, Array(1, 1, 1, 1, 1))
+  }
+
 end IntArrayExtensionSuite

vecxt_re/resources/hillPlot.vl.json

@@ -0,0 +1,33 @@
+{
+  "$schema": "https://vega.github.io/schema/vega-lite/v6.json",
+  "title": "Hill Plot",
+  "width": "container",
+  "height": "container",
+  "data": {
+    "values": [
+      { "k": 10, "estimate": 2.1 },
+      { "k": 20, "estimate": 2.05 },
+      { "k": 30, "estimate": 1.98 },
+      { "k": 40, "estimate": 2.02 },
+      { "k": 50, "estimate": 1.95 }
+    ]
+  },
+  "mark": {
+    "type": "line",
+    "point": { "filled": true, "size": 40 },
+    "tooltip": { "content": "data" }
+  },
+  "encoding": {
+    "x": {
+      "field": "k",
+      "type": "quantitative",
+      "title": "k (number of upper order statistics)"
+    },
+    "y": {
+      "field": "estimate",
+      "type": "quantitative",
+      "title": "Tail Index Estimate α̂(k)"
+    },
+    "color": { "value": "steelblue" }
+  }
+}

vecxt_re/resources/loglogCdf.vl.json

@@ -0,0 +1,56 @@
+{
+  "$schema": "https://vega.github.io/schema/vega-lite/v6.json",
+  "title": "Log-Log CDF Plot",
+  "width": "container",
+  "height": "container",
+  "data": {
+    "values": [
+      { "x": 1, "y": 0.9, "source": "model" },
+      { "x": 10, "y": 0.5, "source": "model" },
+      { "x": 100, "y": 0.1, "source": "model" },
+      { "x": 1000, "y": 0.01, "source": "model" },
+      { "x": 2, "y": 0.85, "source": "empirical" },
+      { "x": 15, "y": 0.55, "source": "empirical" },
+      { "x": 80, "y": 0.15, "source": "empirical" },
+      { "x": 500, "y": 0.05, "source": "empirical" }
+    ]
+  },
+  "layer": [
+    {
+      "transform": [{ "filter": "datum.source === 'model'" }],
+      "mark": "line",
+      "encoding": {
+        "x": {
+          "field": "x",
+          "type": "quantitative",
+          "scale": { "type": "log", "domainMin": 1 },
+          "title": "X (log scale)"
+        },
+        "y": {
+          "field": "y",
+          "type": "quantitative",
+          "scale": { "type": "log" },
+          "title": "Survival Probability S(x)"
+        },
+        "color": { "value": "steelblue" }
+      }
+    },
+    {
+      "transform": [{ "filter": "datum.source === 'empirical'" }],
+      "mark": { "type": "point", "shape": "cross", "size": 100, "tooltip": {"content": "data"} },
+      "encoding": {
+        "x": {
+          "field": "x",
+          "type": "quantitative",
+          "scale": { "type": "log", "domainMin": 1 }
+        },
+        "y": {
+          "field": "y",
+          "type": "quantitative",
+          "scale": { "type": "log" }
+        },
+        "color": { "value": "red" }
+      }
+    }
+  ]
+}

vecxt_re/src-jvm/all.scala

@@ -3,6 +3,7 @@ package vecxt_re
 object all:
   export vecxt_re.Scenario
   export vecxt_re.Scenarr
+  export vecxt_re.Tower
   export vecxt_re.Plots.*
   export vecxt_re.TrendAnalysis.*
   export vecxt_re.TrendFitResult
@@ -16,4 +17,8 @@ object all:
   export vecxt_re.Poisson
   export vecxt_re.Empirical
   export vecxt_re.Pareto
+  export vecxt_re.Mixed
+  export vecxt_re.HillEstimator
+  export vecxt_re.PickandsEstimator
+
 end all

vecxt_re/src-jvm/plots.scala

@@ -2,6 +2,8 @@ package vecxt_re
 
 import io.circe.syntax.*
 import io.github.quafadas.plots.SetupVega.{*, given}
+import vecxt_re.HillEstimator.HillPlotResult
+import vecxt_re.PickandsEstimator.PickandsPlotResult
 
 object Plots:
   // These must be private otherwise scaladoc get crazy.
@@ -13,6 +15,8 @@ object Plots:
   private lazy val rootogram = VegaPlot.fromResource("rootogram.vl.json") // hanging rootogram
   private lazy val pearsonResiduals = VegaPlot.fromResource("pearsonResiduals.vl.json") // residual plot
   private lazy val poissonTrend = VegaPlot.fromResource("poissonTrend.vl.json") // Poisson GLM trend
+  private lazy val logLogPlot = VegaPlot.fromResource("loglogCdf.vl.json") // log-log CDF plot
+  private lazy val hillPlotSpec = VegaPlot.fromResource("hillPlot.vl.json") // Hill plot for tail index
 
   extension (idx: CalendarYearIndex)
     def plotIndex(reportingThreshold: Double)(using viz.LowPriorityPlotTarget) =
@@ -351,6 +355,108 @@ object Plots:
       )
   end extension
 
+  extension (mixed: Mixed)
+    /** Plot log-log comparison of theoretical Mixed distribution vs empirical sample data.
+      *
+      * This plot shows the complementary CDF (1 - CDF) on log-log scale, which is useful for visualizing tail behavior
+      * of heavy-tailed distributions. The x-axis shows log(value) and y-axis shows log(1 - CDF).
+      */
+    def plotLogLogVsSample(samples: IndexedSeq[Double], threshold: Double)(using viz.LowPriorityPlotTarget) =
+      val sortedSamples = samples.sorted
+      val n = sortedSamples.length.toDouble
+
+      // Empirical survival function using Hazen plotting position
+      val empiricalData = sortedSamples.zipWithIndex.collect {
+        case (x, i) if x > 0 =>
+          val survivalProb = (n - i - 0.5) / n
+          if survivalProb > 0 then Some((x = x, y = survivalProb, source = "empirical"))
+          else None
+          end if
+      }.flatten
+
+      // Theoretical survival function (1 - CDF)
+      // For Pareto: S(x) = (xₘ/x)^α, so log(S) = α·log(xₘ) - α·log(x) is linear
+      val minX = sortedSamples.filter(_ > 0).headOption.getOrElse(1.0)
+      val maxX = sortedSamples.last
+      val numPoints = 200
+      val theoreticalData = (0 until numPoints).flatMap { i =>
+        val x = minX * math.pow(maxX / minX, i.toDouble / (numPoints - 1))
+        val survivalProb = 1.0 - mixed.cdf(x)
+        if survivalProb > 1e-10 && x > 0 then Some((x = x, y = survivalProb, source = "model"))
+        else None
+        end if
+      }
+
+      val allData = theoreticalData ++ empiricalData
+
+      logLogPlot.plot(
+        _.title(s"Mixed Distribution Log-Log Plot"),
+        _.data.values := allData.asJson,
+        _.layer._0.encoding.x.scale.domainMin := threshold
+      )
+    end plotLogLogVsSample
+  end extension
+
+  extension (hp: HillPlotResult)
+    /** Plot a Hill plot showing tail index estimates α̂(k) vs k.
+      *
+      * A Hill plot helps identify the optimal number of upper order statistics to use for Pareto tail estimation. Look
+      * for a stable plateau region where the estimate is relatively constant - this indicates the range of k where the
+      * Pareto assumption holds. Too small k gives high variance; too large k includes non-tail observations.
+      */
+    def plotHill(using viz.LowPriorityPlotTarget) =
+      val data = hp.kValues.zip(hp.estimates).map { case (k, est) =>
+        (k = k, estimate = est)
+      }
+
+      hillPlotSpec.plot(
+        _.title("Hill Plot for Pareto Tail Index Estimation"),
+        _.data.values := data.asJson
+      )
+    end plotHill
+  end extension
+
+  extension (pp: PickandsPlotResult)
+    /** Plot a Pickands plot showing tail index estimates α̂(k) = 1/γ̂(k) vs k.
+      *
+      * The Pickands estimator is more robust than Hill to model misspecification but has higher variance. Look for a
+      * stable plateau region. Unlike the Hill plot, the Pickands estimator can give negative γ values for light-tailed
+      * distributions; only positive γ (and thus positive α) indicates heavy tails.
+      */
+    def plotPickands(using viz.LowPriorityPlotTarget) =
+      val data = pp.kValues
+        .zip(pp.alphaEstimates)
+        .filter { case (_, est) => !est.isNaN && est.isFinite && est > 0 }
+        .map { case (k, est) =>
+          (k = k, estimate = est)
+        }
+
+      hillPlotSpec.plot(
+        _.title("Pickands Plot for Pareto Tail Index Estimation"),
+        _.data.values := data.asJson
+      )
+    end plotPickands
+
+    /** Plot the raw extreme value index γ̂(k) from Pickands estimator.
+      *
+      * For heavy-tailed data, γ > 0. For light-tailed data (e.g., exponential), γ = 0. For bounded distributions, γ <
+      * 0.
+      */
+    def plotPickandsGamma(using viz.LowPriorityPlotTarget) =
+      val data = pp.kValues
+        .zip(pp.gammaEstimates)
+        .filter { case (_, est) => !est.isNaN && est.isFinite }
+        .map { case (k, est) =>
+          (k = k, estimate = est)
+        }
+
+      hillPlotSpec.plot(
+        _.title("Pickands Plot for Extreme Value Index γ"),
+        _.data.values := data.asJson
+      )
+    end plotPickandsGamma
+  end extension
+
   // extension (negBin: NegativeBinomial)
   //   inline def plotPdf(using viz.LowPriorityPlotTarget) =
   //     val numPoints = 1000