Skip to content

Commit 6ad8d69

Browse files
[autofix.ci] apply automated fixes
1 parent 75eece2 commit 6ad8d69

10 files changed

Lines changed: 192 additions & 174 deletions

File tree

experiments/src/index.scala

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -8,4 +8,5 @@ import io.github.quafadas.plots.SetupVegaBrowser.{*, given}
88
val calYrIdx = vecxt_re.CalendarYearIndex(2025, idx.year, idx.idx)
99
println(calYrIdx)
1010
calYrIdx.plotIndex(1.0)
11-
println("finished")
11+
println("finished")
12+
end plotIndex

vecxt_re/src-jvm/dist/Empirical.scala

Lines changed: 21 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -6,14 +6,13 @@ import io.circe.syntax.*
66

77
/** Empirical distribution (JVM only).
88
*
9-
* This is a nonparametric distribution built directly from observed samples.
10-
* It supports positive weights $w_i$.
9+
* This is a nonparametric distribution built directly from observed samples. It supports positive weights $w_i$.
1110
*
12-
* The distribution is represented as a discrete measure on the (possibly repeated) sample values:
13-
* $$\mathbb{P}(X = x) = \sum_{i: x_i = x} \frac{w_i}{\sum_k w_k}.$$
11+
* The distribution is represented as a discrete measure on the (possibly repeated) sample values: $$\mathbb{P}(X = x) =
12+
* \sum_{i: x_i = x} \frac{w_i}{\sum_k w_k}.$$
1413
*
15-
* Consequently, the CDF is a right-continuous step function
16-
* $$F(t) = \mathbb{P}(X \le t) = \sum_{x \le t} \mathbb{P}(X=x).$$
14+
* Consequently, the CDF is a right-continuous step function $$F(t) = \mathbb{P}(X \le t) = \sum_{x \le t}
15+
* \mathbb{P}(X=x).$$
1716
*
1817
* Sampling is performed by inverse-transform sampling on the cumulative weights.
1918
*
@@ -44,6 +43,7 @@ case class Empirical(values: IArray[Double], weights: IArray[Double])
4443
j += 1
4544
end while
4645
out
46+
end pairs
4747

4848
scala.util.Sorting.stableSort(pairs, (a: (Double, Double), b: (Double, Double)) => a._1 < b._1)
4949

@@ -79,6 +79,7 @@ case class Empirical(values: IArray[Double], weights: IArray[Double])
7979
j += 1
8080
end while
8181
out
82+
end cdfVals
8283

8384
private val meanVal: Double =
8485
var s = 0.0
@@ -88,6 +89,7 @@ case class Empirical(values: IArray[Double], weights: IArray[Double])
8889
j += 1
8990
end while
9091
s
92+
end meanVal
9193

9294
private val varVal: Double =
9395
var s2 = 0.0
@@ -98,6 +100,7 @@ case class Empirical(values: IArray[Double], weights: IArray[Double])
98100
j += 1
99101
end while
100102
s2
103+
end varVal
101104

102105
def mean: Double = meanVal
103106

@@ -107,11 +110,14 @@ case class Empirical(values: IArray[Double], weights: IArray[Double])
107110
def probabilityOf(x: Double): Double =
108111
val idx = java.util.Arrays.binarySearch(xs, x)
109112
if idx >= 0 then probs(idx) else 0.0
113+
end if
114+
end probabilityOf
110115

111116
/** Draw a sample using inverse CDF sampling over the atomic masses. */
112117
def draw: Double =
113118
val u = rng.nextDouble()
114119
inverseCdf(u)
120+
end draw
115121

116122
/** CDF $F(t)=P(X\le t)$ (right-continuous). */
117123
def cdf(x: Double): Double =
@@ -136,15 +142,18 @@ case class Empirical(values: IArray[Double], weights: IArray[Double])
136142
val ip = java.util.Arrays.binarySearch(cdfVals, p)
137143
val idx = if ip >= 0 then ip else -ip - 1
138144
xs(math.min(idx, xs.length - 1))
145+
end if
146+
end inverseCdf
139147

