ajwheeler / ajwheeler/Korg.jl

`ews_to_stellar_parameters`: `sigma_mean` is computed from `length` of a scalar, so no √N reduction is applied

Open
#571 1 comment 0 reactions 0 assignees View on GitHub
Dominant language
Julia
Stars
67
Forks
15
Avg merge
2d 6h
Merged PRs (30d)
6

Description

Hi, I found an issue in `src/Fit/fit_via_EWs.jl:560-567`:

```julia
estimated_err = std(A[isfinite.(A)])

sigma_mean = estimated_err ./ sqrt(length(estimated_err))
teff_residual_sigma = estimated_err *
get_slope_uncertainty([line.E_lower for line in linelist[neutrals]])
vmic_residual_sigma = estimated_err * get_slope_uncertainty(REWs[neutrals])

[teff_residual_sigma, sigma_mean, vmic_residual_sigma, sigma_mean]
```

`estimated_err` is a scalar, so `length(estimated_err)` is 1 and `sigma_mean`
is just `estimated_err`. Looks like the standard error of a mean was intended.

Reproduction, no Korg call needed:

```julia
using Statistics
A = randn(150) .+ 7.5
estimated_err = std(A[isfinite.(A)])
sigma_mean = estimated_err ./ sqrt(length(estimated_err))

sigma_mean == estimated_err # true
estimated_err / sqrt(length(A)) # 0.0775 vs 0.9493 — factor 12.25
```

`sigma_mean` is returned for equations 2 and 4, both differences of means:

- eq. 2, `mean(A_I) - mean(A_II)`: should be `err * sqrt(1/n_I + 1/n_II)`,
inflated ~4.2x for 150 Fe I + 20 Fe II
- eq. 4, `mean(A) - ([m/H] + A_solar)`: should be `err / sqrt(n_total)`,
inflated ~12x for 150 lines

Both feed `param_uncertainties` at line 426, so reported σ(log g) and σ([m/H])
are too large by those factors. Conservative, not dangerous — but it also makes
them incomparable with σ(Teff) and σ(vmic), which are scaled correctly through
`get_slope_uncertainty`.

The two equations have different denominators, so one `sigma_mean` cannot cover
both:

```julia
n_neutral = count(neutrals .& isfinite.(A))
n_ionised = count((.!neutrals) .& isfinite.(A))
n_total = count(isfinite.(A))

sigma_ionisation = estimated_err * sqrt(1/n_neutral + 1/n_ionised)
sigma_metallicity = estimated_err / sqrt(n_total)

[teff_residual_sigma, sigma_ionisation, vmic_residual_sigma, sigma_metallicity]
```

I can PR this. It shrinks reported σ(log g) and σ([m/H]) by roughly an order of
magnitude, so it may warrant a release note despite being a bug fix.

Contributor guide

No contributing guide indexed for this repository

Research direction

Start at src/Fit/fit_via_EWs.jl:560-567 and inspect how the returned uncertainty values are consumed at line 426. Reproduce the scalar-length behavior with the Julia example in the issue, then verify the ionisation and metallicity terms use their respective finite neutral, ionised, and total line counts; done means the returned uncertainties reflect those denominators and remain consistent with the other terms.

Written by the indexing model from the issue text.

Assessment

Tech stack
julia
Domain
backend
Issue type
Bug
Difficulty
3/5
Estimated time
1-2 days
Activity status
Active
Clarity
Clearly specified
Newbie friendliness
74/100

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.