140148
/** Plot a (weighted) histogram density estimate. */
141149
def plot(using viz.LowPriorityPlotTarget) =
142150
val plot = VegaPlot.fromResource("empiricalPdf.vl.json")
143-
val data = (0 until n).map { i => (x = values(i), w = weights(i)) }
151+
val data = (0 until n).map(i => (x = values(i), w = weights(i)))
144152
plot.plot(
145153
_.data.values := data.asJson,
146154
_ += (title = s"Empirical Distribution (n=$n)").asJson
147155
)
156+
end plot
148157

149158
/** Plot the empirical CDF (step function). */
150159
def plotCdf(using viz.LowPriorityPlotTarget) =
@@ -160,19 +169,21 @@ case class Empirical(values: IArray[Double], weights: IArray[Double])
160169
j += 1
161170
end while
162171
pts.toVector
172+
end points
163173

164174
plot.plot(
165175
_.data.values := points.asJson,
166176
_ += (title = s"Empirical CDF (n=$n)").asJson
167177
)
178+
end plotCdf
168179

169180
end Empirical
170181

171182
object Empirical:
172183
/** Construct an unweighted empirical distribution (all weights equal to $1$).
173184
*
174-
* Note: We intentionally avoid `apply` overloads here because `IArray[Double]` erases to `Array[Double]`
175-
* on the JVM, which can create signature collisions with the case class companion methods.
185+
* Note: We intentionally avoid `apply` overloads here because `IArray[Double]` erases to `Array[Double]` on the JVM,
186+
* which can create signature collisions with the case class companion methods.
176187
*/
177188
inline def equalWeights(values: Array[Double]): Empirical =
178189
Empirical(
@@ -183,4 +194,4 @@ object Empirical:
183194
/** Construct a weighted empirical distribution from arrays. */
184195
inline def weighted(values: Array[Double], weights: Array[Double]): Empirical =
185196
Empirical(IArray.from(values), IArray.from(weights))
186-
197+
end Empirical

vecxt_re/src-jvm/dist/NegativeBinomial.scala

Lines changed: 38 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -124,9 +124,14 @@ object NegativeBinomial:
124124
* @param tol
125125
* convergence tolerance for parameter 'a'
126126
* @return
127-
* Named tuple with `dist`: the fitted NegativeBinomial distribution, and `converged`: whether the optimizer converged within maxIter
127+
* Named tuple with `dist`: the fitted NegativeBinomial distribution, and `converged`: whether the optimizer
128+
* converged within maxIter
128129
*/
129-
def mle(observations: Array[Int], maxIter: Int = 500, tol: Double = 1e-8): (dist: NegativeBinomial, converged: Boolean) =
130+
def mle(
131+
observations: Array[Int],
132+
maxIter: Int = 500,
133+
tol: Double = 1e-8
134+
): (dist: NegativeBinomial, converged: Boolean) =
130135
require(observations.nonEmpty, "observations must not be empty")
131136
require(observations.forall(_ >= 0), "all observations must be non-negative")
132137

@@ -209,6 +214,7 @@ object NegativeBinomial:
209214

210215
if aNew <= 0 || !aNew.isFinite then a = a / 2.0
211216
else a = aNew
217+
end if
212218

213219
converged = math.abs(step * delta) < tol * math.abs(a)
214220
iter += 1
@@ -221,25 +227,21 @@ object NegativeBinomial:
221227

222228
/** Maximum likelihood estimation for the volume-adjusted Negative Binomial.
223229
*
224-
* We observe pairs $(n_j, v_j)$ where $n_j$ is the count and $v_j$ is the volume ratio (historical volume / modeled volume).
225-
* With parameters $(r, \beta)$ and $p = 1/(1+\beta v_j)$ the likelihood is
226-
* $$
227-
* L(r,\beta) = \prod_j \frac{\Gamma(r+n_j)}{\Gamma(r)\,\Gamma(n_j+1)} \left(\frac{\beta v_j}{1+\beta v_j}\right)^{n_j} \left(\frac{1}{1+\beta v_j}\right)^r.
228-
* $$
229-
* The log-likelihood is
230-
* $$
231-
* \ell(r,\beta) = \sum_j \big[\log\Gamma(r+n_j) - \log\Gamma(r) - \log\Gamma(n_j+1) + n_j(\log(\beta v_j) - \log(1+\beta v_j)) - r\,\log(1+\beta v_j)\big].
232-
* $$
233-
* Gradient components:
234-
* $$\partial_\beta \ell = \sum_j \Big( \frac{n_j}{\beta(1+\beta v_j)} - \frac{r v_j}{1+\beta v_j} \Big),\quad
235-
* \partial_r \ell = \sum_j \big[\psi(r+n_j) - \psi(r) - \log(1+\beta v_j)\big],$$
236-
* and Hessian entries:
237-
* $$\partial^2_{\beta\beta} \ell = \sum_j \Big( \frac{r v_j}{(1+\beta v_j)^2} - \frac{n_j(1+2\beta v_j)}{\beta^2(1+\beta v_j)^2} \Big),$$
238-
* $$\partial^2_{rr} \ell = \sum_j \big[\psi'(r+n_j) - \psi'(r)\big],\quad \partial^2_{\beta r} \ell = -\sum_j \frac{v_j}{1+\beta v_j}.$$
230+
* We observe pairs $(n_j, v_j)$ where $n_j$ is the count and $v_j$ is the volume ratio (historical volume / modeled
231+
* volume). With parameters $(r, \beta)$ and $p = 1/(1+\beta v_j)$ the likelihood is $$ L(r,\beta) = \prod_j
232+
* \frac{\Gamma(r+n_j)}{\Gamma(r)\,\Gamma(n_j+1)} \left(\frac{\beta v_j}{1+\beta v_j}\right)^{n_j}
233+
* \left(\frac{1}{1+\beta v_j}\right)^r. $$ The log-likelihood is $$ \ell(r,\beta) = \sum_j \big[\log\Gamma(r+n_j) -
234+
* \log\Gamma(r) - \log\Gamma(n_j+1) + n_j(\log(\beta v_j) - \log(1+\beta v_j)) - r\,\log(1+\beta v_j)\big]. $$
235+
* Gradient components: $$\partial_\beta \ell = \sum_j \Big( \frac{n_j}{\beta(1+\beta v_j)} - \frac{r v_j}{1+\beta
236+
* v_j} \Big),\quad \partial_r \ell = \sum_j \big[\psi(r+n_j) - \psi(r) - \log(1+\beta v_j)\big],$$ and Hessian
237+
* entries: $$\partial^2_{\beta\beta} \ell = \sum_j \Big( \frac{r v_j}{(1+\beta v_j)^2} - \frac{n_j(1+2\beta
238+
* v_j)}{\beta^2(1+\beta v_j)^2} \Big),$$ $$\partial^2_{rr} \ell = \sum_j \big[\psi'(r+n_j) - \psi'(r)\big],\quad
239+
* \partial^2_{\beta r} \ell = -\sum_j \frac{v_j}{1+\beta v_j}.$$
239240
*
240241
* Implementation details:
241242
* - Initialize from method of moments on rates $n_j / v_j$; if underdispersed, start at a small $\beta$.
242-
* - Newton updates solve the $2\times2$ system from the gradient/Hessian; a tiny ridge is added to keep the Hessian invertible.
243+
* - Newton updates solve the $2\times2$ system from the gradient/Hessian; a tiny ridge is added to keep the
244+
* Hessian invertible.
243245
* - Step halving is applied to enforce positivity of $r$ and $\beta$.
244246
*
245247
* @param observations
@@ -253,7 +255,12 @@ object NegativeBinomial:
253255
* @return
254256
* tuple of fitted `NegativeBinomial(r, beta)` and a convergence flag
255257
*/
256-
def volweightedMle(observations: Array[Int], volumes: Array[Double], maxIter: Int = 500, tol: Double = 1e-8): (dist: NegativeBinomial, converged: Boolean) =
258+
def volweightedMle(
259+
observations: Array[Int],
260+
volumes: Array[Double],
261+
maxIter: Int = 500,
262+
tol: Double = 1e-8
263+
): (dist: NegativeBinomial, converged: Boolean) =
257264
require(observations.nonEmpty, "observations must not be empty")
258265
require(observations.length == volumes.length, "observations and volumes must have the same length")
259266
require(observations.forall(_ >= 0), "all observations must be non-negative")
@@ -324,8 +331,7 @@ object NegativeBinomial:
324331
val hrrAdj = hrr + ridge
325332
val det = hbbAdj * hrrAdj - hbr * hbr
326333

327-
if det.isNaN || det.isInfinite || math.abs(det) < 1e-18 then
328-
iter = maxIter
334+
if det.isNaN || det.isInfinite || math.abs(det) < 1e-18 then iter = maxIter
329335
else
330336
val deltaBeta = (gBeta * hrrAdj - gR * hbr) / det
331337
val deltaR = (hbbAdj * gR - hbr * gBeta) / det
@@ -343,17 +349,23 @@ object NegativeBinomial:
343349
if newBeta > 0 && newR > 0 && newBeta.isFinite && newR.isFinite then
344350
beta = newBeta
345351
r = newR
346-
converged =
347-
math.abs(step * deltaBeta) <= tol * math.abs(beta) &&
348-
math.abs(step * deltaR) <= tol * math.abs(r)
352+
converged = math.abs(step * deltaBeta) <= tol * math.abs(beta) &&
353+
math.abs(step * deltaR) <= tol * math.abs(r)
349354
else iter = maxIter
355+
end if
350356
end if
351357

352358
iter += 1
353359
end while
354360

355361
(NegativeBinomial(r, beta), converged)
356-
357-
inline def mleVolumeWeighted(observations: Array[Int], volumes: Array[Double], maxIter: Int = 100, tol: Double = 1e-8): (dist: NegativeBinomial, converged: Boolean) = volweightedMle(observations, volumes, maxIter, tol)
362+
end volweightedMle
363+
364+
inline def mleVolumeWeighted(
365+
observations: Array[Int],
366+
volumes: Array[Double],
367+
maxIter: Int = 100,
368+
tol: Double = 1e-8
369+
): (dist: NegativeBinomial, converged: Boolean) = volweightedMle(observations, volumes, maxIter, tol)
358370

359371
end NegativeBinomial

vecxt_re/src-jvm/dist/Pareto.scala

Lines changed: 3 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -66,9 +66,9 @@ case class Pareto(scale: Double, shape: Double)
6666
def variance: Double = distribution.getVariance()
6767

6868
private def guessMaxXForPlot = shape match
69-
case s if s > 2 => mean + 4 * math.sqrt(variance) // mean and variance are defined
70-
case s if s > 1 => mean + 20 * scale // no well defined variance
71-
case _ => scale * 10 // no well defined mean
69+
case s if s > 2 => mean + 4 * math.sqrt(variance) // mean and variance are defined
70+
case s if s > 1 => mean + 20 * scale // no well defined variance
71+
case _ => scale * 10 // no well defined mean
7272

7373
def plot(using viz.LowPriorityPlotTarget) =
7474

@@ -93,7 +93,6 @@ case class Pareto(scale: Double, shape: Double)
9393
_.layer._1.data.sequence.stop := maxX,
9494
_.layer._1.data.sequence.step := (maxX - scale) / 200,
9595
_.layer._1.transform.head.calculate := pdfExpr,
96-
9796
_ += (title = s"Pareto Distribution PDF (scale=$scale, shape=$shape)").asJson
9897
)
9998
end plot

vecxt_re/src/CalendarYearIndex.scala

Lines changed: 7 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -92,14 +92,13 @@ object CalendarYearIndex:
9292
def plotIndex(reportingThreshold: Double)(using viz.LowPriorityPlotTarget) =
9393
val linePlot2 = VegaPlot.fromResource("index.vl.json")
9494
val cumulative = idx.onLevel(Array.fill(idx.years.length)(1.0), idx.years)
95-
val factors = idx.years.zip(idx.indices).zip(cumulative).map {
96-
case ((year, index), cumulative) =>
97-
(
98-
year = year,
99-
index = index,
100-
missing = 1 / cumulative,
101-
threshold = idx.suggestedNewThreshold(reportingThreshold)
102-
)
95+
val factors = idx.years.zip(idx.indices).zip(cumulative).map { case ((year, index), cumulative) =>
96+
(
97+
year = year,
98+
index = index,
99+
missing = 1 / cumulative,
100+
threshold = idx.suggestedNewThreshold(reportingThreshold)
101+
)
103102
}
104103
linePlot2.plot(
105104
_.data.values := factors.asJson

0 commit comments

Comments
 (0)